Fit Zero-Truncated Poisson Distribution

Count data is often modeled using a Poisson distribution, and you can use the poissfit or fitdist function to fit a Poisson distribution. However, in some situations, counts that are zero are not recorded in the data, so fitting a Poisson distribution is not straightforward because of the missing zeros. In this case, fit a Poisson distribution to zero-truncated data by using the mle function and a custom distribution function.

First, generate some random Poisson data.

rng(18,'twister') % For reproducibility
lambda = 1.75;
n = 75;
x1 = poissrnd(lambda,n,1);

Next, remove all the zeros from the data to simulate the truncation.

x1 = x1(x1 > 0);

Check the number of samples in x1 after truncation.

length(x1)

ans = 65

Plot a histogram of the simulated data.

histogram(x1,0:1:max(x1)+1)

The data looks like a Poisson distribution except it contains no zeros. You can use a custom distribution that is identical to a Poisson distribution on the positive integers, but has no probability at zero. By using a custom distribution, you can estimate the Poisson parameter lambda while accounting for the missing zeros.

You need to define the zero-truncated Poisson distribution by its probability mass function (pmf). Create an anonymous function to compute the probability for each point in x1, given a value for the Poisson distribution's mean parameter lambda. The pmf for a zero-truncated Poisson distribution is the Poisson pmf normalized so that it sums to one. With zero truncation, the normalization is 1–Probability(x1<0).

pf_truncpoiss = @(x1,lambda) poisspdf(x1,lambda)./(1-poisscdf(0,lambda));

For simplicity, assume that all the x1 values given to this function are positive integers, with no checks. For error checking or a more complicated distribution that takes more than a single line of code, you must define the function in a separate file.

Find a reasonable rough first guess for the parameter lambda. In this case, use the sample mean.

start = mean(x1)

start = 2.2154

Provide mle with the data, custom pmf function, initial parameter value, and lower bound of the parameter. Because the mean parameter of the Poisson distribution must be positive, you also need to specify a lower bound for lambda. The mle function returns the maximum likelihood estimate of lambda, and optionally, the approximate 95% confidence intervals for the parameters.

[lambdaHat,lambdaCI] = mle(x1,'pdf',pf_truncpoiss,'Start',start, ...
    'LowerBound',0)

lambdaHat = 1.8760

lambdaCI = 2×1

    1.4990
    2.2530

The parameter estimate is smaller than the sample mean. The maximum likelihood estimate accounts for the zeros not present in the data.

Alternatively, you can specify the truncation bounds by using the TruncationBounds name-value argument.

[lambdaHat2,lambdaCI2] = mle(x1,'Distribution','Poisson', ...
    'TruncationBounds',[0 Inf])

lambdaHat2 = 1.8760

lambdaCI2 = 2×1

    1.4990
    2.2530

You can also compute a standard error estimate for lambda by using the large-sample variance approximation returned by mlecov.

avar = mlecov(lambdaHat,x1,'pdf',pf_truncpoiss);
stderr = sqrt(avar)

stderr = 0.1923

To visually check the fit, plot the fitted pmf against a normalized histogram of the raw data.

histogram(x1,'Normalization','pdf')
xgrid = min(x1):max(x1);
pmfgrid = pf_truncpoiss(xgrid,lambdaHat);
hold on
plot(xgrid,pmfgrid,'-')
xlabel('x1')
ylabel('Probability')
legend('Sample Data','Fitted pmf','Location','best')
hold off

Fit Upper-Truncated Normal Distribution

Continuous data can sometimes be truncated. For example, observations larger than some fixed value might not be recorded because of limitations in data collection.

In this case, simulate data from a truncated normal distribution. First, generate some random normal data.

n = 500;
mu = 1;
sigma = 3;
rng('default') % For reproducibility
x2 = normrnd(mu,sigma,n,1);

Next, remove any observations that fall beyond the truncation point xTrunc. Assume that xTrunc is a known value that you do not need to estimate.

xTrunc = 4;
x2 = x2(x2 < xTrunc);

Check the number of samples in x2 after truncation.

length(x2)

ans = 430

Create a histogram of the simulated data.

histogram(x2)

Fit the simulated data with a custom distribution that is identical to a normal distribution for x2 < xTrunc, but has zero probability above xTrunc. By using a custom distribution, you can estimate the normal parameters mu and sigma while accounting for the missing tail.

Define the truncated normal distribution by its pdf. Create an anonymous function to compute the probability density value for each point in x, given values for the parameters mu and sigma. With the truncation point fixed and known, the pdf for a truncated normal distribution is the pdf truncated and then normalized so that it integrates to one. The normalization is the cdf evaluated at xTrunc. For simplicity, assume that all x2 values are less than xTrunc, without checking.

pdf_truncnorm = @(x2,mu,sigma) ...
    normpdf(x2,mu,sigma)./normcdf(xTrunc,mu,sigma);

Because you do not need to estimate the truncation point xTrunc, it is not included with the input distribution parameters of the custom pdf function. xTrunc is also not part of the data vector input argument. An anonymous function can access variables in the workspace, so you do not have to pass xTrunc to the anonymous function as an additional argument.

