$Matlab入门课程, 数学计算
Jul 202011
MLE(maximum likelihood estimation，最大似然估计)，的基本原理是通过选择参数使似然函数最大化。此处，我们只讲随机数分布估计。 假设随机分布的PDF为 f(x,theta), 其中x为随机数，theta是分布的参数。我们有x1,x2,…,xn共n个取样观察值，似然函数
L(theta)=f(x1,theta)*f(x2,theta)*…*f(xn,theta)
然后选择能让L()最大的theta取值。注意，这个不叫做”极大似然估计”，因为theta可能是有上限或下限，而能够让L()最大的theta值可能恰好是上限或下限，此时的theta就不是极大值点。
下面我们用随机生成的数据实地讲述mle()函数的用法。
语法一：mle(观察值, ‘distribution’, ‘分布名称’), 下面两个命令，第一个是生成100万个服从标准正太分布的随机数，第二句是用mle()函数估计这些随机数服从的正太分布的参数值。
》 testdata=randn纤维蛋白酶(1e6,1);
头脑风暴2011》 [paramhat,paramint]=mle(testdata,'distribution','norm')
paramhat =
-0.000218972940353 0.999431366252859
paramint =
-0.002177825773582 0.998048679206182
0.001739879892876 1.000818918426322
可以获得两个结果，paramhat和paramint，paramhat中有两个数值，第一个是正态分布均值（mu）的点估计，后一个是正态分布的标准差（sigma）。paramint中第1，2两列分布对应mu和sigma的区间估计。从结果来看，mu和sigma的区间估计都很窄，而且很接近我们生成随机数时使用的参数值。 这个语法支持的分布类型列表如下（Matlab2011A):
‘beta’ Beta
‘bernoulli’ Bernoulli
‘binomial’ Binomial
‘birnbaumsaunders’ Birnbaum-Saunders
‘discrete uniform’ or ‘unid’ Discrete uniform
‘exponential’ Exponential
‘extreme value’ or ‘ev’ Extreme value
‘gamma’ Gamma
‘generalized extreme value’ ‘gev’ Generalized extreme value
‘generalized pareto’ or ‘gp’ Generalized Pareto
‘geometric’ Geometric
‘inversegaussian’ Inverse Gaussian
‘logistic’ Logistic
‘loglogistic’ Log-logistic
‘lognormal’ Lognormal
‘nakagami’ Nakagami
‘negative binomial’ or ‘nbin’ Negative binomial
‘normal’ Normal
‘poisson’ Poisson
‘rayleigh’ Rayleigh
‘rician’ Rician
‘tlocationscale’ t location-scale
‘uniform’ Uniform
‘weibull’ or ‘wbl’ Weibull
语法二：mle(观察值, ‘pdf’, 自定义分布pdf, ‘start’, 猜测的分布参数值), pdf是分布的概率密度函数，格式是f(X,theta)， 前面是X，后面跟参数值。
也支持cdf等等，详见帮助。这个语法适合不在上述列表内的随机分布。我们先看一个手工指定正态分布PDF函数的例子。下面的代码：第一句生产10万个标准正态分布随机数，第二句自定义一个正态分布PDF函数，第三句用MLE()
》 testdata=randn(1e5,1);
》 mynormpdf=@(x,mu,sigma)(1/sqrt(2*pi*sigma*sigma)*exp(-(x-mu).^2/2/sigma/sigma));
》 [paramhat,paramint]=mle(testdata,'pdf',mynormpdf,'start',[.1,.5])
paramhat =
0.002066365428022 1.002194631105154
paramint =
-0.004145187769578 0.997802401528087
0.008277918625623 1.006586860682221
结果与第一种语法的差不多。下面看一个实例。数据来自网友提供的问题：/thread-1137305-1-1.html
我们先看这个数据的直方图： hist(testdata,20)
所以觉得可能是混合正态分布，那么我们要估计5个参数(mu1,mu2,sigma1,sigma2,rho)。下面的代码首先定义好这个混合正态分布的PDF函数
霍纳mixedpdf=@(x,mu1,mu2,s1,s2,rho)(rho*normpdf(x,mu1,s1)+(1一切如新-rho)*normpdf(x,mu2,s2));
注：此处normpdf是正态分布的PDF函数。
然后我们先要猜测这5个参数的初始值。猜测办法一是看图，二是凭直觉。从图上我们看到mu1, mu2 一个是-0.7左右，一个是0左右，至于其他三个值我们分别选.15,.15,.5这种比较中庸的值。然后就用mle()函数！
》 [phat1,pci1]=mle(testdata,'pdf',mixedpdf,'start',[-.7,0,.15,.15,.5])
Warning: Maximum likelihood estimation did not converge. Iteration limit exceeded.
> In stats\private\mlecustom at 245
In mle at 232
phat1 =
Columns 守正创新1 through 4
-0.631157726178361 0.040647319628326 0.113402249487918 .0934********
Column 5
0.300855540018213
pci1 =
Columns 1 through 4
-0.664050515555056 .023*********552 .0834******** .0798********
-0.598264936801666 .0574******** 0.143334000192095 0.107071343948353
权益乘数
Column 5
0.234297838028327
0.367413242008099
结果说MLE没有完全收敛，我们可以用得到的结果（phat1）作为新的猜测值代入mle()函数中再做一遍。
》 [phat2,pci2]=mle(testdata,'pdf',mixedpdf,'start',phat1)
phat2 =
Columns 1 through 4
-0.631515418563246 0.040992492767408 0.113519235125114 0.092941280730321
Column 5
0.301236890533836
pci2 =
Columns 1 through 4
-0.664590758235432 .024********* 0.083097143105910 .0793********
-0.598440078891060 .0577******** 0.143941327144319 0.106538716706501
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议