拟合自定义分布：零截断的泊松例

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

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

RNG（18，'twister'）;n = 75;lambda = 1.75;X = Poissrnd（Lambda，N，1）;

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

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

ans = 65.

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

hist（x，[0：1：max（x）+1]）;

第一步是通过其概率函数（PF）来定义零截断的泊松分布。为Poisson分布的平均参数Lambda的值，我们将创建一个函数来计算x中的每个点的概率。零截断泊松的PF只是通常的泊松PF，称为它，使其总和。使用零截断，重新运行仅为1-PR {0}。为PF创建功能的最简单方法是使用匿名功能。

pf_truncpoiss = @（x，lambda）poisspdf（x，lambda）./（1-poisscdf（0，lambda））;

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

下一步是为参数lambda提供一个合理的粗略的初步猜测。这里，我们用样本均值。

start =均值（x）

开始= 2.2154

我们提供m使用数据和匿名功能使用“PDF”参数。（泊松是离散的，所以这真的是一个概率函数，而不是PDF。）因为泊松分布的平均参数必须是正的，我们还指定Lambda的下限。m返回Lambda的最大似然估计，以及可选地，接近参数的95％置信区间。

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

lambdahat = 1.8760.

lambdaci =2×11.4990 - 2.2530

请注意，参数估计小于样本均值。这就像它一样，因为最大的似然估计估计占据数据中缺失的零的估计。

我们还可以使用返回的大样本方差近似来计算Lambda的标准错误估计MLECOV.。

var = mlecov(，'pdf', pf_truncpoiss);stderr =√阿瓦尔人

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 < xTrunc的正态分布相同的分布来拟合它们，但大于xTrunc的概率为零。这样，我们就可以在考虑“缺失的尾巴”的情况下估计正态参数mu和sigma。

SICE（X，（ -  10：.5：4））;

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

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

截断点，xtrunc未被估计，因此它不是PDF函数的输入参数列表中的分发参数。xtrunc也不是数据矢量输入参数的一部分。通过匿名函数，我们可以简单地引用工作区中已定义的变量XTrunc，并且不需要担心将其传递为额外的参数。

我们还需要为参数估计提供粗略的启动猜测。在这种情况下，因为截断不是太端，所以样本均值和标准偏差可能很好。

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

开始=1×20.2735 2.2660

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

