The usual moment estimator of the harmonic mean is 1/mean(1/x_i), the reciprocal of the mean of the reciprocals of the x_i. This is a nonparametric estimator and is produced by the HARMEAN function in the DATA step. Limbrunner et. al. (2000) indicate that this moment estimator can be significantly biased especially for data with large coefficient of variation. They show that the maximum likelihood (ML) estimator "is nearly unbiased and [is] generally more efficient (lower RMSE) than the alternative estimators for lognormal observations."

The ML estimate, exp(μ - ½σ²), is easily produced by fitting a lognormal model in a procedure like PROC LIFEREG or PROC RELIABILITY and using the NLEstimate macro. It can also be done by fitting the model in PROC NLMIXED and using the ESTIMATE statement. In addition to a point estimate, a standard error, confidence interval, and test of the harmonic mean are also given as shown in the following example.

These statements use the HARMEAN function to compute and display the moment estimate of the harmonic mean of the values 1, 2, and 4.

      data a; 
        input y @@; 
        harmean = harmean(1, 2, 4);
        datalines;
      1 2 4
      ;
      
      proc print data=a(obs=1) noobs;
        var harmean;
        title "Moment estimator of Harmonic Mean";
        run;

The following statements use PROC LIFEREG to fit a lognormal model to the data. In the table of parameter estimates from PROC LIFEREG (not shown), the intercept parameter is an estimate of μ and the scale parameter is an estimate of σ. Since the ML estimator is a nonlinear function of the model parameters, it can be estimated using the NLEstimate macro. The NLEstimate macro uses the fitted model information (parameter estimates and covariance matrix) saved by the ODS OUTPUT statement in PROC LIFEREG. It then uses PROC NLMIXED to estimate the function and the delta method to obtain confidence limits. You write the function to be estimated using the parameter names and specify it in the f= macro parameter. A label can be provided in the label= parameter. See the description of the NLEstimate macro for details about displaying parameter names and using the macro.

      proc lifereg data=a; 
        model y=/dist=lognormal covb; 
        ods output parameterestimates=pe covb=covb;
      run;     
      %NLEstimate(inest=pe, incovb=covb, label=ML estimator of Harmonic Mean,
                  f=exp(b_p1-0.5*b_p2**2), df=3)

The estimated harmonic mean is similar to the moment estimate. A test and confidence interval for the harmonic mean are also provided.

PROC NLMIXED can also be used to fit the lognormal model to the data and estimate the harmonic mean. The variable LL is the log likelihood for the lognormal. The mean, MU, and variance, V, are estimated by PROC NLMIXED. The maximum likelihood estimator, exp(μ - ½σ²), is computed in the ESTIMATE statement.

      proc nlmixed data=a;
        ll=-0.5*(((log(y)-mu)**2)/v + log(v));
        model y~general(ll);
        estimate "ML estimator of Harmonic Mean" exp(mu-.5*v);
        run;

The results are the same as from the NLEstimate macro.

_____

Limbrunner J.F., Vogel R.M., and Brown L.C. (2000), "Estimation of Harmonic Mean of a Lognormal Variable", J. Hydrologic Engineering, 5(1), 59-66.

Operating System and Release Information