There are many stochastic parameters to have effect on the inherent vibration of steam turbine blade. To meet with the requirement of blade high reliability operation, it is necessary to take these stochastic parameters into account, determine the statistical characteristics of inherent vibration frequencies, create the performance function to avoid the resonance and further carry out the vibration reliability analysis. A variable cross-section torsion blade is investigated and geometrical parameters, material parameters, installation parameters and rotation speed are considered as random variables, a method which is the combination of finite element method (FEM), response surface method (RSM) and Monte Carlo Simulation (MCS) is put forward to solve the blade vibration reliability analysis. Based on the finite element parametrical model of torsion blade and experiment design, the two kinds of response surface methods, multinomial response surface method (MRSM) and artificial neural network (ANN), are respectively applied to construct the approximation analytic expressions between the blade responses( including static frequencies and dynamic frequencies) and input random variables, which act as a surrogate of finite element solver to drastically reduce the number of solvers call. Then the surrogate is used for most of the samples needed in Monte Carlo method and the statistical parameters and cumulative distribution functions of the static frequencies and dynamic frequencies are obtained by Monte Carlo simulation. Aiming at the dangerous modes of vibration, the performance functions to avoid resonance are determined and the vibration reliability analysis is achieved. Moreover, the analysis results induced by the two different response surface methods MRSM and ANN are compared respectively to the result of Latin Hypercube sampling Monte Carlo simulation (LH-MCS), which is used as relative exact solution method.