In the paper, we consider one of the most popular nonparametric estimators of regression reliability models proposed by R. Beran. Such estimator allows to evaluate the conditional reliability function with the given values of covariates by the following formula: (S) over cap (hn) (t vertical bar x) = Pi(Y(i)<= t) {1-W-n(i) (x,h(n))/1-Sigma W-i=1(j=1)n(j) (x,h(n))} where x is the value of covariate in the reliability function S(t vertical bar x); Y-( i) is the element of variational series; W-n(i) (x;h(n)), i=1,...,n, are the Nadaraya - Watson weights, i.e., W-n(i) (x; h(n)) = K (x-x(i)/h(n)) / Sigma(n)(j=1) K (x-x(j)/h(n)). It is well-known, the quality of the Beran estimator essentially depends on the chosen value of the bandwidth parameter h(n). In our previous paper, the method of selecting the optimal bandwidth parameter was proposed, which is based on the minimization of the distance of failure times with kernel estimation for the inverse reliability function. Here, we consider the modification of this method by solving such optimization problem: h(n)(opt) = arg min h(n) Sigma(n)(i=1) vertical bar (g) over cap((p) over cap (i) vertical bar x(i))-Y-i vertical bar, where (g) over cap((p) over cap vertical bar x(i)) = Sigma(n)(j=1) omega(j) ((p)over cap(i)).Y-j. The probabilities (p) over cap (i) are calculated by using the instrumentality of the Beran estimators, omega(j) ((p) over capi) are certain weights which can be calculated with various weight functions. We investigate the statistical properties of the Beran estimators by Monte Carlo simulations. It is shown that the accuracy of this estimators depend on the sample size, the number of covariates' values, the selection of the weight function's form, the method of smoothing parameter estimation and the type of kernel functions used in the smoothing parameter estimation and the bandwidth parameter calculation. The obtained results allow us to formulate recommendations for estimating the conditional reliability function by the Beran estimator. In our opinion, the most appropriate results are achieved by the Priestley - Chao weight function omega((2))(j) ((p) over cap (i)) = n ((p) over cap ((i)) - (p) over cap ((i-1))) K ((p) over cap (i)-(p) over cap (j)/b(NS)), with the smoothing parameter b(NS) = [8 pi R-1/2(K)/3 mu(2)(K)(2)n](1/5) (sigma) over cap, where mu(2) (K) = integral y(2)K(y)dy, R(K) = integral K-2(y)dy. We recommend the robust standard deviation estimator based on the mixing method using the median absolute deviation and the Hodges - Lehmann estimator: (sigma) over cap (robust) = 1.4826 med (i=1...n)vertical bar(p) over cap (i) - med (j=1...n,k=j+1...n)((p) over cap (j)+(p) over cap (k)/2)vertical bar. Also, we note that the quartic and Epanechnikov kernel functions lead to the most accurate Beran estimators.