拟合零截断泊松分布

计数数据通常使用泊松分布建模，您可以使用poissfit或fitdist函数拟合泊松分布。然而，在某些情况下，数据中没有记录为零的计数，因此由于缺少零，拟合泊松分布并不简单。在这种情况下，通过使用泊松分布拟合零截断数据大中型企业函数和自定义分布函数。

首先，生成一些随机泊松数据。

rng (18,“旋风”）%的再现性λ= 1.75;n = 75;x1 = poissrnd(λ,n, 1);

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

X1 = X1 (X1 > 0);

检查样品的数量x1后截断。

长度(x1)

ans = 65

绘制模拟数据的直方图。

直方图(x1, 0:1:马克斯(x1) + 1)

数据看起来像泊松分布，只是它不包含0。可以使用与正整数上的泊松分布相同的自定义分布，但在0处没有概率。通过使用自定义分布，可以估计泊松参数λ同时解释缺失的零。

您需要通过其概率质量函数(pmf)来定义零截断泊松分布。创建一个匿名函数来计算每个点的概率x1，给定泊松分布均值参数的值λ．零截断泊松分布的pmf是泊松pmf归一化，使其和为1。截断为零时，归一化是1-Probability (x1 < 0)．

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

为简单起见，假设所有的x1给这个函数的值都是正整数，没有校验。对于错误检查或需要多行代码的更复杂的发行版，必须在单独的文件中定义函数。

为参数找到一个合理的粗略的第一猜测λ．在这种情况下，使用样本均值。

开始=意味着(x1)

开始= 2.2154

提供大中型企业用数据，自定义PMF函数，初始参数值，参数下界。因为泊松分布的均值参数必须为正，所以还需要指定的下界λ．的大中型企业函数返回的最大似然估计λ，以及可选的参数的近似95%置信区间。

[lambdaHat, lambdaCI] =大中型企业(x1,“pdf”pf_truncpoiss,“开始”开始,.．.下界的, 0)

lambdaHat = 1.8760

lambdaCI =2×11.4990 - 2.2530

参数估计值小于样本均值。最大似然估计解释了数据中不存在的零。

方法也可以指定截断边界TruncationBounds名称-值参数。

[lambdaHat2, lambdaCI2] =大中型企业(x1,“分布”，“泊松”，.．.“TruncationBounds”[0正])

lambdaHat2 = 1.8760

lambdaCI2 =2×11.4990 - 2.2530

你也可以计算一个标准误差估计λ通过使用返回的大样本方差近似mlecov．

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

stderr = 0.1923

为了直观地检查拟合，将拟合的pmf与原始数据的归一化直方图进行绘制

直方图(x1,“归一化”，“pdf”) xgrid = min(x1):max(x1);pmfgrid = pf_truncpoiss (xgrid lambdaHat);持有在情节(xgrid pmfgrid,“- - -”)包含(x1的) ylabel (“概率”)传说(样本数据的，的安装及，“位置”，“最佳”)举行从

拟合截断上正态分布

连续数据有时会被截断。例如，由于数据收集的限制，大于某个固定值的观测值可能不会被记录。

在本例中，模拟截断正态分布的数据。首先，生成一些随机的正态数据。

n = 500;μ= 1;σ= 3;rng (“默认”）%的再现性x2 = normrnd(μ、σ,n, 1);

接下来，删除任何超出截断点的观察值xTrunc．假设xTrunc是一个您不需要估计的已知值。

xTrunc = 4;x2 = x2(x2 < xTrunc);

检查样品的数量x2后截断。

长度(x2)

ans = 430

创建模拟数据的直方图。

直方图(x2)

用与正态分布相同的自定义分布拟合模拟数据x2 < xTrunc，但上面的概率为0xTrunc．通过使用自定义分布，可以估计正态参数μ而且σ同时解释失踪的尾巴。

定义截断正态分布的pdf。在给定的参数值下，创建一个匿名函数来计算x中每个点的概率密度值μ而且σ．在截断点固定且已知的情况下，截断正态分布的pdf就是截断后的pdf，然后归一化，使其积分为1。归一化是cdf在xTrunc．为了简单起见，假设所有x2值小于xTrunc,没有检查。

pdf_truncnorm = @ (x2,μ、σ).．.normpdf (x2,μ、σ)。/ normcdf (xTrunc、μ、σ);

因为你不需要估计截断点xTrunc，它不包含在自定义PDF函数的输入分布参数中。xTrunc也不是数据向量输入参数的一部分。匿名函数可以访问工作空间中的变量，因此不必传递xTrunc作为附加参数的匿名函数。

为参数估计提供一个粗略的初始猜测。在这种情况下，因为截断不是极值的，所以使用样本均值和标准差。

开始= [(x2),性病(x2)]

开始=1×20.1585 - 2.4125

提供大中型企业使用数据、自定义PDF函数、初始参数值和参数的下界。因为σ必须为正，还需要指定较低的参数边界。大中型企业的最大似然估计μ而且σ作为一个单独的向量，以及两个参数的近似95%置信区间的矩阵。

[paramEsts, paramCIs] =大中型企业(x2,“pdf”pdf_truncnorm,“开始”开始,.．.下界的(从0))

paramEsts =1×21.1298 - 3.0884

paramCIs =2×20.5713 2.7160 1.6882 3.4607

的估计μ而且σ大于样本均值和标准差。模型拟合解释了分布中缺失的上尾。

方法也可以指定截断边界TruncationBounds名称-值参数。

[paramEsts2, paramCIs2] =大中型企业(x2,“分布”，“正常”，.．.“TruncationBounds”(负无穷xTrunc))

paramEsts2 =1×21.1297 - 3.0884

paramCIs2 =2×20.5713 2.7160 1.6882 3.4607

你可以为参数估计计算一个近似的协方差矩阵mlecov．这种近似通常适用于大样本，可以通过对角线元素的平方根近似标准误差。

acov = mlecov (paramEsts x2,“pdf”pdf_truncnorm)

acov =2×20.0812 0.0402 0.0402 0.0361

stderr =√诊断接头(acov))

stderr =2×10.2849 - 0.1900

为了直观地检查适合度，可以根据原始数据的规范化直方图绘制适合度的pdf。

直方图(x2,“归一化”，“pdf”) xgrid = linspace(min(x2)，max(x2));pdfgrid = pdf_truncnorm (xgrid paramEsts (1) paramEsts (2));持有在情节(xgrid pdfgrid,“- - -”)包含(“x2”) ylabel (的概率密度)传说(样本数据的，“安装pdf”，“位置”，“最佳”)举行从

拟合两个正态分布的混合

有些数据集表现为双模态，甚至多模态，因此为这些数据拟合标准分布通常是不合适的。然而，简单的单峰分布的混合通常可以很好地模拟这类数据。

在这种情况下，拟合两个正态分布的混合模拟数据。考虑具有以下建设性定义的模拟数据:

从Student的混合数据中生成一个数据集t分布而不是使用你拟合的相同模型。通过使用不同的分布(类似于蒙特卡洛模拟中使用的技术)，您可以测试拟合方法对偏离拟合模型的假设的鲁棒性。

rng (10)%的再现性X3 = [trnd(20,1,50) trnd(4,1,100)+3];直方图(x3)

通过创建一个计算概率密度的匿名函数来定义拟合模型。两个正态分布的混合物的pdf是两个正态组分的pdf的加权和，由混合概率加权。匿名函数接受6个输入:一个用于评估pdf的数据向量和5个分布参数。每个分量都有其均值和标准差的参数。

pdf_normmixture = @ (x3, p, mu1、mu2 sigma1, sigma2).．.p * normpdf (x3, mu1 sigma1) + (1 - p) * normpdf (x3, mu2 sigma2);

您还需要对参数进行初步猜测。随着模型参数数量的增加，定义起点变得更加重要。先从等量的混合物开始(p= 0.5)的正态分布，以数据的两个四分之一为中心，标准差相等。标准差的起始值来自于混合物的方差公式，由每个成分的均值和方差表示。

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

muStart =1×20.3351 - 3.3046

sigmasterart =√(var(x3) - .25*diff(muStart).^2)

sigmaStart = 1.1602

start = [pStart muStart sigmasterart];

为混合概率指定0和1的界限，为标准差指定0的下界。将边界向量的剩余元素设置为+正或负，表示没有限制。

lb = [0 -Inf -Inf 0 0];ub = [1 Inf Inf Inf];paramEsts =大中型企业(x3,“pdf”pdf_normmixture,“开始”开始,.．.下界的磅,“UpperBound”乌兰巴托)

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

paramEsts =1×50.3273 -0.2263 2.9914 0.9067 1.2059

警告消息表示函数不与默认迭代设置收敛。显示默认选项。

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:[]衍生step: 6.0555e-06 FunValCheck: 'on'健壮:[]RobustWgtFun: [] WgtFun: [] Tune: [] UseParallel: [] UseSubstreams: [] Streams: {} OutputFcn: []

自定义发行版的默认最大迭代次数是200。方法创建的选项结构，覆盖默认值以增加迭代次数statset函数。另外，增加函数的最大求值。

选择= statset (“麦克斯特”, 300,“MaxFunEvals”, 600);paramEsts =大中型企业(x3,“pdf”pdf_normmixture,“开始”开始,.．.下界的磅,“UpperBound”乌兰巴托,“选项”选项)

paramEsts =1×50.3273 -0.2263 2.9914 0.9067 1.2059

收敛的最后迭代只在结果的最后几位数有意义。然而，最佳实践是始终确保达到收敛。

为了直观地检查拟合性，可以根据原始数据的概率直方图绘制拟合密度。

直方图(x3,“归一化”，“pdf”)举行在xgrid = linspace(1.1 *分钟(x3), 1.1 *马克斯(x3), 200);pdfgrid = pdf_normmixture (xgrid,.．.paramEsts paramEsts (1) (2), paramEsts (3), paramEsts (4), paramEsts (5));情节(xgrid pdfgrid,“- - -”)举行从包含(“x3”) ylabel (的概率密度)传说(样本数据的，“安装pdf”，“位置”，“最佳”）

或者，对于正态分布的混合，您可以使用fitgmdist函数。由于初始估计和迭代算法的设置不同，估计可能不同。

Mdl = fitgmdist (x3 ', 2)

Mdl =一维2组分的高斯混合分布。组分1:混合比例:0.329180平均值:-0.2200组分2:混合比例:0.670820平均值:2.9975

Mdl。σ

Ans = Ans (:，:，1) = 0.8274 Ans (:，:，2) = 1.4437

加权正态分布拟合精度不等的数据

假设您有10个数据点，其中每个点实际上是1到8个观测值的平均值。没有原始观测数据，但每个数据点的观测数是已知的。每个点的精度取决于相应的观测次数。你需要估计原始数据的平均值和标准差。

X4 = [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. /米;

在这个模型中，您可以通过它的pdf定义分布。但是，使用pdf的对数更合适，因为正常的pdf具有这种形式

C .* exp(-0.5 .* z.²)，

而且大中型企业获取PDF的对数，以计算对数可能性。因此，相反，创建一个函数直接计算pdf的对数。

pdf函数的对数必须计算每个点的概率密度的对数x，给定正态分布参数μ而且σ．它还需要考虑到不同的方差权重。定义一个名为helper_logpdf_wn1在一个单独的文件中helper_logpdf_wn1.m．

函数呆呆的= helper_logpdf_wn1 (x, m,μ、σ)权重正态分布的pdf的对数此函数只支持示例Fit自定义分布金宝app% (customdist1demo.mlx)，并可能在未来的版本中更改。v =σ。^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - .5.*log(2.*pi.*v);结束

为参数估计提供一个粗略的第一猜测。在这种情况下，使用未加权样本均值和标准差。

开始= [(x4),性病(x4)]

开始=1×21.0280 - 1.5490

因为σ必须为正，您需要指定较低的参数边界。

[paramEsts1, paramCIs1] =大中型企业(x4,“logpdf”，.．.@ (x,μ、σ)helper_logpdf_wn1 (x, m,μ、σ),.．.“开始”开始,下界的(从0))

paramEsts1 =1×20.6244 - 2.8823

paramCIs1 =2×2-0.2802 1.6191 1.5290 4.1456

的估计μ小于样本均值估计的三分之二。这一估计受到最可靠的数据点的影响，即根据最多原始观测所得的数据点。在这个数据集中，这些点倾向于将估计值拉低至未加权样本均值。

利用参数变换拟合正态分布

的大中型企业函数计算参数的置信区间，如果没有精确的方法，则使用大样本正态逼近估计量的分布。对于小样本，可以通过转换一个或多个参数来改进正态近似。在这种情况下，将一个正态分布的尺度参数转换为它的对数。

首先，定义一个新的日志pdf函数helper_logpdf_wn2它使用转换后的参数σ．

函数呆呆的= helper_logpdf_wn2 (x, m,μ,logsigma)的权重正态分布的pdf的对数%日志(σ)参数化此函数只支持示例Fit自定义分布金宝app% (customdist1demo.mlx)，并可能在未来的版本中更改。v = exp (logsigma)。^ 2。/ m;呆呆的= - (xμ)。^2 ./ (2.*v) - .5.*log(2.*pi.*v);结束

使用转换为新参数化的相同起点σ，即样本标准差的对数。

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

开始=1×21.0280 - 0.4376

因为σ可以是任何正值，日志(σ)是无界的，您不需要指定下界或上界。

[paramEsts2, paramCIs2] =大中型企业(x4,“logpdf”，.．.@ (x,μ、σ)helper_logpdf_wn2 (x, m,μ、σ),.．.“开始”,开始)

paramEsts2 =1×20.6244 - 1.0586

paramCIs2 =2×2-0.2802 0.6203 1.5290 1.4969

因为参数化使用日志(σ)的估计值和置信区间σ．

sigmaHat = exp (paramEsts2 (2))

sigmaHat = 2.8823

sigmaCI = exp (paramCIs2 (:, 2))

sigmaCI =2×11.8596 - 4.4677

对两者的估计μ而且σ与第一次拟合相同，因为最大似然估计对参数化是不变的。的置信区间σ稍微不同于paramCIs1 (: 2)．