拟合自定义分布:一个零截断泊松的例子

计数数据使用的是泊松分布常为蓝本，并且可以使用统计和机器学习工具箱功能poissfit来拟合泊松模型。然而，在某些情况下，为零的计数不会被记录到数据中，因此，由于那些“丢失的零”，拟合泊松分布并不是一件简单的事情。本例将展示如何使用该函数将泊松分布适合于零截断的数据MLE。

在这个例子中，我们会从零截断泊松分布使用模拟数据。首先，我们生成一些随机数据泊松。

rng (18,“扭腰”);n = 75;λ= 1.75;x = poissrnd(λ,n, 1);

接下来，我们从数据中删除所有的零来模拟截断。

x = x(x > 0);长度（x）的

ans = 65

这是这些模拟数据的柱状图。请注意，数据看似合理像泊松分布，但没有零。我们将适合他们分布是相同的正整数泊松，但有在零没有可能性。通过这种方式，我们可以在占“零缺失”估计泊松参数拉姆达。

嘘(x,[0:1:马克斯(x) + 1]);

第一步是定义零截断泊松分布的概率函数(PF)。我们将创建一个函数来计算x中每个点的概率，给定泊松分布的均值参数。零截尾的泊松函数的PF就是通常的泊松函数，经过重新正规化后，其和为1。在截断为0的情况下，重正化过程仅为1-Pr{0}。为PF创建函数的最简单方法是使用匿名函数。

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

为简单起见，我们假设所有给这个函数的x值将是正整数，不检查。错误检查，或更复杂的分布，可能会需要超过一个单一的代码行，这表明函数应在一个单独的文件中定义。

下一步是提供参数拉姆达合理粗糙的第一个猜测。在这种情况下，我们只用样本平均值。

开始=意味着(x)

开始= 2.2154

我们提供MLE对于数据和匿名函数，使用'pdf'参数。(泊松是离散的，所以这是一个概率函数，而不是PDF)因为泊松分布的平均参数必须是正的，所以我们也为指定一个下界。MLE返回lambda的最大似然估计值，以及参数的大约95%置信区间(可选)。

[lambdaHat，lambdaCI] = MLE（X，“pdf”，pf_truncpoiss，“开始”开始,“低”,0)

lambdaHat = 1.8760

lambdaCI =2×11.4990 2.2530

注意，参数估计值小于样本均值。这是应该的，因为最大似然估计解释了数据中不存在的缺失零。

我们还可以计算lambda的标准误差估计，使用返回的大样本方差近似mlecov。

mlecov(lambdaHat, x，“pdf”,pf_truncpoiss);stderr =√阿瓦尔人

stderr = 0.1923

供给其他值的分布函数：截短的普通实施例

它有时也碰巧连续数据被截断。例如，比一些固定值较大的意见可能不会因为这样的数据被收集的局限性记录。这个例子将展示如何适应一个正态分布截断数据，使用功能MLE。

在这个例子中，我们模拟从截断正态分布的数据。首先，我们生成一些随机的正常数据。

n = 75;亩= 1;西格玛= 3;X = normrnd（亩，SIGMA，N，1）;

接下来，我们会删除落在超出截断点，xTrunc任何意见。在本例中，我们假设xTrunc是已知的，并且不需要进行估计。

xTrunc = 4;X = X（X 
              

ans = 64

这是这些模拟数据的柱状图。我们将用一个与x < xTrunc的正常情况相同的分布来匹配它们，但在xTrunc之上的概率为零。这样，我们可以估计正常参数mu和sigma，同时考虑“缺尾”。

嘘(x, (-10: .5:4));

如前面的例子中，我们将通过它的PDF定义截断正态分布，并且创建一个函数来计算用于以x的每个点的参数mu和sigma的概率密度，给定的值。随着截断点固定的和已知的，对于截断正常PDF只是平时的正常PDF，截断，然后重新归一化，使其集成为一个。重整化只是在xTrunc评估的CDF。为简单起见，我们假设所有的x值都小于xTrunc，不检查。我们使用一个匿名函数来定义PDF。

