最大似然估计gamma分布参数(附MATLAB完整代码)

如何使用最大似然估计法（MLE）来估计伽马分布的参数:

伽马分布是一个连续型概率分布，其概率密度函数为：

f(x|α,β) = (β^α / Γ(α)) x^(α-1) e^(-βx) ，其中, Γ(α) 是伽马函数，x > 0, α, β > 0.

假设我们有一个由n个独立同分布的随机样本X1, X2, …, Xn构成的样本，它们都服从参数为α和β的伽马分布。那么，这n个样本的联合概率密度函数将是每个样本概率密度函数的乘积，即似然函数L(α,β)为：

L(α,β) = Π (β^α / Γ(α)) x^(α-1) e^(-βx)

为了方便求解，我们通常使用对数似然函数，即：

lnL(α,β) = Σ [αlnβ + (α - 1)lnx - βx - lnΓ(α)]

我们的目标是找到α和β的值，使得这个函数取最大值。这就需要我们对α和β分别求偏导，并令其等于0，得到一组联立方程，解这组方程我们可以得到α和β的估计值。

然而，这个问题的解析解是无法直接求得的，需要依赖于数值算法（例如：牛顿迭代法）来进行求解。在实际的应用中，我们通常会使用一些MATLAB或Python来完成这个任务。

MATLAB代码如下:

clc;close all;clear all;warning off;%清除变量
format long g;

adata= gamrnd(5,2,1000,1)

% [phat,pci] = gamfit(adata,0.95);
% phat
x0=[1,1];
y=adata;
[parmhat1,ML]=fminsearch(@(x) mygammapdffun(y,x),x0);
disp('最大似然估计得到的gamma分布参数');
parmhat1

%% 画概率图
[counts,centers]=hist(y,100);
figure;
bar(centers,counts/sum(counts)); %画出概率密度分布图
xlabel('额度','fontname','宋体');
ylabel('概率','fontname','宋体');
title('分布图','fontname','宋体');

a1=gamcdf(sort(y),parmhat1(1),parmhat1(2));

[counts,centers]=hist(y,100);
a0=counts/sum(counts);
g=cumsum(a0);

figure;
plot(centers,g,'b*'); %画出概率密度分布图
hold on;
plot(sort(y),a1,'r');
legend({'样本','拟合'},'fontname','宋体');
xlabel('数值','fontname','宋体');
ylabel('累积概率','fontname','宋体');
title('gamma分布累积概率','fontname','宋体');

function y=mygammapdffun(data,x)
% gamma分布
a=x(1);% 分布的参数1
b=x(2);% 分布的参数2
p=1/(b^a*gamma(a)).*data.^(a-1).*exp(-data./b);% gamma分布的概率密度函数
y=-sum(log(p));% 似然函数值最大化
 

数据如下:

adata =

          3.18488078866483
          5.03069162867442
          7.51749037559498
          10.2835778803303
           21.603895294973
          14.1441113119225
          14.7865293223264
          8.14067810458262
          5.06282895822804
          7.65819802379196
          13.5270169969674
          17.0370414842674
          12.3091911437045
          11.1321241815143
          9.59563779171205
          7.05518598131705
           7.3971689150796
          6.08394395856918
          6.81302478109756
          14.3572119616325
          17.8091493287761
          20.1847205029288
          2.75323227564708
          2.64232307308316
          22.8297342219077
          11.5996016372587
          7.80526023064686
          7.80557373707557
          5.65136312002752
          8.13348512874601
          9.72119202730636
          3.89003555558571
          6.27485075615423
          9.55033639102666
          14.2806247237316
          4.60212965022015
          7.90983340245262
          6.50889311430597
          10.4755128619394
          6.57234682832661
          8.59327563483049
          4.33742254479794
          6.39067595216208
          14.2185641379183
           6.5362012959033
          2.59374764667069
          3.76193744697243
          18.2185228490581
          7.51517323567695
          13.3130107547276
          24.8864878455546
          5.14189489625271
           5.6518406684059
          13.8269601152627
          9.63129826361994
          4.65958302250466
          3.80216423522717
          5.97345782586643
          11.5473765605036
          4.75931300049055
           11.364680079473
          6.93790347627799
          7.58815633610866
          11.7172324209223
          6.91995739570363
          8.59110905917763
          7.71543620382142
           6.5021785418607
          19.5628272958143
          12.0078420961938
          6.97438879245941
           15.123718686986
          5.18891661218098
          17.3260811939267
          16.1996830361599
          11.5501158462896
          3.34981290836644
          9.50304758087526
          7.15333184666909
          11.7674131865076
          12.1311270077904
          8.34965382173295
          6.16462118334621
          5.89701477494658
          6.82509503581429
          3.64097940761756
          11.8950632662039
          10.3211403069445
          13.4413727934656
          17.2964717011917
          11.8610426129828
          5.91299332526704
          4.57398306654572
          4.77286770722602
          9.63135879012397
          8.39949078456355
          12.3286922969972
          10.4168057540566
          12.2717592257667
          8.19011381460088

程序结果如下:

最大似然估计得到的gamma分布参数

parmhat1 =

          4.24093001674139          2.24584755857828