如何用matlab编最大似然法的程序

命令

利用mle函数进行参数估计

函数

mle

格式

phat=mle('dist′

,X)%返回用dist指定分布的最大似然估计值

[phat,

pci]=mle('dist′

,X)

%置信度为95

[phat,

pci]=mle('dist′

,X,alpha)%置信度由alpha确定

[phat,

pci]=mle('dist′

,X,a;pha,pl)%仅用于二项分布，pl为试验次数。

说明

dist为分布函数名，如：beta(β分布)、bino（二项分布）等，X为数据样本，alpha为显著水平α，1-a为置信度。

phat为估计值，pci为置信区间

希望能对你有所帮助！

数值积分法是求定积分的近似值的数值方法。即用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。

求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解。对微积分学作出杰出贡献的数学大师，如I牛顿、L欧拉、CF高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础。

数值积分法也是计算机仿真中常用的一种方法。在已知函数的微分方程时，求解函数下一时刻的值，我们主要有欧拉法、梯形法和龙格库塔法。

欧拉法，这些方法中精度最低的，程序相对简单。欧拉法的表达式可以写成下面的形式：

我们用欧拉法近似替代则有：

y'=f(t,y)

其中y’(k) 函数y(x)在k时刻的导数，h为积分的步长(也可以说是采样周期)。

欧拉法的主要思想是用当前时刻的值y(k)  + 当前时刻  y’(k)h（积分步长）来近似替代y(k+1)时刻的值。这种方法主要的误差来源于，在h区间内y’(t)是变化的，而这里我们都用y(k)替代了。

梯形法，是在欧拉法的基础上进行改进的算法，也称为改进的欧拉法。

在用欧拉法求出y(k+1)之后，我们再求出y(k+1)处的导数f (k+1,y(k+1) )。再用求梯形面积的方法求解这个积分的大小，得到上面的公式。

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。这里主要介绍工程中常用的四阶龙格-库塔法。

该算法是构建在数学支持的基础之上的。对于一阶精度的拉格朗日中值定理有：

对于微分方程：

y'=f(x,y)

y(i+1)=y(i)+hK1

K1=f(xi,yi)

当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K的近似值，那么就会得到二阶精度的改进拉格朗日中值定理：

y(i+1)=y(i)+[h( K1+ K2)/2]

K1=f(xi,yi)

K2=f(x(i)+h,y(i)+hK1)

依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：

y(i+1)=y(i)+h( K1+ 2K2 +2K3+ K4)/6

K1=f(x(i),y(i))

K2=f(x(i)+h/2,y(i)+hK1/2)

K3=f(x(i)+h/2,y(i)+hK2/2)

K4=f(x(i)+h,y(i)+hK3)

k1是时间段开始时的斜率；

k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；

k3也是中点的斜率，但是这次采用斜率k2决定y值；

k4是时间段终点的斜率，其y值用k3决定。

上面三种方法，龙格-库塔法精度最高，改进的梯形法次之，欧拉法精度最低。对于精度，我们引入截断误差的概念。

所谓截断误差就是，近似解和它真值之间的差值。数值积分对的真值做泰勒展开。而欧拉法在数值积分时，取的是泰勒级数的第一项，梯形法取的是泰勒级数的前两项，四阶龙格-库塔法取的是泰勒级数前三项。所以他们的截断误差是泰勒级数剩下的部分。

数值积分法在控制领域主要应用于对状态空间的输出进行求解。下面是我做过的一个matlab仿真

已知开环传递函数

5 s + 100

-------------------------------------------

s^4 +8 s^3 + 3199 s^2 + 7521 s

给定单位阶跃，求解闭环函数的输出曲线。

首先我们要把传递函数转换成状态方程。在

主M文件：

clear

t=15;h=01;

num = [0 0 0 5 100];

den = [1 8 3199 7521 0];

den = den + num;

Tsys1 =tf(num,den);

[A,B,C,D] =tf2ss(num,den);

x0 = [0 0 0 0]';

u = 1;

i=0:h:t+h;

Y0 = EEluer(h,t,A,B,C,D,x0,u);

Y1 = SecondRunge(h,t,A,B,C,D,x0,u);

Y2 = FourRunge(h,t,A,B,C,D,x0,u);

plot(i,Y0,i,Y1,'-',i,Y2,'g-');

legend('euler','SecondRunge','FourRunge');

四阶龙格库塔法M  function：

function save=FourRunge(h,t,A,B,C,D,x0,u)

count = 1;

y0 =C x0 + D u;

len = size(y0,1);

wide=length(0:h:t);

save=zeros(wide,len);

save(count,len) = y0;

for i=0:h:t

k1 = h My_function(A,B,x0,u) ;

k2 = h My_function(A,B,x0+05k1,u);

k3 = h My_function(A,B,x0+05k2,u);

k4 = h My_function(A,B,x0+k3,u);

x1 = x0 + (1/6)(k1 + 2k2 +2k3 +k4);

y0 = C x1 + D u;

count = count + 1;

save(count,len) = y0;

x0 = x1;

end

end

欧拉法 M function ：

function save=EEluer(h,t,A,B,C,D,x0,u)

count = 1;

y0 =C x0 + D u;

len = size(y0,1);

wide=length(0:h:t);

save=zeros(wide,len);

save(count,len) = y0;

for i=0:h:t

x1 = h My_function(A,B,x0,u) +x0;

y0 = C x1 + D u;

count = count + 1;

save(count,len) = y0;

x0 = x1;

end

end

梯形法  M function：

function save=SecondRunge(h,t,A,B,C,D,x0,u)

count = 1;

y0 =C x0 + D u;

len = size(y0,1);

wide=length(0:h:t);

save=zeros(wide,len);

save(count,len) = y0;

for i=0:h:t

k1 = h My_function(A,B,x0,u) ;

k2 = h My_function(A,B,x0+k1,u);

x1 = x0 + (1/2)(k1 + k2);

y0 = C x1 + D u;

count = count + 1;

save(count,len) = y0;

x0 = x1;

end

end

微分方程 M function

function [ yy ] = My_function(A,B,x0,u)

yy = Ax0 +  Bu;

end

实验结果

图中振荡最大的是欧拉法，梯形法和四阶龙格库塔法精度差不多。四阶龙格库塔法精度最高。

x=[10,20,30,35,40,50,60,70,80,90,100]';

y=[21681,22030,22482,22783,23075,23687,24364,25053,25882,26663,27611]';

X=[ones(size(x)) x x^2];

coe=X\y

MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

勉强弄了3种，你试试，不行再说。

clear all;clc;

%---------------------------------

s1=sum((1/(1:1000))^2);

%---------------------------------

n=1000;

s2=0;

for i=1:n

a(i)=1/i^2;

s2=s2+a(i);

end

%---------------------------------

s3=(1/(1:1000))(1/(1:1000))';

以上就是关于如何用matlab编最大似然法的程序全部的内容，包括:如何用matlab编最大似然法的程序、工程数值方法的MATLAB程序编写、如何用MATLAB编写程序（最小二乘法)曲线拟合等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！