normfit

Generate 1000 normal random numbers from the normal distribution with mean 3 and standard deviation 5.

rng('default') % For reproducibility
x = normrnd(3,5,[1000,1]);

Find the parameter estimates and the 99% confidence intervals.

[muHat,sigmaHat,muCI,sigmaCI] = normfit(x,0.01)

muHat = 
2.8368

sigmaHat = 
4.9948

muCI = 2×1

    2.4292
    3.2445

sigmaCI = 2×1

    4.7218
    5.2989

muHat is the sample mean, and sigmaHat is the square root of the unbiased estimator of the variance. muCI and sigmaCI contain the 99% confidence intervals of the mean and standard deviation parameters, respectively. The first row is the lower bound, and the second row is the upper bound.

Find the MLEs of a data set with censoring by using normfit. Use statset to specify the iterative algorithm options that normfit uses to compute MLEs for censored data, and then find the MLEs again.

Load the sample data.

load lightbulb

The first column of the data contains the lifetime (in hours) of two types of bulbs. The second column contains the binary variable indicating whether the bulb is fluorescent or incandescent. 1 indicates that the bulb is fluorescent, and 0 indicates that the bulb is incandescent. The third column contains the censorship information, where 0 indicates the bulb is observed until failure, and 1 indicates the item (bulb) is censored.

Find the indices for fluorescent bulbs.

idx = find(lightbulb(:,2) == 0);

Assume that the lifetime follows the normal distribution, and find the MLEs of the normal distribution parameters. The second input argument of normfit specifies the confidence level. Pass in [] to use its default value 0.05. The third input argument specifies the censorship information.

censoring = lightbulb(idx,3) == 1;
[muHat1,sigmaHat1] = normfit(lightbulb(idx,1),[],censoring)

muHat1 = 
9.4966e+03

sigmaHat1 = 
3.0640e+03

Display the default algorithm parameters that normfit uses to estimate the normal distribution parameters.

statset('normfit')

ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 200
          MaxIter: 100
           TolBnd: 1.0000e-06
           TolFun: 1.0000e-08
       TolTypeFun: []
             TolX: 1.0000e-08
         TolTypeX: []
          GradObj: []
         Jacobian: []
        DerivStep: []
      FunValCheck: []
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

Save the options using a different name. Change how the results are displayed (Display) and the termination tolerance for the objective function (TolFun).

options = statset('normfit');
options.Display = 'final';
options.TolFun = 1e-10;

Alternatively, you can specify algorithm parameters by using the name-value pair arguments of the function statset.

options = statset('Display','final','TolFun',1e-10);

Find the MLEs with the new algorithm parameters.

[muHat2,sigmaHat2] = normfit(lightbulb(idx,1),[],censoring,[],options)

Successful convergence: Norm of gradient less than OPTIONS.TolFun

muHat2 = 
9.4966e+03

sigmaHat2 = 
3.0640e+03

normfit displays a report on the final iteration.

The function normfit finds the sample mean and the square root of the unbiased estimator of the variance with no censoring. The sample mean is equal to the MLE of the mean parameter, but the square root of the unbiased estimator of the variance is not equal to the MLE of the standard deviation parameter.

Find the normal distribution parameters by using normfit, convert them into MLEs, and then compare the negative log likelihoods of the estimates by using normlike.

Generate 100 normal random numbers from the standard normal distribution.

rng('default') % For reproducibility
n = 100;
x = normrnd(0,1,[n,1]);

Find the sample mean and the square root of the unbiased estimator of the variance.

[muHat,sigmaHat] = normfit(x)

muHat = 
0.1231

sigmaHat = 
1.1624

Convert the square root of the unbiased estimator of the variance into the MLE of the standard deviation parameter.

sigmaHat_MLE = sqrt((n-1)/n)*sigmaHat

sigmaHat_MLE = 
1.1566

The difference between sigmaHat and sigmaHat_MLE is negligible for large n.

Alternatively, you can find the MLEs by using the function mle.

phat = mle(x)

phat = 1×2

    0.1231    1.1566

phat(1) and phat(2) are the MLEs of the mean and the standard deviation parameter, respectively.

Confirm that the log likelihood of the MLEs (muHat and sigmaHat_MLE) is greater than the log likelihood of the unbiased estimators (muHat and sigmaHat) by using the normlike function.

logL = -normlike([muHat,sigmaHat],x)

logL = 
-156.4424

logL_MLE = -normlike([muHat,sigmaHat_MLE],x)

logL_MLE = 
-156.4399