拟合自定义分布：一个零截断泊松分布示例

计数数据通常使用泊松分布建模，您可以使用统计信息和计算机学习工具箱功能Poissfit.以拟合泊松模型。然而，在某些情况下，零的计数不会记录在数据中，因此拟合泊松分布不是直接的，因为这些“缺失的零”。此示例将演示如何使用函数将泊松分布拟合到零截断数据m.

对于此示例，我们将使用来自零截断的泊松分布的模拟数据。首先，我们生成一些随机泊松数据。

rng（18，'twister'）;n = 75;拉姆达= 1.75;X = poissrnd（拉姆达，N，1）;

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

x=x（x>0）；长度（x）

ans=65

这是这些模拟数据的直方图。请注意，数据看起来非常像泊松分布，不同之处在于没有零。我们将符合与正整数上泊松的分布相同，但这并无零概率。通过这种方式，我们可以估计泊松参数lambda，同时考虑“缺少零”。

HIST（X，[0：1：最大（X）1]）;

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

pf_truncpoiss = @（X，拉姆达）poisspdf（X，拉姆达）./（1- poisscdf（0，拉姆达））;

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

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

开始=平均值（x）

start = 2.2154.

我们提供m使用数据和匿名功能使用“PDF”参数。（泊松是离散的，所以这是一个真正的概率函数，而不是一个PDF）。因为泊松分布的均值参数必须是积极的，我们也可以指定一个下界拉姆达。m返回lambda的最大似然估计值，以及参数的近似95%置信区间（可选）。

[lambdahat，lambdaci] = mle（x，“pdf”，pf_truncpoiss，'开始'，开始，'降低'，0）

lambdaHat = 1.8760

lambdaCI =2×11.4990 2.2530

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

我们还可以使用MLECOV..

avar = mlecov（lambdahat，x，“pdf”，pf_truncpois）；标准偏差=sqrt（avar）

stderr = 0.1923

向分发功能提供额外的值：截断的正常示例

有时也会发生连续数据被截断的情况。例如，由于数据收集方式的限制，可能无法记录大于某个固定值的观测值。此示例将演示如何使用函数将正态分布拟合到截断数据m.

对于此示例，我们从截短的正态分布模拟数据。首先，我们生成一些随机的正常数据。

n = 75;mu = 1;西格玛= 3;X = NORMRND（MU，SIGMA，N，1）;

接下来，我们删除掉落到截断点XTrunc的任何观察结果。在此，我们假设Xtrunc是已知的，不需要估计Xtrunc。

xtrunc = 4;x = x（x 
              

ans=64

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

HIST（X，（ -  10：0.5：4））;

与前面的示例一样，我们将通过其PDF定义截断的正态分布，并创建一个函数来计算X中每个点的概率密度，给定参数mu和sigma的值。使用截断点固定和已知，截断正常的PDF仅是通常的正常PDF，截断，然后重新成型，使其集成到一个。重整化只是在XTrunc评估的CDF。为简单起见，我们假设所有x值都小于xtrunc而不检查。我们将使用匿名函数来定义PDF。

pdf_truncnnorm=@（x，mu，sigma）normpdf（x，mu，sigma）。/normcdf（xTrunc，mu，sigma）；

该截断点，xTrunc，没有被估计，所以它是不是在PDF函数的输入参数列表中的分布参数之一。xTrunc也不是数据矢量输入参数的一部分。随着一个匿名函数，我们可以简单地指的是已在工作区中定义的变量xTrunc，也没有必要担心它传递中作为一个额外的参数。

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

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

开始=1×20.2735 2.2660

我们提供m使用数据和匿名功能使用“PDF”参数。因为Sigma必须是正面的，所以我们还指定了更低的参数界限。m作为单个向量返回mu和sigma的最大似然估计值，以及两个参数的近似95%置信区间矩阵。

[参数，参数]=mle（x，“pdf”，pdf_truncnorm，'开始'，开始，'降低'，[-Inf 0]）

paramests =1×20.9911 2.7800

帕拉米西斯=2×2-0.1605 1.9680 2.1427 3.5921

请注意，MU和SIGMA的估计比样本均值和标准偏差大得多。这是因为模型拟合已经占“缺失”的分布上尾。

