This example shows how to use some more advanced techniques with the Statistics and Machine Learning Toolbox™ function
mle to fit custom distributions to univariate data. The techniques include fitting models to censored data, and illustration of some of the numerical details
This example shows how to use some more advanced techniques with the Statistics and Machine Learning Toolbox™ function
mle to fit custom distributions to univariate data. The techniques include fitting models to censored data, and illustration of some of the numerical details of fitting with custom distributions.
See Fitting Custom Univariate Distributions for additional examples of fitting custom distributions to univariate data.
The extreme value distribution is often used to model failure times of mechanical parts, and experiments in such cases are sometimes only run for a fixed length of time. If not all of the experimental units have failed within that time, then the data are right-censored, that is, the value of some failure times are not known exactly, but only known to be larger than a certain value.
The Statistics and Machine Learning Toolbox includes the function
evfit, which fits an extreme value distribution to data, including data with censoring. However, for the purposes of this example, we will ignore
evfit, and show how to use
mle and custom distributions to fit a model to censored data, using the extreme value distribution.
Because the values for the censored data are not known exactly, maximum likelihood estimation becomes more difficult. In particular, both the PDF and the CDF are needed to compute the log-likelihood. Therefore, you must provide
mle with functions for both of those in order to fit censored data. The Statistics and Machine Learning Toolbox includes the functions
evcdf, so for this example, the work of writing the code has already been done.
We'll fit the model to simulated data. The first step is to generate some uncensored extreme value data.
rng(0,'twister'); n = 50; mu = 5; sigma = 2.5; x = evrnd(mu,sigma,n,1);
Next, we censor any values that are larger than a predetermined cutoff, by replacing those values with the cutoff value. This is known as Type II censoring, and with the cutoff at 7, about 12% of the original data end up being censored.
c = (x > 7); x(c) = 7; sum(c)/length(c)
ans = 0.1200
We can plot a histogram of these data, including a stacked bar to represent the censored observations.
[uncensCnts,binCtrs] = hist(x(~c)); censCnts = hist(x(c),binCtrs); bar(binCtrs,[uncensCnts' censCnts'],'stacked');
Although there is censoring, the fraction of censored observations is relatively small, and so the method of moments can provide a reasonable starting point for the parameter estimates. We compute the values of mu and sigma that correspond to the observed mean and standard deviation of the uncensored data.
sigma0 = sqrt(6)*std(x(~c))./pi
sigma0 = 2.3495
mu0 = mean(x(~c))-psi(1).*sigma0
mu0 = 3.5629
In addition to passing the data, x, and handles to the PDF and CDF functions into
mle, we also use the 'censoring' parameter to pass in the censoring vector, c. Because the scale parameter, sigma, must be positive, we specify lower parameter bounds.
mle returns the maximum likelihood estimates of the two extreme value distribution parameters, mu and sigma, as well as approximate 95% confidence intervals.
[paramEsts,paramCIs] = mle(x, 'censoring',c, 'pdf',@evpdf, 'cdf',@evcdf,...'start',[mu0 sigma0], 'lower',[-Inf,0])
paramEsts = 1×2 4.5530 3.0215
paramCIs = 2×2 3.6455 2.2937 5.4605 3.7494
Fitting a custom distribution requires an initial guess for the parameters, and it's often difficult to know how good or bad a starting point is a priori. In the previous example, if we had picked a starting point that was farther away from the maximum likelihood estimates, some of the observations could have been very far out in the tails of the extreme value distribution corresponding to the starting point. One of two things might then have happened.
First, one of the PDF values might have become so small that it underflowed to zero in double precision arithmetic. Second, one of the CDF values might have become so close to 1 that it rounded up in double precision. (It's also possible that a CDF value might have become so small as to underflow, but that turns out not to be a problem.)
Either of these conditions causes problems when
mle computes the log-likelihood, because they lead to log-likelihood values of -Inf, and the optimization algorithm in
mle cannot normally be expected to step out of such regions.
Knowing what the maximum likelihood estimates are, let's see what happens with a different starting point.
start = [1 1]; try [paramEsts,paramCIs] = mle(x, 'censoring',c, 'pdf',@evpdf, 'cdf',@evcdf,...'start',start, 'lower',[-Inf,0]) catch ME disp(ME.message) end
The CDF function returned values g