一个含定积分非线性方程组的数值解与统计手册一处错误.doc

一个含定积分非线性方程组的数值解与《统计手册》一处错误 为寻求彩电平均寿命的贝叶斯估计需要确定的先验分布。由经验选择倒伽玛分布作为的先验分布。由大量寿命试验数据得到如下含有定积分的非线性方程组 （1） 其中，为倒伽玛分布的密度函数 《统计手册》给出的计算结果为，可惜，这个结果显然是错误的，因为此时，与第一个方程矛盾。下面用MATLAB给出正确结果。 1直接法 对于非线性方程组，可以直接用fsolve函数求其数值解。fsolve函数的调用格式为：x = fsolve(fun, x0, options)，其中x为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，x0是求根过程的初值，options为最优化工具箱的选项设定。 (1) 建立函数文件pai.m。 由于非线性方程组内有函数，先建立函数文件 %函数 function fun = pai(theta,alpha,beta) if theta0 fun = beta^alpha./gamma(alpha).*1./(theta.^(alpha+1)).*exp(-beta./theta); else fun = 0; end (2)建立非线性方程组文件myequs.m。 准备了非线性方程组可能用的文件后，建立非线性方程组文件 function F = myequs(x) alpha = x(1); beta = x(2); F = [beta/(alpha-1) - 30000; quadl(@(theta)pai(theta,alpha,beta),0,11250)-0.1]; % quadl为数值积分命令 (1) 在给定的初值，下，调用fsolve函数求方程的根。 clear,clc format long g x0 = [3;2000]; % 初值 xsol = fsolve(myequs,x0,optimset(Display,off)) alpha = xsol(1); beta = xsol(2); format long g x0 = [3, 2000]; % 初值 xsol = fsolve(myequs,x0,optimset(Display,off)) % alpha = xsol(1); beta = xsol(2); 运行结果为 xsol = 2.99138028363385 59741.4085090153 即，，代入原来方程 errors=myequs(xsol); 得到errors = -3.27418092638254e-011 4.16333634234434e-017 误差平方和 mis1=errors*errors mis1 = 1.07202607387046e-021 误差非常小，可见得到了较高精度的数值解。 2 化为无约束优化方法 设非线性方程组为，即 (2) 为了应用无约束方法求解上面问题(2)，通常将其化为最小值问题 (3) 或 (4) MATLAB中有许多优化命令可以求解，如用于无约束优化的fminunc，无约束不可微优化的fminsearch。它们的调用格式为：[x, fval] = fminunc(fun,x0,options)或[x, fval] = fminsearch(fun,x0,options)，其中x为返回的解，fval为返回解x对应的极小值，fun是用于定义需求解的多元函数文件名，x0是求解过程的初值，options为最优化工具箱的选项设定。 最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中‘off’为不显示，‘iter’表示每步都显示，‘final’只显示最终结果。optimset(‘Display’, ‘off’)将设定Display选项为‘off’。 一般地，fminunc速度更快，但不稳定，fminsearch速度稍慢，但对初值不如fminunc敏感。我们这里使用fminsearch。 按式(3)建立与非线性方程组等价的多元函数mymfun.m % 非线性方程组平方和 function f = mymfun(x) alpha = x(1); beta = x(2); %不一定是整数 F = [beta/(alpha-1)