参数估计值使用，我们可以计算的近似协方差矩阵MLECOV..近似通常较大样品中合理，估计的标准误差可以通过平方根对角元素来近似。

acov=mlecov（参数，x，“pdf”，pdf_truncnorm）

acov=2×20.3452 0.1660 0.1660 0.1717

stderr = sqrt（diag（acov））

标准错误=2×10.5876 0.4143

拟合更复杂的分布：两种法线的混合物

一些数据集表明双极性，甚至多模，并拟合到这些数据的标准分布通常是不合适的。然而，简单的单峰分布的混合物通常可以非常良好地模拟这些数据。实际上，甚至可以基于特定于应用的知识来对混合物中每个组分的来源进行解释。

在本例中，我们将对一些模拟数据拟合两个正态分布的混合。可使用以下生成随机值的构造性定义来描述该混合物：

首先，翻转偏置硬币。如果它落在头部，则从具有平均MU_1和标准偏差Sigma_1的正常分布随机挑选值。如果硬币降落尾，则从平均MU_2和标准偏差Sigma_2随机挑选值。

在这个例子中，我们将生成从学生的t分布的混合数据，而不是使用相同的模型，因为我们装修。这是那种你可能会在蒙特卡罗模拟做测试拟合方法如何健壮从模型拟合被的假设出发的事情。然而在这里，我们将只适合一个模拟数据集。

x=[trnd（20,1,50）trnd（4,1100）+3]；历史（x，-2.25:.5:7.25）；

与前面的示例一样，我们将通过创建一个计算概率密度的函数来定义模型。两个法线的混合的PDF只是两个法线分量的PDF的加权和，由混合概率加权。此PDF非常简单，可以使用匿名函数创建。该函数接受六个输入：用于计算PDF的数据向量和分布的五个参数。每个成分都有其平均值和标准偏差的参数；混合概率总计为五。

pdf_normmixture = @（x，p，mu1，mu2，sigma1，sigma2）...p * normpdf（x，mu1，sigma1）+（1-p）* normpdf（x，mu2，sigma2）;

我们还需要对参数进行初步猜测。模型的参数越多，合理的起点就越重要。对于本例，我们将从法线的相等混合（p=0.5）开始，以数据的两个四分位数为中心，具有相等的标准偏差。标准偏差的起始值来自混合物方差公式，即各成分的平均值和方差。

pStart=0.5；muStart=分位数（x，[25.75]）

穆斯塔特=1×20.5970 3.2456

sigmaStart = SQRT（VAR（X） -  0.25 * DIFF（muStart）^ 2）。

sigmaStart = 1.2153

start = [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 =.结构与字段：显示：'OFF'MAXFUNevals：400 MAXITER：200 TOLBND：1.0000E-06 TOLFUN：1.0000E-06 TOLTYPEFUN：[] TOLX：1.0000E-06 TOLTYPEX：[] GRAVOBJ：'OFF'jacobian：[]德语：6.0555e-06 funvalcheck：'on'robust：[] robustwgtfun：[] wgtfun：[] tune：[]使用adplelial：[]使用subsubstreams：[] streams：{} outputfcn：[]

使用使用创建的选项结构覆盖该默认值statset功能。还增加了功能评价限制。

选项= statset（'MAXITER'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=钢筋（料仓，历史（x，料仓）/（长度（x）*.5），'histc'）;h.FaceColor = [0.9 0.9 0.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（'概率密度')

使用嵌套功能：具有不等精度的正常示例

当收集数据时，它有时会发生每个观察，以不同的精度或可靠性。例如，如果几个实验者每个实验者都制作相同数量的许多独立测量，但每个报告只有其测量的平均值，每个报告的数据点的可靠性都将取决于进入它的原始观测的数量。如果原始原始数据不可用，则其分发的估计必须考虑到可用的数据，平均值，每个都具有不同的差异。该模型实际上具有用于最大似然参数估计的显式解决方案。但是，出于插图的目的，我们将使用m来估计参数。

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

x=[0.25-1.241.381.39-1.432.793.520.921.441.26]；m=[82 13 84 2 5 2 4]；

每个数据点的方差与进入它的观察数成反比，因此我们将使用1 / m重量最大似然符合的每个数据点的方差。

w = 1 ./ m

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），

和m不管怎样，都必须获取PDF的日志，以计算日志可能性。因此，我们将创建一个直接计算日志PDF的函数。

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

我们为一个名为wgtnormfit.m. 此函数包含数据初始化、加权正态模型中日志PDF的嵌套函数以及对m功能实际拟合模型。由于西格玛必须是积极的，我们必须指定较低的参数界限。呼吁m返回单个载体mu和sigma的最大似然估计。

类型wgtnormfit.m

功能paramEsts = wgtnormfit％WGTNORMFIT件例如用于加权正态分布。％版权所有1984-2012 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]';功能迟缓= logpdf_wn（X，μ，西格玛）V =Σ-^ 2 ./米。logy =  - （x-mu）。^ 2 ./（2. * v） -  .5。* log（2. * pi。* v）;结束paramEsts = MLE（X， 'logpdf'，@ logpdf_wn， '启动'，[平均（x）中，STD（x）]的， '下'，[ - 天道酬勤，0]）;结尾

