在做论文的最后阶段需要用到，我的具体程序如下：

1：
function  y=jie_cheng(n)
  y=1;
  for j=1:1:n;
    y=y*j;
  end;
return

2：
function mt(a0,a,b)
syms a0 b a
ff=((1+b*t)^(i+1)-1)/((i+1)*jie_cheng(i+1));
Y=symsum(ff,i,0,inf);
mt=a0*(1+a*t-(b*t+1)/exp(b*t))-(a0*a*(1+b*t))/(b*exp(b*t+1))*(log(b*t+1)+Y)
3：
%syms b a
%f=((1+b*0.5)^(i+1)-1)/((i+1)*jie_cheng(i+1));
%symsum(f,i,0,inf);
%mt=a0*(1+a*t-(b*t+1)/exp(b*t))-(a0*a*(1+b*t))/(b*exp(b*t+1))*(log(b*t+1)+symsum(f,i,0,inf));
syms a0 b a
ff=((1+b*t)^(i+1)-1)/((i+1)*jie_cheng(i+1));
Y=symsum(ff,i,0,inf);
mt=a0*(1+a*t-(b*t+1)/exp(b*t))-(a0*a*(1+b*t))/(b*exp(b*t+1))*(log(b*t+1)+Y);
tt=[5.5 7.33 10.08 80.97 84.91 99.89 103.36 113.32 124.71 144.59 152.40 166.99 178.41 197.35 262.65 262.69 388.36 471.05 471.50 503.11 632.42 680.02];
yy=subs(mt,t,tt);
yk=[1:1:22];
aaa=yk(1)*log(yy(1))-yy(1);

%aaa=(yk(k)-yk(k-1))*log(yy(k)-yy(k-1))-yy(k)-yy(k-1);
%bbb=log(jie_cheng(y[k]-y[k-1]))
bbb=jie_cheng(yk(1));
lnL=0;
for j=2:1:22
  lnL1=(yk(j)-yk(j-1))*log(yy(j)-yy(j-1))-yy(j)+yy(j-1)-jie_cheng(yk(j)-yk(j-1));
  lnL=lnL+lnL1;
end
lnL=lnL+aaa-bbb;
ddnn1=diff(lnL,a0,1);
ddnn2=diff(lnL,a,1);
ddnn3=diff(lnL,b,1);
 
最后四步分别是最大似然函数lnL求关于a0,a,b的偏导数，其中约束条件为：
   0<=a0<=1024
    0<a<1
    0<b<1

因为对最大似然函数lnL求关于a0,a,b的偏导数用人工手算的方法无法表达出来，所以只有用计算机运算了， 具体用牛顿法解出a0,b,a的算法本人实在不知道怎末写了，请教大家帮助一下！！