mlecov

Asymptotic covariance of maximum likelihood estimators

Syntax

acov = mlecov(params,data,'pdf',pdf)

acov = mlecov(params,data,'pdf',pdf,'cdf',cdf)

acov = mlecov(params,data,'logpdf',logpdf)

acov = mlecov(params,data,'logpdf',logpdf,'logsf',logsf)

acov = mlecov(params,data,'nloglf',nloglf)

acov = mlecov(___,Name,Value)

Description

acov = mlecov(params,data,'pdf',pdf) returns an approximation to the asymptotic covariance matrix of the maximum likelihood estimators of the parameters for a distribution specified by the custom probability density function pdf.

mlecov computes a finite difference approximation to the Hessian of the log-likelihood at the maximum likelihood estimates params, given the observed data, and returns the negative inverse of that Hessian.

acov = mlecov(params,data,'pdf',pdf,'cdf',cdf) returns acov for a distribution specified by the custom probability density function pdf and cumulative distribution function cdf.

acov = mlecov(params,data,'logpdf',logpdf) returns acov for a distribution specified by the custom log probability density function logpdf.

acov = mlecov(params,data,'logpdf',logpdf,'logsf',logsf) returns acov for a distribution specified by the custom log probability density function logpdf and custom log survival function logsf.

acov = mlecov(params,data,'nloglf',nloglf) returns acov for a distribution specified by the custom negative loglikelihood function nloglf.

acov = mlecov(___,Name,Value) specifies options using name-value pair arguments in addition to any of the input arguments in previous syntaxes. For example, you can specify the censored data and frequency of observations.

Examples

Custom Probability Density Function

Load the sample data.

load carbig

The vector Weight shows the weights of 406 cars.

In the MATLAB^® Editor, create a function that returns the probability density function (pdf) of a lognormal distribution. Save the file in your current working folder as lognormpdf.m.

function newpdf = lognormpdf(data,mu,sigma)
newpdf = exp((-(log(data)-mu).^2)/(2*sigma^2))./(data*sigma*sqrt(2*pi));

Estimate the parameters, mu and sigma, of the custom-defined distribution.

phat = mle(Weight,'pdf',@lognormpdf,'start',[4.5 0.3])

phat =

    7.9600    0.2804

Compute the approximate covariance matrix of the parameter estimates.

acov = mlecov(phat,Weight,'pdf',@lognormpdf)

acov =

   1.0e-03 *

    0.1937   -0.0000
   -0.0000    0.0968

Estimate the standard errors of estimates.

se = sqrt(diag(acov))

se =

    0.0139
    0.0098

The standard error of the estimates of mu and sigma are 0.0139 and 0.0098, respectively.

Custom Log Probability Density Function

In the MATLAB Editor, create a function that returns the log probability density function of a beta distribution. Save the file in your current working folder as betalogpdf.m.

function logpdf = betalogpdf(x,a,b)
logpdf = (a-1)*log(x)+(b-1)*log(1-x)-betaln(a,b);

Generate sample data from a beta distribution with parameters 1.23 and 3.45 and estimate the parameters using the simulated data.

rng('default')
x = betarnd(1.23,3.45,25,1);
phat = mle(x,'dist','beta')

phat =

    1.1213    2.7182

Compute the approximate covariance matrix of the parameter estimates.

acov = mlecov(phat,x,'logpdf',@betalogpdf)

acov =

    0.0810    0.1646
    0.1646    0.6074

Custom Log pdf and Survival Function

Load the sample data.

load(fullfile(matlabroot,'examples','stats','readmissiontimes.mat'));

The sample data includes ReadmissionTime, which has readmission times for 100 patients. The column vector Censored has the censorship information for each patient, where 1 indicates a censored observation, and 0 indicates the exact readmission time is observed. This is simulated data.

Define a custom log probability density and survival function.

custlogpdf = @(data,lambda,k) log(k)-k*log(lambda)...
             +(k-1)*log(data)-(data/lambda).^k;
custlogsf = @(data,lambda,k) -(data/lambda).^k;

Estimate the parameters, lambda and k, of the custom distribution for the censored sample data.

phat = mle(ReadmissionTime,'logpdf',custlogpdf,...
'logsf',custlogsf,'start',[1,0.75],'Censoring',Censored)

phat = 1×2

    9.2090    1.4223

The scale and shape parameters of the custom-defined distribution are 9.2090 and 1.4223, respectively.

Compute the approximate covariance matrix of the parameter estimates.

acov = mlecov(phat,ReadmissionTime,...
'logpdf',custlogpdf,'logsf',custlogsf,'Censoring',Censored)

acov = 2×2

    0.5653    0.0102
    0.0102    0.0163

Custom Log Negative Likelihood Function

Load the sample data.

load(fullfile(matlabroot,'examples','stats','readmissiontimes.mat'));

The sample data includes ReadmissionTime, which has readmission times for 100 patients. This is simulated data.

Define a negative log likelihood function.

custnloglf = @(lambda,data,cens,freq) -length(data)*log(lambda)...
+ nansum(lambda*data);

Estimate the parameters of the defined distribution.

phat = mle(ReadmissionTime,'nloglf',custnloglf,'start',0.05)

phat = 0.1462

Compute the variance of the parameter estimate.

acov = mlecov(phat,ReadmissionTime,'nloglf',custnloglf)