在wgtnormfit.m，我们通过m嵌套功能的句柄logpdf_wn，使用“logpdf”参数。嵌套功能是指在计算加权日志PDF的计算中的观察计数m。因为矢量m在其父函数中定义，logpdf_wn可以访问它，并且不必担心将m作为显式输入参数传入。

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

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

开始=1×21.0280 1.5490

为了适应模型，我们运行拟合功能。

paramests = wgtnormfit.

paramests =1×20.6244 2.8823

请注意，MU的估计率小于样本意味着的三分之二。这就像它应该是：估计受到最可靠的数据点，即基于最大原始观察数量的数据点受到影响。在该数据集中，这些点倾向于将估计从未加权的样本平均拉下来。

使用参数变换：正常实施例（续）

在最大似然估计，对于参数置信区间使用大样品法线近似为估计量的分布通常计算。这通常是一个合理的假设，但用小样本大小，有时有利的是通过将一个或多个参数，以改善该正态近似。在这个例子中，我们有一个位置参数和尺度参数。缩放参数常常转化为自己的日志中，我们将做到这一点这里西格玛。首先，我们将创建一个新的日志PDF功能，然后使用该参数重新计算的估计。

新的日志PDF函数在函数中创建为嵌套功能wgtnormfit2.m.如在第一个拟合中，此文件包含数据初始化，在加权正常模型中的日志PDF的嵌套功能，以及对此的呼叫m功能实际拟合模型。因为Sigma可以是任何正值，所以log（sigma）是无限的，所以我们不再需要指定较低或上限。此外，呼叫m在这种情况下，返回参数估计和置信区间。

类型wgtnormfit2.m

函数[参数，参数]=wgtnormfit2%wgtnormfit2加权正态分布拟合示例（对数（西格玛）参数化）。%版权所有1984-2012 The MathWorks，Inc.x=[0.25-1.24 1.38 1.39-1.43 2.79 3.52 0.92 1.44 1.26]；m=[82 13 84 2 5 2 4]；函数logy=logpdf_wn2（x，mu，logsigma）v=exp（logsigma）。^2./m；逻辑=-（x-mu）。^2./（2.*v）-.5.*log（2.*pi.*v）；结束[参数，参数]=mle（x，'logpdf'，@logpdf_wn2，'start'，[均值（x），对数（标准（x））]；结束

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

开始=[平均值（x），对数（标准值（x））]

开始=1×21.0280 0.4376

[paramEsts，paramCIs] = wgtnormfit2

paramests =1×20.6244 1.0586

帕拉米西斯=2×2-0.2802 0.6203 1.5290 1.4969

由于参数化使用日志（Sigma），因此我们必须转换回原始比例以获得Sigma的估计和置信区间。请注意，MU和SIGMA的估计与第一次适合相同，因为最大似然估计是参数化的不变性。

Muhat = Paramests（1）

Muhat = 0.6244.

sigmaHat = EXP（paramEsts（2））

sigmaHat=2.8823

muCI=paramCIs（：，1）

粘液=2×1-0.2802 1.5290

sigmaci = exp（paramcis（:,2））

sigmaci =2×11.8596 4.4677