Provide a rough starting guess for the parameter estimates. In this case, because the truncation is not extreme, use the sample mean and standard deviation.

start = [mean(x2),std(x2)]

start = 1×2

    0.1585    2.4125

Provide mle with the data, custom pdf function, initial parameter values, and lower bounds of the parameters. Because sigma must be positive, you also need to specify lower parameter bounds. mle returns the maximum likelihood estimates of mu and sigma as a single vector, as well as a matrix of approximate 95% confidence intervals for the two parameters.

[paramEsts,paramCIs] = mle(x2,'pdf',pdf_truncnorm,'Start',start, ...
    'LowerBound',[-Inf 0])

paramEsts = 1×2

    1.1298    3.0884

paramCIs = 2×2

    0.5713    2.7160
    1.6882    3.4607

The estimates of mu and sigma are larger than the sample mean and standard deviation. The model fit accounts for the missing upper tail of the distribution.

Alternatively, you can specify the truncation bounds by using the TruncationBounds name-value argument.

[paramEsts2,paramCIs2] = mle(x2,'Distribution','Normal', ...
    'TruncationBounds',[-Inf xTrunc])

paramEsts2 = 1×2

    1.1297    3.0884

paramCIs2 = 2×2

    0.5713    2.7160
    1.6882    3.4607

You can compute an approximate covariance matrix for the parameter estimates using mlecov. The approximation typically works well for large samples, and you can approximate the standard errors by the square roots of the diagonal elements.

acov = mlecov(paramEsts,x2,'pdf',pdf_truncnorm)

acov = 2×2

    0.0812    0.0402
    0.0402    0.0361

stderr = sqrt(diag(acov))

stderr = 2×1

    0.2849
    0.1900

To visually check the fit, plot the fitted pdf against a normalized histogram of the raw data.

histogram(x2,'Normalization','pdf')
xgrid = linspace(min(x2),max(x2));
pdfgrid = pdf_truncnorm(xgrid,paramEsts(1),paramEsts(2));
hold on
plot(xgrid,pdfgrid,'-')
xlabel('x2')
ylabel('Probability Density')
legend('Sample Data','Fitted pdf','Location','best')
hold off

Fit Mixture of Two Normal Distributions

Some data sets exhibit bimodality, or even multimodality, and fitting a standard distribution to such data is usually not appropriate. However, a mixture of simple unimodal distributions can often model such data very well.

In this case, fit a mixture of two normal distributions to simulated data. Consider simulated data with the following constructive definition:

Generate a data set from a mixture of Student's t distributions instead of using the same model that you are fitting. By using different distributions, similar to a technique used in a Monte-Carlo simulation, you can test how robust a fitting method is to departures from the assumptions of the model being fit.

rng(10) % For reproducibility
x3 = [trnd(20,1,50) trnd(4,1,100)+3];
histogram(x3)

Define the model to fit by creating an anonymous function that computes the probability density. The pdf for a mixture of two normal distributions is a weighted sum of the pdfs of the two normal components, weighted by the mixture probability. The anonymous function takes six inputs: a vector of data at which to evaluate the pdf and five distribution parameters. Each component has parameters for its mean and standard deviation.

pdf_normmixture = @(x3,p,mu1,mu2,sigma1,sigma2) ...
    p*normpdf(x3,mu1,sigma1) + (1-p)*normpdf(x3,mu2,sigma2);

You also need an initial guess for the parameters. Defining a starting point becomes more important as the number of model parameters increases. Here, start with an equal mixture (p = 0.5) of normal distributions, centered at the two quartiles of the data, with equal standard deviations. The starting value for the standard deviation comes from the formula for the variance of a mixture in terms of the mean and variance of each component.

pStart = .5;
muStart = quantile(x3,[.25 .75])

muStart = 1×2

    0.3351    3.3046

sigmaStart = sqrt(var(x3) - .25*diff(muStart).^2)

sigmaStart = 1.1602

start = [pStart muStart sigmaStart sigmaStart];

Specify bounds of zero and one for the mixing probability, and lower bounds of zero for the standard deviations. Set the remaining elements of the bounds vectors to +Inf or –Inf, to indicate no restrictions.

lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
    'LowerBound',lb,'UpperBound',ub)

Warning: Maximum likelihood estimation did not converge.  Iteration limit exceeded.

paramEsts = 1×5

    0.3273   -0.2263    2.9914    0.9067    1.2059

The warning message indicates that the function does not converge with the default iteration settings. Display the default options.

statset('mlecustom')

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

The default maximum number of iterations for custom distributions is 200. Override the default to increase the number of iterations, using an options structure created with the statset function. Also, increase the maximum function evaluations.

options = statset('MaxIter',300,'MaxFunEvals',600);
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
    'LowerBound',lb,'UpperBound',ub,'Options',options)

paramEsts = 1×5

    0.3273   -0.2263    2.9914    0.9067    1.2059

The final iterations to convergence are significant only in the last few digits of the result. However, a best practice is to always make sure that convergence is reached.

To visually check the fit, plot the fitted density against a probability histogram of the raw data.