pdf_truncnorm = @（X，μ，西格马）normpdf（X，μ，西格马）./ normcdf（xTrunc，μ，西格马）;

截断点xTrunc没有得到估计，因此不在PDF函数的输入参数列表中的分布参数中。xTrunc也不是数据向量输入参数的一部分。使用匿名函数，我们可以简单地引用工作区中已经定义的变量xTrunc，而无需担心将其作为附加参数传入。

我们还需要为参数估计提供一个粗略的初始猜测。在这种情况下，由于截断不是太极端，样本均值和标准差可能会工作得很好。

开始= [(x),性病(x))

开始=1×20.2735 2.2660

我们提供MLE对于数据和匿名函数，使用'pdf'参数。由于西格玛必须是积极的，我们也指定较低参数范围。MLE返回mu和sigma作为单个向量的最大似然估计值，以及两个参数的约95%置信区间的矩阵。

[paramEsts, paramCIs] =大中型企业(x,“pdf”，pdf_truncnorm，“开始”开始,“低”(从0))

paramEsts =1×20.9911 - 2.7800

paramCIs =2×2-0.1605 1.9680 2.1427 3.5921

注意，mu和sigma的估计是比样品平均值和标准偏差较大的相当多。这是因为该模型拟合已占到分布的“失踪”上尾。

我们可以为参数估计计算一个近似的协方差矩阵mlecov。在大样本中，近似通常是合理的，估计的标准误差可以用对角线元素的平方根近似。

ACOV = mlecov（paramEsts，X，“pdf”pdf_truncnorm)

acov =2×20.3452 0.1660 0.1660 0.1717

标准错误= SQRT（DIAG（ACOV））

stderr =2×10.5876 - 0.4143

装修更复杂的分布：两条法线的混合物

一些数据集显示出双峰性，甚至是多模性，对这样的数据进行标准分布的拟合通常是不合适的。然而，混合使用简单的单峰分布通常可以很好地对这些数据进行建模。事实上，甚至可以基于特定于应用程序的知识对混合中的每个组件的源给出解释。

在本例中，我们将把两个正态分布的混合物拟合到一些模拟数据中。这种混合可以用以下产生随机值的构造定义来描述:

首先，抛一枚有偏差的硬币。如果是正面，从均值为mu_1、标准差为sigma_1的正态分布中随机选择一个值。如果硬币反面朝上，从均值为mu_2、标准差为sigma_2的正态分布中随机选择一个值。

在本例中，我们将从学生的t分布混合生成数据，而不是使用与我们拟合的相同的模型。这是你在蒙特卡罗模拟中可能会做的事情，来测试一个拟合方法在偏离模型拟合假设的情况下有多健壮。然而，在这里，我们只拟合一个模拟数据集。

X = [TRND（20,1,50）TRND（4,1,100）3];HIST（X，-2.25：0.5：7.25）;

正如在前面的例子中，我们将通过创建一个计算概率密度函数定义模型来拟合。为两条法线的混合物PDF是两个正常组分的PDF文件，通过将混合物概率加权的只是一个加权和。这PDF是很简单的创建使用匿名函数。该函数有六个输入：在该评估的PDF数据的矢量和分布的五个参数。每个组件具有用于它的平均值和标准偏差参数;该混合物概率使得总共五个。

pdf_normmixture = @（X，P，MU1，MU2，sigma1，σ-2）…P * normpdf（X，MU1，sigma1）+（1-P）* normpdf（X，MU2，σ-2）;

我们还需要对这些参数的初始猜测。在多个参数的模型，更合理的出发点很重要。对于这个例子，我们将与法线的等量混合物（P = 0.5），在数据的两个25为中心，以相等的标准偏差开始。对于标准偏差的起始值来自公式混合物的方差的平均值和各成分的方差的条款。