acov = 2.1374e-04

Compute the standard error.

sqrt(acov)

ans = 0.0146

Input Arguments

params — Parameter estimates

scalar value | vector

Parameter estimates, specified as a scalar value or vector of scalar values. These parameter estimates must be maximum likelihood estimates. For example, you can specify parameter estimates returned by mle.

Data Types: single | double

data — Sample data

vector

Sample data mle uses to estimate the distribution parameters, specified as a vector.

Data Types: single | double

pdf — Custom probability density function

function handle

Custom probability distribution function, specified as a function handle created using @.

This custom function accepts the vector data and one or more individual distribution parameters as input parameters, and returns a vector of probability density values.

For example, if the name of the custom probability density function is newpdf, then you can specify the function handle in mlecov as follows.

Example: @newpdf

Data Types: function_handle

cdf — Custom cumulative distribution function

function handle

Custom cumulative distribution function, specified as a function handle created using @.

This custom function accepts the vector data and one or more individual distribution parameters as input parameters, and returns a vector of cumulative probability values.

You must define cdf with pdf if data is censored and you use the 'Censoring' name-value pair argument. If 'Censoring' is not present, you do not have to specify cdf while using pdf.

For example, if the name of the custom cumulative distribution function is newcdf, then you can specify the function handle in mlecov as follows.

Example: @newcdf

Data Types: function_handle

logpdf — Custom log probability density function

function handle

Custom log probability density function, specified as a function handle created using @.

This custom function accepts the vector data and one or more individual distribution parameters as input parameters, and returns a vector of log probability values.

For example, if the name of the custom log probability density function is customlogpdf, then you can specify the function handle in mlecov as follows.

Example: @customlogpdf

Data Types: function_handle

logsf — Custom log survival function

function handle

Custom log survival function, specified as a function handle created using @.

This custom function accepts the vector data and one or more individual distribution parameters as input parameters, and returns a vector of log survival probability values.

You must define logsf with logpdf if data is censored and you use the 'Censoring' name-value pair argument. If 'Censoring' is not present, you do not have to specify logsf while using logpdf.

For example, if the name of the custom log survival function is logsurvival, then you can specify the function handle in mlecov as follows.

Example: @logsurvival

Data Types: function_handle

nloglf — Custom negative loglikelihood function

function handle

Custom negative loglikelihood function, specified as a function handle created using @.

This custom function accepts the following input arguments.

params	Vector of distribution parameter values
data	Vector of data
cens	Boolean vector of censored values
freq	Vector of integer data frequencies

nloglf must accept all four arguments even if you do not use the 'Censoring' or 'Frequency' name-value pair arguments. You can write 'nloglf' to ignore cens and freq arguments in that case.

nloglf returns a scalar negative loglikelihood value and optionally, a negative loglikelihood gradient vector (see the 'GradObj' field in 'Options').

If the name of the custom negative log likelihood function is negloglik, then you can specify the function handle in mlecov as follows.

Example: @negloglik

Data Types: function_handle

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'Censoring',cens,'Options',opt specifies that mlecov reads the censored data information from the vector cens and performs according to the new options structure opt.

'Censoring' — Indicator for censoring

array of 0s (default) | array of 0s and 1s

Indicator for censoring, specified as the comma-separated pair consisting of 'Censoring' and a Boolean array of the same size as data. Use 1 for observations that are right censored and 0 for observations that are fully observed. The default is all observations are fully observed.

For censored data, you must use cdf with pdf, or logsf with logpdf, or nloglf must be defined to account for censoring.

For example, if the censored data information is in the binary array called Censored, then you can specify the censored data as follows.

Example: 'Censoring',Censored

Data Types: logical

'Frequency' — Frequency of observations

array of 1s (default) | vector of nonnegative integer counts

Frequency of observations, specified as the comma-separated pair consisting of 'Frequency' and an array containing nonnegative integer counts, which is the same size as data. The default is one observation per element of data.

For example, if the observation frequencies are stored in an array named Freq, you can specify the frequencies as follows.

Example: 'Frequency',Freq

Data Types: single | double

'Options' — Numerical options

structure

Numerical options for the finite difference Hessian calculation, specified as the comma-separated pair consisting of 'Options' and a structure returned by statset.

You can set the options under a new name and use it in the name-value pair argument. The applicable statset parameters are as follows.

Parameter	Value
'GradObj'	Default is 'off'. 'on' or 'off', indicating whether or not the function provided with the nloglf input argument can return the gradient vector of the negative log-likelihood as a second output.
'DerivStep'	Default is eps^(1/4). Relative step size used in finite difference for Hessian calculations. It can be a scalar, or the same size as params. A smaller value than the default might be appropriate if 'GradObj' is 'on'.

Example: 'Options',statset('mlecov')

Data Types: struct

Output Arguments

acov — Approximation to asymptotic covariance matrix

p-by-p matrix

Approximation to asymptotic covariance matrix, returned as a p-by-p matrix, where p is the number of parameters in params.

More About

Survival Function

The survival function is the probability of survival as a function of time. It is also called the survivor function. It gives the probability that the survival time of an individual exceeds a certain value. Since the cumulative distribution function, F(t), is the probability that the survival time is less than or equal to a given point in time, the survival function for a continuous distribution, S(t), is the complement of the cumulative distribution function: S(t) = 1 – F(t).