histogram(x3,'Normalization','pdf')
hold on
xgrid = linspace(1.1*min(x3),1.1*max(x3),200);
pdfgrid = pdf_normmixture(xgrid, ...
    paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5));
plot(xgrid,pdfgrid,'-')
hold off
xlabel('x3')
ylabel('Probability Density')
legend('Sample Data','Fitted pdf','Location','best')

Alternatively, for a mixture of normal distributions, you can use the fitgmdist function. The estimates can be different due to initial estimates and settings of the iterative algorithm.

Mdl = fitgmdist(x3',2)

Mdl = 

Gaussian mixture distribution with 2 components in 1 dimensions
Component 1:
Mixing proportion: 0.329180
Mean:   -0.2200

Component 2:
Mixing proportion: 0.670820
Mean:    2.9975

Mdl.Sigma

ans = 
ans(:,:,1) =

    0.8274


ans(:,:,2) =

    1.4437

Fit Weighted Normal Distribution to Data with Unequal Precisions

Assume that you have 10 data points, where each point is actually the average of anywhere from one to eight observations. The original observations are not available, but the number of observations for each data point is known. The precision of each point depends on its corresponding number of observations. You need to estimate the mean and standard deviation of the raw data.

x4 = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]';
m = [8 2 1 3 8 4 2 5 2 4]';

The variance of each data point is inversely proportional to its corresponding number of observations, so use 1/m to weight the variance of each data point in a maximum likelihood fit.

w = 1./m;

In this model, you can define the distribution by its pdf. However, using a logarithm of pdf is more suitable, because the normal pdf has the form

  c .* exp(-0.5 .* z.^2),

and mle takes the log of the pdf to compute the loglikelihood. So, instead, create a function that computes the logarithm of pdf directly.

The logarithm of pdf function must compute the logarithm of the probability density for each point in x, given normal distribution parameters mu and sigma. It also needs to account for the different variance weights. Define a function named helper_logpdf_wn1 in a separate file helper_logpdf_wn1.m.

function logy = helper_logpdf_wn1(x,m,mu,sigma)
%HELPER_LOGPDF_WN1 Logarithm of pdf for a weight normal distribution
% This function supports only the example Fit Custom Distributions 
% (customdist1demo.mlx) and might change in a future release.
v = sigma.^2 ./ m;
logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
end

Provide a rough first guess for the parameter estimates. In this case, use the unweighted sample mean and standard deviation.

start = [mean(x4),std(x4)]

start = 1×2

    1.0280    1.5490

Because sigma must be positive, you need to specify lower parameter bounds.

[paramEsts1,paramCIs1] = mle(x4,'logpdf', ...
    @(x,mu,sigma)helper_logpdf_wn1(x,m,mu,sigma), ...
    'Start',start,'LowerBound',[-Inf,0])

paramEsts1 = 1×2

    0.6244    2.8823

paramCIs1 = 2×2

   -0.2802    1.6191
    1.5290    4.1456

The estimate of mu is less than two-thirds of the estimate of the sample mean. The estimate is influenced by the most reliable data points, that is, the points based on the largest number of raw observations. In this data set, those points tend to pull the estimate down from the unweighted sample mean.

Fit Normal Distribution Using Parameter Transformation

The mle function computes confidence intervals for the parameters using a large-sample normal approximation for the distribution of the estimators if an exact method is not available. For small sample sizes, you can improve the normal approximation by transforming one or more parameters. In this case, transform the scale parameter of a normal distribution to its logarithm.

First, define a new log pdf function named helper_logpdf_wn2 that uses a transformed parameter for sigma.

function logy = helper_logpdf_wn2(x,m,mu,logsigma)
%HELPER_LOGPDF_WN2 Logarithm of pdf for a weight normal distribution with 
% log(sigma) parameterization
% This function supports only the example Fit Custom Distributions 
% (customdist1demo.mlx) and might change in a future release.
v = exp(logsigma).^2 ./ m;
logy = -(x-mu).^2 ./ (2.*v) - .5.*log(2.*pi.*v);
end

Use the same starting point transformed to the new parameterization for sigma, that is, the log of the sample standard deviation.

start = [mean(x4),log(std(x4))]

start = 1×2

    1.0280    0.4376

Because sigma can be any positive value, log(sigma) is unbounded, and you do not need to specify lower or upper bounds.

[paramEsts2,paramCIs2] = mle(x4,'logpdf', ...
    @(x,mu,sigma)helper_logpdf_wn2(x,m,mu,sigma), ...
    'Start',start)

paramEsts2 = 1×2

    0.6244    1.0586

paramCIs2 = 2×2

   -0.2802    0.6203
    1.5290    1.4969

Because the parameterization uses log(sigma), you have to transform back to the original scale to get an estimate and confidence interval for sigma.

sigmaHat = exp(paramEsts2(2))

sigmaHat = 2.8823

sigmaCI = exp(paramCIs2(:,2))

sigmaCI = 2×1

    1.8596
    4.4677

The estimates for both mu and sigma are the same as in the first fit, because maximum likelihood estimates are invariant to parameterization. The confidence interval for sigma is slightly different from paramCIs1(:,2).