pStart = 5;muStart =分位数(x, [。25。)

muStart =1×20.5970 - 3.2456

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

sigmaStart = 1.2153

开始= [PSTART muStart sigmaStart sigmaStart];

最后，我们需要指定混合概率的界限为0和1，以及标准差的下界为0。边界向量的其余元素设置为+Inf或-Inf，表示没有限制。

lb = [0 -Inf -Inf 0 0];ub = [1 Inf Inf Inf Inf Inf];paramEsts = MLE（X，“pdf”，pdf_normmixture，“开始”开始,…“低”，磅，“上”乌兰巴托)

警告:极大似然估计不收敛。迭代超过限制。

paramEsts =1×50.3523 0.0257 3.0489 1.0546 1.0629

定制发行版的默认值是200次迭代。

statset（'mlecustom')

ans =结构体字段:显示： '断开' MaxFunEvals：400 MAXITER：200 TolBnd：1.0000e-06 TolFun：1.0000e-06 TolTypeFun：[] TolX：1.0000e-06 TolTypeX：[] GradObj： '断开' 雅可比：[] DerivStep：6.0555e-06 FunValCheck： '上' 鲁棒：[] RobustWgtFun：[] WgtFun：[]调谐：[] UseParallel：[] UseSubstreams：[]流：{} OutputFcn：[]

方法创建的选项结构重写该默认值statset函数。也增加了功能的评估限制。

选择= statset (“麦克斯特”,300,“MaxFunEvals”，600）;paramEsts = MLE（X，“pdf”，pdf_normmixture，“开始”开始,…“低”，磅，“上”乌兰巴托,“选项”选项)

paramEsts =1×50.3523 0.0257 3.0489 1.0546 1.0629

看来，收敛的最终迭代只在结果的最后几个数字中起作用。尽管如此，确保达到一致始终是一个好主意。

最后，我们可以根据原始数据的概率直方图绘制拟合密度图，直观地检验拟合。

垃圾箱= -2.5:.5:7.5;h =酒吧(垃圾箱,histc (x,垃圾箱)/(长度(x) * 5),“histc”);h。FaceColor = [。9。9。9);xgrid = linspace(1.1 *分钟(x), 1.1 *马克斯(x), 200);pdfgrid = pdf_normmixture (xgrid paramEsts (1) paramEsts (2), paramEsts (3), paramEsts (4), paramEsts (5));持有上情节(xgrid pdfgrid,“- - -”)举行从包含(“x”)ylabel (的概率密度)

使用嵌套函数:精度不相等的普通例子

在收集数据时，有时会发生这样的情况，即每次观测的精度或可靠性都不同。例如，如果几个实验者各自进行了相同数量的独立测量，但每个人只报告测量的平均值，那么每个报告数据点的可靠性将取决于原始观测数据的数量。如果原始数据不可用，对其分布的估计必须考虑这样一个事实，即可用的数据，即平均值，每一个都有不同的方差。该模型对极大似然参数估计有显式求解。但是，出于演示的目的，我们将使用MLE估计参数。

假设我们有10个数据点，每个点实际上是1到8个观察值的平均值。这些原始的观察结果是不可用的，但是我们知道每个数据点的数量。我们需要估计原始数据的均值和标准差。

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

每个数据点的方差与观测数据的个数成反比，因此我们将使用1/m来加权每个数据点的方差，使其符合最大似然拟合。

W = 1 ./米

W =1×100.1250 0.5000 1.0000 0.3333 0.1250 0.2500 0.5000 0.2000 0.5000 0.2500

在我们这里拟合的模型中，我们可以通过其PDF定义分布，但是使用日志PDF更自然一些，因为普通的PDF是这种形式的

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

和MLE无论如何都要取PDF的log来计算log-likelihood。因此，我们将创建一个直接计算日志PDF的函数。

在给定mu和的情况下，log PDF函数必须计算x中每个点的概率密度的log。它还需要考虑不同的方差权重。与前面的示例不同，这个分布函数比一行程序稍微复杂一些，并且在它自己的文件中作为单独的函数编写。因为log PDF函数需要将观察值作为额外的数据，所以最直接的方法就是使用嵌套函数来实现这种匹配。