[paramests，paramcis] = mle（x，'pdf'，pdf_truncnorm，'开始'，开始，“低”，[ -  INF 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

stderr = sqrt（diag（acov））

stderr =2×10.5876 0.4143

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

一些数据集表现出双模态，甚至多模态，对这些数据拟合标准分布通常是不合适的。然而，简单的单峰分布的混合通常可以很好地模拟这样的数据。事实上，甚至可以根据特定于应用程序的知识来解释混合物中每个组件的来源。

在这个例子中，我们将把两个正常分布的混合符合到一些模拟数据。可以使用以下构造定义来描述该混合物，用于生成随机值：

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

对于此示例，我们将从学生的T发行版的混合生成数据，而不是使用与我们拟合相同的模型。这是您在Monte-Carlo仿真中可能做的事情，以测试拟合方法如何从模型的假设偏离拟合的假设。但是，在这里，我们将仅适用一个模拟数据集。

X = [trnd(20,1,50) trnd(4, 1100)+3];嘘(x, -2.25: .5:7.25);

与前面的示例一样，我们将通过创建一个计算概率密度的函数来定义模型。两个正态分量的混合PDF只是两个正态分量的PDF的加权和，由混合概率加权。这个PDF非常简单，可以使用匿名函数创建。该函数接受6个输入:一个用于评估PDF的数据向量，以及分布的5个参数。每个分量都有其均值和标准差的参数;混合概率是5。

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

我们还需要对参数进行初步猜测。一个模型的参数越多，合理的起点就越重要。对于这个例子，我们将从相同的正态数混合物(p = 0.5)开始，集中在数据的两个四分之一处，具有相同的标准差。标准偏差的起始值来自混合物的方差公式，该公式由每个成分的均值和方差表示。

pstart = .5;mustart = smartile（x，[。25 .75]）

MustArt =.1×20.5970 3.2456

sigmastart = sqrt（var（x） -  .25 * diff（mustart）。^ 2）

sigmastart = 1.2153.

start = [pstart mustart sigmastart sigmastart];

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

lb = [0 -Inf -Inf 0 0];ub = [1 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 =.结构与字段：显示：'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 ('maxiter', 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 = BAR（箱，HISTC（x，箱）/（长度（x）*。5），'histc'）;h.facecolor = [.9 .9 .9];XGrid = Linspace（1.1 * min（x），1.1 * max（x），200）;pdfgrid = pdf_normmixture（Xgrid，Paramests（1），Paramests（2），Paramests（3），Paramests（4），Paramests（5））;抓住在情节(xgrid pdfgrid,' - ')举行离开包含('X') ylabel (的概率密度）

使用嵌套函数:精度不等的正常例子

有时，当收集数据时，每次观察的精确度和可靠性都不一样。例如，如果几个实验者各自进行了一些相同数量的独立测量，但每个人只报告了测量值的平均值，那么每个报告的数据点的可靠性将取决于原始观测数据的数量。如果原始的原始数据无法获得，那么对它们分布的估计必须考虑到这样一个事实:可用的数据，平均值，每个都有不同的方差。这个模型实际上对极大似然参数估计有一个显式解。然而，为了说明的目的，我们将使用m估计参数。

假设我们有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 ./ 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的值。它还需要考虑不同的方差权重。与前面的示例不同的是，这个分布函数比一行代码稍微复杂一些，并且很清楚地作为单独的函数编写在它自己的文件中。因为log PDF函数需要将观察计数作为额外的数据，所以实现这种匹配的最直接的方法是使用嵌套函数。

我们为调用的函数创建了一个单独的文件wgtnormfit.m.。此功能包含数据初始化，在加权正常模型中的日志PDF嵌套功能，以及对此的呼叫m函数来适应模型。因为Sigma必须是正的，所以我们必须指定更低的参数界限。呼吁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，mu，sigma）v = sigma。^ 2 ./ m;logy =  - （x-mu）。^ 2 ./（2. * v） -  .5。* log（2. * pi。* v）;End Paramests = MLE（x，'logpdf'，@ logpdf_wn，'start'，[均值（x），std（x），'dight'，[ -  inf，0]）;结尾

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

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

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

开始=1×21.0280 - 1.5490

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

paramEsts = wgtnormfit

paramEsts =1×20.6244 2.8823

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

使用参数转换：正常示例（续）

在最大似然估计中，通常使用用于分布估计器的大样本正常近似来计算参数的置信区间。这通常是合理的假设，而是具有小的样本尺寸，有时是有利的，通过改变一个或多个参数来改善正常近似。在此示例中，我们有一个位置参数和比例参数。比例参数通常会转换为他们的日志，我们将在这里用西格玛进行。首先，我们将创建一个新的日志PDF函数，然后使用该参数化重新计算估计值。

新的日志PDF函数在函数中创建为嵌套功能wgtnormfit2.m。如在第一个拟合中，此文件包含数据初始化，在加权正常模型中的日志PDF的嵌套功能，以及对此的呼叫m函数来适应模型。因为可以是任何正数，log()是无界的，我们不再需要指定上下限。还有，号召m在本例中，返回参数估计值和置信区间。

类型wgtnormfit2.m

功能[paramests，paramcis] =加权正态分布的WGTNORMFIT2％WGTNORMFIT2装配示例（LOG（SIGMA）参数化）。％版权所有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_wn2（x，mu，logsigma）v = exp（logsigma）。^ 2 ./ m;logy =  - （x-mu）。^ 2 ./（2. * v） -  .5。* log（2. * pi。* v）;结束[paramests，paramcis] = mle（x，'logpdf'，@ logpdf_wn2，'start'，[均值（x），log（std（x））]）;结尾

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

start = [均值（x），log（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），因此我们必须转换回原始比例以获得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