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

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

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

rng (18,“旋风”）;n = 75;λ= 1.75;x = poissrnd(λ,n, 1);

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

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

ans = 65

这是这些模拟数据的直方图。请注意，除了没有零之外，这些数据看起来很像泊松分布。我们将用一个分布来拟合它与正整数上的泊松分布相同，但在0处没有概率。这样，我们可以估计泊松参数lambda，同时考虑“缺失零”。

嘘(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)

开始= 2.2154

我们提供大中型企业使用'pdf'参数对数据和匿名函数进行处理。(泊松是离散的，所以这实际上是一个概率函数，而不是PDF。)因为泊松分布的平均参数必须是正的，我们也指定了的下界。大中型企业返回lambda的最大似然估计，并可选地返回参数的近似95%置信区间。

[lambdaHat, lambdaCI] =大中型企业(x,“pdf”pf_truncpoiss,“开始”开始,“低”, 0)

lambdaHat = 1.8760

lambdaCI =2×11.4990 2.2530

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

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

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

stderr = 0.1923

为分布函数提供附加值:截断的正态示例

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

对于这个例子，我们模拟的数据来自截断的正态分布。首先，我们生成一些随机正态数据。

n = 75;μ= 1;σ= 3;x = normrnd(μ、σ,n, 1);

接下来，我们删除所有超出截断点xTrunc的观测值。在整个示例中，我们假设xTrunc是已知的，不需要估计。

xTrunc = 4;x = x(x < x trunc);长度(x)

ans = 64

这是这些模拟数据的直方图。我们将用一个与x < xTrunc的正态分布相同的分布来拟合它们，但大于xTrunc的概率为零。这样，我们就可以在考虑“缺失的尾巴”的情况下估计正态参数mu和sigma。

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

与前面的示例一样，我们将通过其PDF定义截断的正态分布，并创建一个函数来计算x中每个点的概率密度，给定参数mu和sigma的值。在截断点固定和已知的情况下，被截断的法线的PDF只是普通的法线PDF，被截断，然后重新规格化，以便它集成为1。重正化就是在xTrunc处计算的CDF。为了简单起见，我们假设所有的x值都小于xTrunc，而不进行检查。我们将使用一个匿名函数来定义PDF。

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

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

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

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

开始=1×20.2735 - 2.2660

我们提供大中型企业使用'pdf'参数对数据和匿名函数进行处理。因为必须是正的，我们也指定了参数的下界。大中型企业返回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（参数，x，“pdf”pdf_truncnorm)

acov =2×20.3452 0.1660 0.1660 0.1717

stderr =√诊断接头(acov))

stderr =2×10.5876 - 0.4143

拟合一个更复杂的分布:两个正态的混合

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

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

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

对于本例，我们将从混合的Student的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 = 5;muStart =分位数(x,[。25。)

muStart =1×20.5970 - 3.2456

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

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大中型企业的(x) =“pdf”pdf_normmixture,“开始”开始,．..“低”，磅，“上”乌兰巴托)

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

paramEsts =1×50.3523 0.0257 3.0489 1.0546 1.0629

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

statset (“mlecustom”）

ans =结构体字段:Display: 'off' MaxFunEvals: 400 MaxIter: 200 TolBnd: 1.0000e-06 TolFun: 1.0000e-06 TolTypeFun: [] TolX: 1.0000e-06 TolTypeX: [] GradObj: 'off' Jacobian:[]衍生步骤:6.0555e-06 FunValCheck: 'on' Robust: []

属性创建的选项结构覆盖该默认值statset函数。还增加了函数的求值限制。

选择= statset (“麦克斯特”, 300,“MaxFunEvals”, 600);paramEsts大中型企业的(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 (的概率密度）

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

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

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

x=[0.25-1.241.381.39-1.432.793.520.921.441.26]；m=[8213824]；

每个数据点的方差与进入该数据点的观测数成反比，因此我们将使用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)，

和大中型企业需要对PDF取对数，来计算对数可能性。因此，我们将创建一个直接计算日志PDF的函数。

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

我们已经为一个名为wgtnormfit.m．该函数包含数据初始化、加权正常模型中用于日志PDF的嵌套函数和对大中型企业函数来适应模型。因为必须是正的，我们必须指定参数的下界。调用大中型企业返回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,我们通过大中型企业嵌套函数的句柄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的嵌套函数以及对大中型企业函数来适应模型。因为可以是任何正数，log()是无界的，我们不再需要指定上下限。还有，号召大中型企业在本例中，返回参数估计值和置信区间。

类型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)，我们必须转换回原始尺度，以得到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