我们为一个名为的函数创建了一个单独的文件wgtnormfit.m。此函数包含数据初始化、用于加权普通模型中的日志PDF的嵌套函数和对MLE函数实际适合模型。因为必须是正的，我们必须指定较低的参数界限。调用MLE返回单个向量中mu和的最大似然估计值。

类型wgtnormfit.m

加权正态分布的函数参数= wgtnormfit % wgtnormfit拟合示例。% MathWorks, Inc. x = [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]';函数logy = logpdf_wn(x,mu,sigma) v = sigma。^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - 5.*log(2.* .*v);结束paramEsts =大中型企业(x, logpdf, @logpdf_wn,‘开始’,(意思是(x)性病(x)],“低”,[负0]);结束

在wgtnormfit.m,我们通过MLE一个句柄嵌套函数logpdf_wn，使用“logpdf”参数。这嵌套函数是指观察计数，男，在加权日志PDF的计算。因为向量m被其父函数定义logpdf_wn可以访问它，并且不需要担心将m作为显式输入参数传入。

我们需要为参数估计提供一个粗略的初步猜测。这种情况下，非加权样本均值和标准差应该是可以的wgtnormfit.m用途。

开始= [(x),性病(x))

开始=1×21.0280 - 1.5490

为了适应模型，我们运行的拟合函数。

paramEsts = wgtnormfit

paramEsts =1×20.6244 - 2.8823

注意，亩的估计值小于三分之二样品平均的。这只是因为它应该是：估计是最可靠的数据点的影响最大，即是基于原始观测人数最多的人。在此数据集，这些点往往从加权样本均值拉下来估计。

使用参数转换:正常示例(续)

在极大似然估计中，参数的置信区间通常是使用估计器分布的大样本正态逼近来计算的。这通常是一个合理的假设，但对于小样本，有时通过转换一个或多个参数来改进正态逼近是有利的。在本例中，我们有一个location参数和一个scale参数。尺度参数经常被转换成它们的对数，我们在这里用sigma来做。首先，我们将创建一个新的log PDF函数，然后使用该参数重新计算估计值。

新的日志PDF功能创建为函数内的嵌套函数wgtnormfit2.m。作为第一嵌入，该文件包含的数据初始化，在加权正常模式日志PDF的嵌套函数，并以呼叫MLE函数实际适合模型。因为可以是任何正值，log()是无界的，我们不再需要指定下界或上界。还有，打电话给MLE在本例中，返回参数估计和置信区间。

类型wgtnormfit2.m

函数[paramEsts,paramCIs] = wgtnormfit2 % wgtnormfit2加权正态分布(log(sigma)参数化)拟合示例。% MathWorks, Inc. x = [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]';函数logy = logpdf_wn2(x,mu,logsigma) v = exp(logsigma)^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - 5.*log(2.* .*v);end [paramEsts,paramCIs] = mle(x， 'logpdf'，@logpdf_wn2， 'start'，[mean(x)，log(std(x))];结束

注意wgtnormfit2.m使用相同的起点，转换为新的参数，即，样本标准差的对数。

开始= [(x),日志(std (x)))

开始=1×21.0280 - 0.4376

[paramEsts, paramCIs] = wgtnormfit2

paramEsts =1×20.6244 1.0586

paramCIs =2×2-0.2802 0.6203 1.5290 1.4969

由于参数的使用日志（SIGMA），我们必须转变回原来的规模得到西格玛估计和置信区间。注意这两个mu和sigma估计是一样的首次适应，因为最大似然估计是不变的参数。

muHat = paramEsts (1)

muHat = 0.6244

sigmaHat = exp (paramEsts (2))

sigmaHat = 2.8823

muCI = paramCIs (: 1)

muCI =2×1-0.2802 - 1.5290

sigmaCI = EXP（paramCIs（：，2））

sigmaCI =2×11.8596 - 4.4677