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

计数数据使用的是泊松分布常为蓝本，并且可以使用统计和机器学习工具箱功能Poissfit.来拟合泊松模型。然而，在某些情况下，数据中不会记录为零的计数，因此，由于这些“缺失的零”，拟合泊松分布并不简单。这个例子将展示如何使用函数将泊松分布拟合到零截尾数据MLE．

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

rng (18,'twister'）;n = 75;λ= 1.75;x = poissrnd(λ,n, 1);

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

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

ans = 65

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

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

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

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

为简单起见，我们假设给出的所有x值都是正整数，没有检查。错误检查或更复杂的分发，可能需要多行代码，表明该函数应在单独的文件中定义。

下一步是为参数lambda提供合理的粗略猜测。在这种情况下，我们只需使用样本意味着。

开始=意味着(x)

start = 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．

avar = mlecov（lambdahat，x，“pdf”，pf_truncpoiss）;stderr = sqrt（avar）

stderr = 0.1923

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

它有时也会发生连续数据被截断。例如，由于收集数据方式的限制，可能无法记录大于某些固定值的观察。此示例将展示如何使用该函数拟合正常分布以截断的数据MLE．

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

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

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

xTrunc = 4;X = X（X 
              

ans = 64

下面是这些模拟数据的直方图。我们将使用与X

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

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

pdf_trungnorm = @（x，mu，sigma）normpdf（x，mu，sigma）./ normcdf（xtrunc，mu，sigma）;

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

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

start = [均值（x），std（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随机挑选值。

对于本例，我们将从混合的Student的t分布生成数据，而不是使用我们拟合的相同模型。这是你在蒙特卡罗模拟中可能会做的事情，以测试拟合方法对偏离拟合模型假设的稳健程度。然而，在这里，我们只适合一个模拟数据集。

X = [TRND（20,1,50）TRND（4,1,100）+3];stay（x，-2.25：.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）开始，以数据的两个四分位数为中心，具有相同的标准偏差。标准偏差的起始值来自于混合物方差的公式，就每个组分的平均值和方差而异。

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

muStart =1×20.5970 - 3.2456

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

sigmaStart = 1.2153

开始= [PSTART muStart sigmaStart sigmaStart];

最后，我们需要指定零的范围，一个用于混合概率的范围，以及标准偏差的零边界。界限向量的剩余元素设置为+ INF或-INF，以表示没有限制。

lb = [0 -inf -inf 0 0];UB = [1个IM INF INF];paramEsts = MLE（X，“pdf”，pdf_normmixture，“开始”开始,．..'降低'，磅，'上'，UB）

警告：最大可能性估计没有收敛。超出迭代限额。

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，“开始”开始,．..'降低'，磅，'上'，UB，'选项'选项)

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，“- - -”） 抓住从Xlabel（“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];m = [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取对数，来计算对数可能性。因此，我们将创建一个直接计算日志PDF的函数。

Log PDF功能必须计算X中每个点的概率密度的日志，给定mu和sigma的值。它还需要考虑不同的方差权重。与前面的示例不同，此分发功能比单套接在一起更复杂，并且最清晰地写在自己文件中的单独函数。由于日志PDF功能需要观察计数作为额外数据，所以实现这一合身的最直接方式是使用嵌套功能。

我们已经为一个名为wgtnormfit.m．该函数包含数据初始化、加权正常模型中用于日志PDF的嵌套函数和对MLE函数来适应模型。因为必须是正的，我们必须指定参数的下界。调用MLE返回mu和sigma在单个向量中的最大似然估计。

类型wgtnormfit.m

加权正态分布的拟合示例。% Copyright 1984-2012 The MathWorks, Inc. x = [0.25 -1.24 1.38 1.39 2.79 3.52 0.92 1.44 1.26]';M = [8 2 1 3 8 4 2 5 2 4]';函数= logpdf_wn(x,mu,sigma) v = sigma。^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - 5.*log(2.*pi.*v);结束paramEsts =大中型企业(x, logpdf, @logpdf_wn,‘开始’,(意思是(x)性病(x)],“低”,[负0]);结束

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

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

start = [均值（x），std（x）]

开始=1×21.0280 1.5490

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

paramests = wgtnormfit.

paramests =1×20.6244 - 2.8823

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

使用参数转换:普通示例(续)

在最大似然估计中，参数的置信区间通常是用估计量分布的大样本正态近似来计算的。这通常是一个合理的假设，但在小样本容量的情况下，有时通过转换一个或多个参数来改进正态逼近是有利的。在这个例子中，我们有一个位置参数和一个比例参数。尺度参数通常转换成对数，我们这里用。首先，我们将创建一个新的日志PDF函数，然后使用该参数重新计算估计。

新的日志PDF功能创建为函数内的嵌套函数WGTNORMFIT2.M.．作为第一嵌入，该文件包含的数据初始化，在加权正常模式日志PDF的嵌套函数，并以呼叫MLE函数来适应模型。因为Sigma可以是任何正值，所以log（sigma）是无限的，所以我们不再需要指定较低或上限。此外，呼叫MLE在这种情况下，返回参数估计和置信区间。

类型WGTNORMFIT2.M.

加权正态分布(log(sigma)参数化)的拟合示例。% Copyright 1984-2012 The MathWorks, Inc. x = [0.25 -1.24 1.38 1.39 2.79 3.52 0.92 1.44 1.26]';M = [8 2 1 3 8 4 2 5 2 4]';函数logpdf_wn2(x,mu,logsigma) v = exp(logsigma)。^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - 5.*log(2.*pi.*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