38.利用matlab解 有约束无约束的参数估计对比（matlab程序）

1.简述

      

1.离散型随机变量的极大似然估计法:
(1) 似然函数
若X为离散型, 似然函数为

(2) 求似然函数L(θ)的最大值点 θ, 则θ就是未知参数的极大似然估计值.
2.连续型随机变量的极大似然估计法:
(1) 似然函数
若 X 为连续型, 似然函数为

(2) 求似然函数L(θ)的最大值点θ, 则θ就是未知参数 的极大似然估计值.

一、矩估计
设总体X的均值、方差均存在，样本（X1，X2，……，X n）,则不管总体服从什么分布，总体均值的矩估计均为样本均值，方差的矩估计均为样本二阶中心矩。
matlab中提供了下列函数来实现总体均值的矩估计值与方差的矩估计值的计算，如下：
mu_ju=mean(X) % 返回样本X的均值
sigma2_ju =moment(X,2) % 返回样本X的2阶中心矩
例：来自某总体X的样本值如下：
232.50, 232.48, 232.15, 232.52, 232.53, 232.30, 232.48, 232.05, 232.45, 232.60, 232.47, 232.30，求X的均值与方差的矩估计。

>> x=[232.50,232.48,232.15,232.52,232.53,232.30,232.48,232.05,232.45,232.60,232.47,232.30]
 mu_ju=mean(x)
sigma2_ju= moment(x,2)
x =
  232.5000  232.4800  232.1500  232.5200  232.5300  232.3000  232.4800  232.0500  232.4500  232.6000  232.4700  232.3000
mu_ju =
  232.4025
sigma2_ju =
    0.0255

二、单个总体极大似然估计与区间估计（参数均未知）
命令: [a,b]=namefit (X, ALPHA) % 返回总体参数的极大似然估计a与置信度为100(1- ALPHA)的置信区间 [a,b]，若参数为多个，ab也是多个，若省略ALPHA，则置信度为0.95
下表列出了几种常用分布的参数估计函数：

函数名    调 用 形 式    函 数 说 明
binofit    PHAT= binofit(X, N)；[PHAT, PCI] = binofit(X,N)；[PHAT, PCI]= binofit (X, N, ALPHA)    二项分布的概率的最大似然估计置信度为95%的最大似然估计和置信区间返回水平α的参数估计和置信区间
poissfit    Lambdahat=poissfit(X)；[Lambdahat, Lambdaci] = poissfit(X)；[Lambdahat, Lambdaci]= poissfit (X, ALPHA)    泊松分布的参数的最大似然估计置信度为95%的最大似然估计和置信区间返回水平α的λ参数和置信区间
normfit    [muhat,sigmahat,muci,sigmaci] = normfit(X)；[muhat,sigmahat,muci,sigmaci] = normfit(X, ALPHA)    正态分布的最大似然估计，置信度为95%返回水平α的期望、方差值和置信区间
betafit    PHAT =betafit (X)；[PHAT, PCI]= betafit (X, ALPHA)    返回β分布参数a和 b的最大似然估计返回最大似然估计值和水平α的置信区间
unifit    [ahat,bhat] = unifit(X)；[ahat,bhat,ACI,BCI] = unifit(X)；[ahat,bhat,ACI,BCI]=unifit(X, ALPHA)    均匀分布参数的最大似然估计置信度为95%的参数估计和置信区间返回水平α的参数估计和置信区间
expfit    muhat =expfit(X)；[muhat,muci] = expfit(X)；[muhat,muci] = expfit(X,alpha)    指数分布参数的最大似然估计置信度为95%的参数估计和置信区间返回水平α的参数估计和置信区间
gamfit    phat =gamfit(X)；[phat,pci] = gamfit(X)；[phat,pci] = gamfit(X,alpha)    γ分布参数的最大似然估计置信度为95%的参数估计和置信区间返回最大似然估计值和水平α的置信区间
weibfit    phat = weibfit(X)；[phat,pci] = weibfit(X)；[phat,pci] = weibfit(X,alpha)    韦伯分布参数的最大似然估计置信度为95%的参数估计和置信区间返回水平α的参数估计及其区间估计
Mle    phat = mle(‘dist’,data)；[phat,pci] = mle(‘dist’,data)；[phat,pci] = mle(‘dist’,data,alpha)；[phat,pci] = mle(‘dist’,data,alpha,p1)    分布函数名为dist的最大似然估计置信度为95%的参数估计和置信区间返回水平α的最大似然估计值和置信区间仅用于二项分布，pl为试验总次数
对于上表函数，以均匀分布的参数估计为例说明：
例、产生 100 行1列服从区间(1, 10)上的均匀分布的随机数, 计算区间端点“a”和“b”的极大似然估计值, 求出置信度为0.95的这两个参数的置信区间。

三、单个正态总体参数估计
设X1,X2,……Xn,为来自正态总体N(u，sigma^2)的一个样本，求u，sigma ^2的极大似然估计与区间估计。此处会较多用到数理统计的知识，参见下表。
只需要按照表中所对应的估计函数计算即可。只不过在计算之前，有必要去了解一下matlab中如何去计算某一分布的临界值。

2.代码

function my_ga_fit_nonline()
x=[8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00 14.00 14.00 16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.00 20.00 22.00 22.00 24.00 24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00 34.00 36.00 36.00 38.00 38.00 40.00 42.00]';

y=[0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43 0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41 0.40 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39]';

A=[];
B=[];
AE=[];
BE=[];

 lb =[ ]; 
 ub =[ ];        
 options=optimset('display','iter','MaxFunEvals',1e3,'MaxIter',1e10,'TolFun',1e-8,'TolX',1e-8);
numberOfVariables = 2;%需拟合参数的个数

opt= gaoptimset('PlotFcns',{@gaplotbestf,@gaplotmaxconstr},'Display','iter',...
    'PopulationSize',50, 'Generations',50,'StallGenLimit',50,'StallTimeLimit',50,'CrossoverFraction',0.85,'MutationFcn',@mutationadaptfeasible);
[x_unc, fval_unc, exitflag]=ga(@myfuntest,numberOfVariables,A,B,AE,BE,lb,ub,[]);
[x_ga_fit, fval_error, exitflag]=ga(@myfuntest,numberOfVariables,A,B,AE,BE,lb,ub,@myconfun,opt); 

x_unc
x_ga_fit
y_48=(x_ga_fit(1)+(0.49-x_ga_fit(1)).*exp(-x_ga_fit(2).*(10-8)))

y_unc=x_ga_fit(1)+(0.49-x_ga_fit(1)).*exp(-x_ga_fit(2).*(x-8));
 yy=x_unc(1)+(0.49-x_unc(1)).*exp(-x_unc(2).*(x-8));
 err_unc=sum(abs(y_unc-y))
 err=sum(abs(yy-y))
figure
 plot(x,y,'r*',x,yy,'b-^',x,y_unc,'--gs')
 
function f = myfuntest(z)
  f=zeros(1,1);
for i=1:length(x)
  f=f+((z(1)+(0.49-z(1)).*exp(-z(2).*(x(i)-8))-y(i)).^2);

end

end

function [c, ceq] =myconfun(z)
        c=[];
        ceq=[(z(1)+(0.49-z(1)).*exp(-z(2).*(10-8))-0.48)];
end
end
 

3.运行结果