参数估计软件（精选十篇）
参数估计软件 篇1
1972年, Z.Jelinski和P.Moranda 提出的JM模型[4]是软件可靠性模型的代表, 它是面向时间的模型, 属于出错计数模型, 既可以估计软件的可靠性, 又可以给出程序中残留错误的数目, 但是它的假设排错是完全的, 不能适应由于排错不完全引起的新错的情况, 所以, 这并不符合实际情况。1986年出现多故障模型, 它允许在同一时刻, 查出、改正或引入的错误个数可以多于一个[5]。但这些假设与处理还没有完全说明实际情况, 在改错过程当中, 引入错误个数与剩余错误个数并不一定存在正比关系, 改正不同的错误会对软件产生不同的影响。所以, 有必要针对故障率与剩余错误个数成其他关系的情形进行研究。
对故障率与剩余错误个数成指数关系的情形进行了研究, 在排错不完全的假设下对软件的初始错误个数和指数系数给出了一个合理的估计方法。
1 模型的建立以及参数估计
1.1 模型的假设
1) 软件中初始错误个数为一个未知但固定的常数, 用X表示;
2) 将错误引入数看成负的消除数, 设软件错误消除数服从参数为θ的指数分布;
3) 设μ表示第i次改正错误个数的数学期望, σ表示改正错误个数的方差。
其中:
μ>0表明测试过程是顺利进行的, 软件总错误个数在减少, 软件可靠性随之提高。
μ=0表明测试工作没有收效, 软件总错误个数基本不变, 软件可靠性基本不变。
μ<0表明测试工作是失败的, 软件总错误个数在增加, 软件可靠性随之降低。
σ值的大小表示测试工作过程中改正错误数波动的大小[7]。
μ, σ的值可以由测试记录数据估计得出。设ai是每次改正的消除数, 则
undefined
4) 错误一旦查出, 立即改正错误。分三种情况考虑, 这三种情况发生的概率相同, 服从离散型均匀分布。
情况一:发现错误就能顺利改正错误, 于是每次改正之后, X就要减去μ;
情况二:发现错误改正错误的同时又产生一个错误, 于是每次改正之后, X就要加上μ;
情况三:发现错误改正错误的同时又产生多个错误, 于是每次改正之后, X就要加上 (μ+a) 。
5) 在任何时候的故障率都与软件中的剩余错误个数成指数关系, 指数系数用θ表示。
1.2 模型的建立
根据假设, 在第一个错误被排除之后, 故障率会发生变化。
情况一:故障率由undefined变为undefined;
情况二:故障率由undefined变为undefined;
情况三:故障率由undefined变为undefined。
设t1, t2, …, tn表示相继出现的并被改正的错误之间的时间区间样本。根据假设, 在每两个错误之间只有唯一的故障率, 关于ti的密度即为:
1.3 参数的估计
模型中含有两个未知参数X和θ, 所以, 对模型的参数估计就是X和θ的值。
由 (1) 式可得到对数似然函数为:
undefined
令undefined, 对 (2) 式中的X和θ求偏导, 并令结果为0, 得
undefined
undefined (4)
方程 (3) 式中都不含θ, 而由测试收集的数据, 可计算出undefined与undefined的值。将它们代入 (3) 式中, 则可以解出X的估计值undefined:
undefined
再将X代入 (4) 式中, 可以解出θ的估计值undefined:
undefined
2 模型评价
模型 (下称优化模型) 在以下几个方面有所改进:
1) 软件所改正的错误数与软件错误数的改正的次数两个概念完全分开, 并且软件所改正的错误数大于软件错误数的改正的次数。
2) 在任何情况下, 改正错误时既可能减少错误, 也可能增加错误。
以纵轴表示改错过程中软件剩余的错误个数, 横轴表示时间, 某次测试的优化模型剩余的错误个数与时间关系, 如图1所示。
3) 从图1可以看出, 优化模型的剩余错误量都是在减少的, 但改正错误并不一定会使错误总数减少, 反而可能出现错误数增多的现象, 这与实际情况是比较相吻合的。
4) 当时, 它就是原来的JM模型。
在研究中, 仅考虑了三种情况为离散型均匀分布, 即等可能分布, 在以后的研究中, 可以进一步研究不是等可能分布的情况。
摘要：20世纪70年代提出的JM模型是最早代表软件可靠性的数学模型, 其假设排错是完全的, 这与实际情况不相符合。因此, 有必要对JM模型假设进行改进, 从而提高软件可靠性。在排错是不完全的假设下建立了一个软件可靠性模型, 然后运用最大似然估计法进行参数估计, 找出了一个比较合理的估计初始错误个数和指数系数的方法。
关键词：软件可靠性,JM模型,参数估计,故障率
参考文献
[1]尚珊珊, 赵轶群.软件可靠性综述[J].软件导刊, 2006 (8) :3-5.
[2]邹丰忠, 李传湘.软件可靠性模型理论分析[J].武汉水利电力大学学报, 1998, 31 (2) :76-80.
[3]朱东鸣.MATHEMATICA软件建立软件可靠性模型[J].成都大学学报 (自然科学版) , 1996, 15 (2) :32-37.
[4]徐仁佐, 谢云, 郑人杰.软件可靠性模型及应用[M].北京:清华大学出版社, 1994.
[5]谭民, 疏松桂.软件可靠性模型的发展[J].计算机学报, 1990 (5) :382-390.
[6]王峰, 陈杰, 喻林.计算机软件测试[M].北京:机械工业出版社, 2004.
考研数学压轴必考题型：参数估计 篇2
参数估计是考研概率的最后一个考点，近几年参数估计一直是数一和数三的必考题目，必出现在整张试卷的最后一道大题，压轴出场，分值11分。
虽然考研数学一和数学三最后一道题均未考查，但16年数学一填空题考查了区间估计，分值4分，但数一和数三均考查了一道大题，分值11分，迄今参数估计这个考点的重要地位仍不可撼动。跨考教育数学教研室田晓辉老师来为大家解析。
参数估计这章，数一和数三公共考点为点估计，包括矩估计和极大似然估计，另外数一还考查区间估计，包括单个正态总体的均值和方差的区间估计、两个正态总体的均值差和方差比的区间估计。
本章考研主要题型为：
(1)参数的点估计：矩估计、极大似然估计估计量的评选标准(数一考查)
(2)参数的区间估计：正态总体的区间估计(数一考查)
矩估计的基本思想：由大数定律可知样本矩、样本矩的连续函数依概率收敛于相应的总体矩、总体矩的连续函数，由此可建立总体分布中未知参数满足的方程(组)，解之可得总体未知参数的点估计。这种构造点估计量的方法称为矩估计法，求得的点估计称为矩估计量(值)其方法步骤如下：
构建未知参数的方程，通过总体的`原点矩来构造
解方程，解出未知参数
用样本矩代替总体矩，得未知参数的矩估计量(值)
极大似然估计法的基本思想：样本发生的可能性最大原则――即对未知参数进行估计时，在未知参数的变化范围内选取使“样本取此观测值”的概率最大的参数值作为未知参数的点估计。这样得到的矩估计值为最大似然估计值，相应的量为最大似然估计量。其方法步骤为：“造似然”求导数，找驻点得估计。
构造自然函数，注意，离散总体和连续总体的似然函数不同
取对数
求导数找驻点得估计。
注意，若似然方程无解，则必有导数大于或小于零，此时只要在未知参数的变化范围内找其右边界点或左边界点即可。
估计量的评选标准：无偏性、有效性、一致性，掌握其概念即可。无偏估计考查较多。
参数的区间估计：了解区间估计概念、掌握求置信区间的方法。求置信区间的一般方法步骤为：
第一步，选枢轴量定分布;
第二步，造大概率事件得不等式;
第三步，解不等式得置信区间。
以上是数一和数三对参数估计部分的全部考点，期望大家能熟练理解其思想和熟练掌握方法步骤，多练习，已达到熟练解题的要求。
参数估计软件 篇3
一、项目反应理论模型发展
用来描述IRT模型有很多种。本文主要讨论单维二值反应模型。
（一）正态拱形曲线模型
1952年，Lord提出了双参数正态拱形曲线模型 ( two-parameter normal ogive model)，其表达式为：Pt(H)=ai(H-bi)edz。其中Pi(H)表示能力为H的被试者在项目i上的正确作答概率。项目特征曲线为 S形， 因为它形同累积分布函数且是中心对称的，自然假设曲线可用正态分布的形式来描绘它。bi与能力定义在同一量表上，为难度参数，也称为曲线的定位参数；而ai则被称为区分度参数，是曲线拐点处的切线斜率的函数，ai越大，曲线越陡峭，反之亦然。该模型对于帮助理解项目反应理论模型的特性以及模型参数的均有重要意义。但由于采用积分函数形式，不利于进一步的分析计算。
（二）Rasch模型
Rasch模型是丹麦学者 Rasch于 1960年提出的测量模型。该模型的特点是完全根据被试的能力水平和项目难度的关系而推导出的正确作答概率公式，其表达式为：Pi(H)=。式中符号代表的意义同正态拱形曲线模型相同。Rasch模型中只有一个项目参数，因而该模型又称作单参数模型，简称 1- PL模型。
（三）Logistic模型
针对被试者在客观多项选择题和是非判断题上可能出现的猜测现象， Birnbaum 增加了猜测参数c，并于 1957年到1958年将 Lord的双参数正态拱形曲线模型改为 Logistic模型，其表达式为：Pi(H)=ci+。其中D为1.7，参数ai、bi与正态拱形曲线模型中的定义相同。这就是三参数 Logistic模型 [简称3-PL模型 ]。易知，若答题过程中没有猜测因素，即c=0时，该模型则为双参数 Logistic模型 [简称 2-PL模型 ]；若双参数 Logistic模型中ai=1，该模型则为单参数 Logisti c模型，其形式同Rasch模型。
二、参数估计方法研究
应用IRT模型对项目的不同参数进行估计是连接项目反应理论与应用的最关键的环节。所谓参数估计是指根据被试者的作答反应矩阵估计出被试的能力参数和每个项目的参数。
（一）极大似然估计法
极大似然估计法是根据被试者的作答反应矩阵，在局部独立性的条件下，导出参数估计的似然函数，然后通过求取似然函数的极大值，估计项目参数和被试能力参数。当参数变化时，概率的非极大值可能不止一个，但极大值一般只有一个。极大似然估计法又可细分为如下三种：
(1)联合极大似然估计法 ( JMLE )
在测量实践中，一般的情况是既不知道被试者的能力H，也不知道项目参数，因此只能同时对这两个参数进行联合极大似然估计。该方法是 Birnbaum首先采用的，公式为：L(u1,u2,u3…,un|H) =。其中uij是考生 i对反应模式为 ( 1,2 ,3,….,n)的第 j题的反应，pij表示能力为Hi的被试者答对第 j题的概率，Q1-表示能力为Hi的被试答错第j题的概率。
(2)条件极大似然估计法 (CMLE)
该方法的前提是得到能力参数的充分统计量。因为 Rasch模型的实得分数是被试能力参数的充分统计量，所以对于 Rasch模型可以采用条件极大似然法进行参数估计。但是对于二值计分的双参数模型和三参数模型以及多值计分模型都不能运用这一种方法。
(3)边际极大似然估计法 (MMLE)
对 Rasch模型来说，条件极大似然估计的效果基本上和边际极大似然估计相近。但边际极大似然估计的一个最大缺点是运算量太大，需要进行大量的积分运算，因而只要项目数稍多，则无法实现。
综上，虽然极大似然估计具有一致性、渐进正态性和有效性等基本性质，成为一种应用最为广泛的参数估计方法，但它也有两条明显的缺点：(1)没有利用被试者能力的先验知识； (2)对于满分和零分的被试者无法进行参数估计。
(二)贝叶斯估计法
贝叶斯估计方法是指利用贝叶斯原理，确定项目参数和被试者能力参数的先验分布，建立联合极大似然函数，然后通过求取联合极大似然函数的极大值，估计出项目参数和被试者能力参数。对被试者相应能力的先验概率分布同测验使用的时间长短是正向关系，有可能做出较为客观的估计，但是对项目参数先验概率分布的估计则纯粹是主观的。
三、结语
各参数的先验概率分布确定之后，贝叶斯估计和极大似然估计的方法大致相同。贝叶斯估计方法运用于二值计分模型也是可行的。对于等级计分模型，由于每题有多个难度值，而且这多个难度值有可能是逐渐增大，也有可能是没有变化规律的，项目参数的先验分布难以确定。因此贝叶斯估计方法运用到模型参数当中，仍有许多理论上和技术上的问题未解决，而且关于估计方法的稳定性还缺乏证据。
参数估计软件 篇4
马氏链蒙特卡洛方法( Markov Chain Monte Carlo,MCMC) 提供了另一种计算可靠性参数的分层贝叶斯方法[1,4],其利用计算机模拟从后验分布抽样然后用样本估计值推断总体特征。
1 马氏链蒙特卡罗方法( MCMC) 及Gibbs抽样
在贝叶斯统计分析中,统计模型可概括为先验分布 π( θ) 与样本分布p( x | ( θ) ,其中x = ( x1,…,xn) 为容量为n的样本,θ = ( θ1,…,θp) 为参数,样本和参数都是随机的。由贝叶斯定理,更新后的后验分布表示为[4]
因此,形式上参数 θ 的贝叶斯估计为
更一般参数的函数g( θ) 的贝叶斯估计为
可以论证式( 3) 能够采用下面的平均值近似求得
式( 4) 中,θ( 1),…,θ( m)为来自后验分布 π(θ|x) 的容量为m的样本。如果此样本是独立的,且容量m足够大,由大数定理,样本均值将收敛到期望值g( θ) 。此估计就是蒙特卡罗估计。实际中,从 π(θ | x) 中抽取独立样本比较困难,通常采用马氏链方法获得 π(θ|x)的非独立样本,由于马氏链具有马氏性、收敛性等好的性质,与从 π(θ|x)抽取独立样本作用一样。这种马氏链抽样方法称为马氏链蒙特卡罗方法。马氏链生成链{ θ( 0),θ( 1),θ( 2)…}中,θ( t + 1)只依赖 θ( t),而不依赖{ θ( 0),θ( 1),θ( 2)…,θ( t - 1)} 。
Gibbs抽样是一种特殊的马氏链蒙特卡洛算法,具体是从某个初始点出发,通过全条件分布的循环抽样产生马氏链,算法如下:
1)给定参数的任意初始值θ(0)1,…,θ(0)p;
2)对t=0,1,2,…进行下面的迭代更新:
a)从分布π(θ1|θ(t)2,…,θ(t)p,x)中产生θ(t+1)1;
b)从分布π(θ2|θ(t+1)1,θ(t)3…,θ(t)p,x)中产生θ(t+1)2;
…
c)从分布π(θp|θ(t+1)1,θ(t+1)2…,θ(t+1)p-1,x)中产生θ(t+1)p;
3)迭代至t+1=n,当n→∞,(θ(t)1,…,θ(t)p)的联合分布将收敛到后验分布π(θ1,…,θp|x)。
2 Win BUGS软件
BUGS ( Bayesian Inference Using Gibbs Sampling) 是采用MCMC算法解决贝叶斯估计的软件,其先由英国MRC生物统计协会研发,后由英国圣玛丽皇家医学院协作开发并提供免费下载。WinBUGS是BUGS软件在Windows操作下的版本。
2. 1 软件的下载安装
Win BUGS软件可在http: / / www. mrc-bsu. cam.ac. uk / bugs / 免费下载,然后按如下步骤安装:
1)下载并安装Win BUGS14.exe;
2)下载并安装1.4.3的补丁;
3)下载key并解码( Decode ) 。 若不安装key,运行时迭代次数较少。
2. 2 建模计算的基本步骤
Win BUGS软件提供了窗口图形界面,支持鼠标点击建模,建模计算基本步骤如下[5,6]:
1)逻辑建模: 根据分析对象,在菜单栏Doodle下,建立Doodle模型或编写模型代码,然后Model→Sepcification→check model检验模型;
2)样本数据输入: 在File→New中新建文本,输入样本数据,用Sepcification→load data加载;
3)编译: 用Sepcification→compile编译模型;
4)初始值设置: 在文本中继续设置参数初值,然后Sepcification→load inits加载,并激活gen inits;
5)变量监控: 在Inference → Sample Monitor Tool中的node中,设置监测变量;
6)模型迭代: 在Model菜单中Update Tool中设置迭代次数后,Update数据;
7)结果输出: 迭代完成后,再次选中Sample Monitor Tool,输入监测的节点( 参数) ,可获得节点的轨迹、历史数据、核密度、参数均值/误差/百分位数、直方图、收敛性诊断、自相关等信息。
3 实例计算
3. 1 核电站需求失效可靠性参数估计
核电站设备需求失效概率通常服从参数为 α、β的Beta分布,概率密度函数为
式( 5) 中,p为失效概率,1 - p为成功概率,α、β 为超参数。采用Win BUGS软件计算68 个核电站辅助给水系统电动部分启动失效的参数,样本数据见表1[1],其中x表示失效次数,n表示需求次数。
3. 2 Doodle模型的建立
根据式( 5) ,在Doodle中建立相应的节点模型,如图1,其中假设超参数 α 服从均值和标准差都为1的指数超先验分布,β 服从gamma( 1. 0,0. 035) 超先验分布,p( i) 为第i个核电站辅助给水系统电动部分启动失效概率,n( i) 描述需求次数,x( i) 描述失效次数,i = 1,…,68。
3. 3 样本数据及参数初始值的输入
在文本中输入失效数据及需求次数数据如下:
list( x = c( 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,2,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) ,n = c( 14,9,24,43,13,24,11,26,57,12,15,41,89,66,14,18,36,16,46,30,34,54,5,28,98,24,32,26,23,45,44,11,54,20,18,18,18,12,13,7,12,9,8,16,3,7,28,24,32,13,17,17,30,41,69,87,35,21,24,26,3,6,103,45,38,51,13,8) ,M = 68) ,然后编译,在显示“mode l compiled”后加载超参数的初始值: list( α = 1,β = 1) 。
3. 4 变量监测和结果输出
在Sample Monitor Tool的节点( node) 中设置监测变量 α、β 及p[i]( i = 1,…,68) ,在Update中设置迭代次数,本文设定100 000 次迭代,并剔除前10 000次迭代结果( burn-in) 以减小MCMC方法初始迭代值的影响。迭代完成后,在node中输入alpha、beta,查看stats获得参数估计结果,如表2。
同时可以计算得到68 个核电站的后验失效概率p[i],如图2。部分核电站后验需求失效概率见表3。
3. 5 收敛性判断
Gelman-Rubin方法多用于判断马氏链的收敛性: 通过不同链的充分混合的程度考查链的收敛性———达到平衡后的链的表现是充分混合的,达到平衡后链内方差与链间方差应该是( 近似) 相等的[4]。Win BUGS软件采用Gelman-Rubin法,方差近似用80% 的置信区间代替,计算多条初始值不同的马氏链。软件计算平均链内80% 的区间W( within,蓝线表示) 和所有链数据混合在一起后的链间80% 的区间B ( between,绿线表示) ,并计算R ( 红线) = B /W,马氏链趋于稳定要求[5,7]:
1)链内W的80%区间趋于稳定;
2)链间B的80%区间趋于稳定;
3)比值R趋于1,并设R < 1. 05,表明链内和链间充分混合,马氏链近似收敛。
计算了两组初始值下的马氏链: list( alpha = 1,beta = 1) 和list( α = 10,β = 10) 。
运行完成后,选择bgr diag,生成参数收敛性诊断,见图3。从图上可以看出,在迭代次数大约达到10 000 以后,链内和链间这两条链趋于一致,并且比值R趋于1,因此文中burn-in设置在10 000 次。
4 结语
本文研究了马氏链蒙特卡洛原理及Gibbs抽样原理,简述了Win BUGS软件的下载安装及建模步骤,针对68 个核电站辅助给水系统电动部分启动需求失效数据进行建模计算,计算得到超参数 α 为0. 192 5,β 为46. 77,同时可以获得各个核电站的后验失效概率等有用信息。需要指出的是,Win BUGS
在收敛性判断上,设置burn in的迭代次数上有一定的人为因素,需要考虑比值趋于1 及链间和链内方差都趋于稳定。此外超参数的超先验分布设置上也有人为因素。综上,Win BUGS软件及其信息诊断需进一步深入研究并期望其在国内可靠性参数估计,尤其是稀少样本参数估计中开展应用。
摘要：可靠性数据是核电站概率安全评价的基础,目前国内核行业通常使用经验贝叶斯方法估计可靠性参数。然而极大似然函数对于稀少样本的估计存在求解困难的问题。马氏链蒙特卡洛方法提供了另一种可用于可靠性参数估计的分层贝叶斯方法,其从后验分布大量抽样进而推断总体特征。采用基于马氏链蒙特卡洛方法的Win BUGS软件计算了核电站需求失效型稀少样本的超参数。针对68个核电站辅助给水电动部分启动需求失效数据,计算得到启动失效的超参数α为0.192 5,β为46.77,以及各个核电站辅助给水电动部分启动需求失效的后验失效概率。最后对Win BUGS软件中马氏链的收敛性进行讨论。
关键词：马氏链蒙特卡洛方法,WinBUGS,核电站可靠性参数估计
参考文献
[1] Atwood C L,La Chance J L,Martz H F,et al.Handbook of parameter estimation for probabilistic risk assessment.NUREG/CR-6823,2003
[2] Eide S A,Wierman T E,Gentillon C D,et al.Industry-average performance for components and initiating events at U.S.commercial nuclear power plants.NUREG/CR-6928,2007
[3] 何劼,张彬彬,PSA通用数据库中超参数的估计方法研究.核动力工程,2012;33(5):83—86He Jie,Zhang Binbin.On hyperparameter estimation methods for generic database of PSA.Nuclear Power Engineering,2012;33(5):83 —86
[4] 茆诗松.贝叶斯统计.北京:中国统计出版社,2013Mao Shisong.Bayesian statistics.Beijing:China Statistics Press,2013
[5] Spiegelhalter D,Thomas A,Best N,et al.,Win BUGS user manual.Version 1.4.3,2007
[6] 刘乐平,张美英,李姣娇.基于Win BUGS软件的贝叶斯计量经济学.东华理工学院学报(社会科学版),2007;26(2):101—107Liu Lerping,Zhang Meiying,Li Jiaojiao.Bayesian econometrics based on Win BUGS.Journal of East China Institute of Technology,2007;26(2):101—107
参数估计软件 篇5
摘要：为提高电网状态估计准确性和查找状态估计问题速度，以临沂电网EMS系统为对象，研究了以动态参数维护及观测为基础的电力系统状态估计维护方法，实现了对拓扑、阻抗等多种参数变化情况的全面迅速掌握，并建立了实用的维护辅助决策系统。
关键词：状态估计;参数动态变化;辅助决策
中图分类号：TM711 文献标识码：A 文章编号：1007-0079（2014）32-0199-02
随着电网规模的不断扩大及电网自动化水平的不断提高，作为电网能量管理系统（energy management system，EMS）高级应用软件实用化的基础[1-3]，电力系统状态估计发挥着越来越重要的作用。而日趋复杂的电网结构和庞大电网规模为电网状态估计维护带来了巨大的困难[4-5]。本文从电网参数的动态变化入手，建立针对不同类型参数的分类管理监测机制，通过实用辅助系统的的建立，大大降低状态估计维护难度，提高状态估计维护效率。
一、状态估计影响因素分类
状态估计即根据SCADA提供的实时信息和网络拓扑的分析结果及其它相关数据，实时地给出电网内各母线电压（幅值和相位），各线路、变压器等支路的潮流，各母线的负荷和各发电机出力。影响状态估计的主要因素有拓扑连接关系、元件模型、元件电气参数、实时遥信和遥测，相互间关系如图1所示。
其数学如公式1所示：
（1）
式中，z为状态估计结果;为状态估计结果与参数的函数;m为原件模型;x为实时遥测与遥信;p为电气参数;s为拓扑链接关系。
在四种影响因素中，拓扑连接关系及元件模型决定了电网基本的结构，在建立维护中，出现错误的可能性较小，但一旦出错，对状态估计的影响极为巨大，很可能直接导致状态估计计算不收敛。元件电气参数数据庞大，人工填写工作量大，由于数字无规律且为多位小数，出错可能性相对较大，但出现错误后对状态估计影响较小，尤其是现有状态估计软件普遍具有参数辨识及自动修正功能，可有效减小参数出错情况下对状态估计数据的影响，与其他相比参数，主变档位数及运行档位对状态估计影响较大，在计算其影響时需区别对待。由于电压等级多，量测元件数目庞大，实时遥信/遥测是否对应在各种影响因素中是最容易出现问题的因素，尤其是开关、刀闸位置遥信位置错误对状态估计的准确性具有重大的影响，高电压等级开关、刀闸遥信位置的错误将会导致大量下级变电站状态估计的计算，某些枢纽变电站关键开关、刀闸位置的错误将会直接导致整个系统状态估计计算出现不收敛现象。
各种因素出错可能性及对状态估计影响总结如表1所示。
表1 状态估计影响因素总结
影响因素 出错可能性 对状态估计影响
拓扑连接关系 小 大
元件模型 小 大
元件电气参数 中 小（中）
实时遥信遥测 大 大
二、各影响因素主要问题分析
现代能量管理系统软件普遍采用图模库一体化功能，自动化维护人员在建立电网模型的过程中可有效避免图形元件同数据库的对应问题，从而大大减少了维护人员的工作量，但由于大电网的复杂性、时间的紧迫性及自动化人员水平差异等原因，往往会出现许多错误，从而影响状态估计的准确性。状态估计影响因素参数的主要错误具有不同的特点，每项影响因素对状态估计产生影响的主要问题的如表2所示。
表2 各影响因素主要问题
影响因素 主要问题
拓扑连接关系 各电气元件端口连接是否正确
元件模型 元件模型是否与实际一致
元件电气参数 阻抗与导纳等参数填写是否正确
实时遥信与遥测 是否正确关联元件，所取数据是否正确
三、评价体系建立
在复杂的电网结构下，电网某一元件出现参数错误往往影响有电气连接的多个元件的状态估计计算出现偏差，甚至导致全网状态估计计算错误，为排除冗余信息干扰，快速定位问题所在，需根据各影响因素的特点，建立科学的评价体系，对各个不合格的状态估计量测点进行排序，从而降低状态估计维护消缺工作的难度。
本文针对不同影响因素出现问题时的影响，对全网模型建立了影响值指标，以此表示元件自动化信息出现问题对状态估计指标降低的影响程度，其具体的表达式为：
式中：N为元件影响值;n1为拓扑连接关系影响值，其数值根据元件所连接的元件数及所连接的元件电压等级确定;n2为元件模型影响值，其数值根据元件类型（如母联开关、变压器等）及元件的电压等级确定;n3为元件电气参数影响值，其数值根据元件类型、是否填写及参数是否异常等因素确定;n4为实时遥信与遥测影响值，其数值跟军遥信/遥测是否对应、元件不平衡度及遥测信息是否突变等因素确定;ai为各因素相应权重，数值受由各种因素出错对状态估计影响决定，在确定权重值时，应突出元件模型及拓扑连接关系的先决性，其相应权重值较大。
当出现状态估计错误时，相应元件进行影响值计算，并根据影响值大小进行排序，自动化维护人员即可根据影响值排序快速定位问题所处元件。需要注意的是，由于状态估计影响的区域性，在判断问题潜在元件时，应适当扩大元件范围。
四、以动态参数变化为依据的状态估计维护方案
目前，各级电网对电网状态估计重视程度逐渐增大，指标要求不断提高，但相关维护辅助工具缺乏，仅采用断面提示信息的方式，该种方式具有很大的局限性，冗余干扰信息较多，尤其在状态估计出现计算不收敛时，相关提示不具有参考性。本文以评价体系及动态参数对比为基础，制定了一套针对动态参数对比的状态估计维护方案，如图2所示。以某一合格率较高的状态估计断面为基础，通过系统自动比对相关数据出现的变化，如元件的新增、拓扑关系的改变、电气参数的变动等，通过闭环的循环控制，使状态估计的合格率不断提高，可以有效的提高状态估计维护工作效率，减少自动化人员维护工作量。
在流程中，初始断面的选择影响着整个程序运行的效率，电网拓扑完善、元件模型恰当、电气参数填写正确、实时数据良好的断面将大大缩短系统进入稳定运行的时间，从而尽量少的引入人工维护处理，故建议使用最为接近的状态估计收敛的断面作为初始短命，并应将一些质量良好的断面作为备案资料进行保存，以便以后维护工作使用。
在判断环节，为排除负荷短时间波动对状态估计维护的影响，在判断合格率是否高于初始断面及是否下降时可设定一定的阈值，当上升/下降的数值高于阈值时判断为是，否则忽视其变动，正常运行。为保证其余电力系统高级应用的需要，当出现状态估计合格率下降后，状态估计不应停止计算，应在保存变动数据后，根据实时信息继续进行计算，从而为维护人员提供具有时标的电网信息变动记录，为电网数据维护提供类似事故追忆的参考资料，以发展的角度对比电网状态的变化，以便更为迅速查找问题的关键点，提高维护处理环节的效率。
基于该流程，开发了相应的状态估计辅助决策系统，该系统以图形可视的方式列出了状态估计计算出现不收敛等问题时参数的动态变化，如图3所示维护工作人员可根据维护需要，选择需查看的参数变化类型，如参数错误、主变实时档位错误等。在维护界面中，列出了前后状态的变化，并给出了出现问题的分析结果，有力地支持了状态估计维护工作的进行。
五、结束语
在建立状态估计元件评价体系的基础上，本文提出了基于动态参数变化的电力系统状态估计维护方法，并开发了相应的状态估计辅助决策系统。该辅助决策系统实现了潜在问题元件参数评价、参数状态变化等功能，解决了当前状态估计辅助维护工具缺乏、问题元件定位困难的现状，有助于提高网基础数据质量，为电力高级应用提供了坚实的基础。
参考文献：
[1]李碧君，薛禹胜，顾锦汶，等.电力系统状态估计问题的研究现状和展望[J].电力系統自动化，1998，22（11）：53-60.
[2]董树锋，何光宇，孙英云，等.以合格率最大为目标的电力系统状态估计新方法[J].电力系统自动化，2009，32（16）：40-43.
[3]顾全，陆杏全.电力系统实用状态估计中两个问题的处理[J].电力系统自动化，1998，22（1）：32-35.
[4]于尔铿.电力系统状态估计[M].北京：水利水电出版社，1985.
[5]何光宇，董树锋.基于测量不确定度的电力系统状态估计：（一）结果评价[J].电力系统自动化，2009，33（19）：21-24.
参数估计软件 篇6
在完成概率论部分及统计基本概念的教学后, 就进入参数估计。在本校使用的教材中主要是点估计和区间估计两种, 中间还有点估计的评判标准[1]。点估计一般课本中介绍两种方法:矩估计和极大似然估计。这两种方法需要进行积分、级数和求导, Excel不能解决这类问题, 所以这部分仍然需要运用高等数学的知识进行理解和解答。本文主要是利用Excel中的函数简化区间估计的计算过程。
正态总体有两个参数, 均值μ和方差σ2。
1正态总体均值的区间估计
来自正态总体的样本为x1, x2, …xn在数理统计的基本概念中关于样本均值、样本方差的计算还有分位数的查找都已进行了说明, 但是如果样本容量较大的话计算比较麻烦。在理解原理的基础上, 利用软件计算可以简化过程, 减少学生学习的抵触情绪。
大部分实际情况下σ2都是未知的, 主要介绍这部分。
例:设某种砖头的抗压强度X~N (μ, σ2) , 今随机抽取20块砖头, 测得数据 (单位:kg·cm-2) 如下:64, 69, 49, 92, 55, 97, 41, 84, 88, 99, 84, 66, 100, 98, 72, 74, 87, 84, 48, 81。求μ置信度为的95%置信区间, α=0.05。
分析:首先需要计算样本均值x軃和样本标准差s, 在Excel中输入20个数据, 利用统计函数AVERAGE, STDEV[2]计算出x軃=76.6, s=18.14 097。然后根据函数TINV (α, n-1) 得出本题的分位数t0.025 (19) =2.093 024, 再就是按照区间估计的公式计算相应的结果。利用计算公式76.6-2.093 024*18.140 97/SQRT (20) , 算出置信区间的下限为68.109 77, 将上面公式里的减号改成加号得到上限85.090 23, 所以μ置信度为95%的置信区间是 (68.11, 85.09) 。
如果是σ2已知的情况, 就不用计算样本标准差, 使用标准正态分布的分位数, 使用的函数为NORMSINV (α/2) 。当样本容量大于46时, T分布近似于标准正态分布, 也就是说这个时候公式里的T分布分位数由标准正态分布的分位数代替。
2正态总体方差σ2的区间估计
使用1中例题的数据, 计算σ2置信度为95%的置信区间, α=0.05。
分析:按照公式, 需要计算样本方差, 也就是样本标准差的平方, 这部分已经计算得s2=329.094 7。再就是需要χ2分布的两个分位数, 利用函数CHIINV (α/2或1-α/2, 自由度) , 得出χ20.025 (19) =32.852 33, χ20.975 (19) =8.906 517。再利用计算公式得出置信区间下限190.330 5, 上限702.047 8, 所以σ2置信度为95%的置信区间是 (190.330 5, 702.047 8) 。如果是要σ的置信区间, 只需要将上限下限开方就可以。
统计学部分的公式其实比较容易理解, 但是通常计算起来比较麻烦, 利用软件计算, 可以简化过程, 省掉大量的计算时间, 那么就可以在课堂上增加更多的知识点, 或者给出更为实际的例题给学生练习, 提高学生的学习兴趣。数理统计对大部分的理工科专业来说是一种工具, 如果学习过于枯燥, 会降低学生学习的积极性, 这样也不利于后面知识的掌握。课堂中使用Excel教学, 也可以提高学生的应用能力。对于数学专业的学生, 需要学习的内容会更多, 计算量也会更大, 使用软件同样可以提高学习的效率, 而不会使他们觉得统计学就是不停地简单计算, 一点实用性都没有。
摘要：本文在说明了区间估计的基本原理后, 使用EXCEL完成数理统计中相应的例题, 并说明这类教学方式的优点。
关键词：参数,置信度,置信区间
参考文献
[1]韩旭里, 谢永钦.概率论与数理统计[M].复旦大学出版社, 2010.
寿命分布的参数估计 篇7
关键词：寿命分布,生存函数,危险函数,指数分布,极大似然估计,区间估计
在保险精算及医学等的研究中, 需要对人的寿命作统计分析。工业生产中, 需要考察产品的可靠性, 可靠性一般用产品的寿命来衡量。寿命数据的统计模型可以用参数模型、非参数模型和半参数模型等。参数模型因其理论成熟浅易、计算简单、在实践中操作性强, 更易为人们接受和选用。所谓的参数模型, 是指寿命分布为某种已知类型的分布, 但有一些未知参数, 当这些参数确定后, 寿命分布就完全确定了, 因此, 对寿命分布的推断就可以转化为对参数的推断。常见的寿命分布模型, 有指数分布, 韦布尔分布, 伽玛分布, 对数正态分布等。指数分布是非常重要的一类概率分布, 将其作为寿命数据的统计模型是再恰当不过的, 它背景良好、理论易于理解, 相比其它几个分布, 有更好的应用价值。我们将讨论指数分布下寿命参数的估计。
1 寿命分布及有关的参数
2 指数分布e (λ) 的寿命参数
3 参数估计
设寿命X的分布为指数分布e (λ) , λ>0为未知参数。很容易得到λ的极大似然估计为 。令 , 则θ的极大似然估计为 。由于寿命数据的特殊性, 对其作观测时, 不大可能观测到样本的准确值, 例如, 要考察一种延寿药的效果, 不可能等到样本 (服药人群) 都正常死亡后才作统计分析。为此, 常见的寿命数据多半为截尾数据, 常用的截尾数据有两种类型:定数截尾数据和定时截尾数据。仍以考察一种延寿药的效果为例, 设服药人群的样本容量为n, 在做截尾试验前, 指定服药的样本中死亡个体的个数 , 当观测到第k个个体死亡时试验结束, 这样的试验叫定数截尾试验。如果在做截尾试验前, 指定一个时间长度c, 从0时刻起, 当观测时间到c时刻时试验结束, 这样的试验叫定时截尾试验。
3.1 定数截尾数据下对参数的估计
设样本容量为n, 截尾个数为 , 记k个观察到的寿命的顺序统计量为 , 似然函数为:
延拓Vasicek模型参数估计 篇8
近几年来, 经济和金融的领域中, 短期利率的模型被广泛地应用。这是因为短期利率模型能够很好的刻画瞬时短期利率的变动规律。正是这个原因, 目前有很多的短期利率模型的被提出。同时这类的模型很好地被应用到一些利率的衍生品的定价, 如:债券、利率互换、利率的远期合约等利率衍生品。
Vasicek (1977) 建立Gauss扩散的短期利率模型。但这个模型有一个不如人意, 瞬时的利率有可能是负值。Coxetal. (1985) (CIR) 建立根均值的短期利率模型。相对Vasicek, 在参数满足Feller条件下, CIR模型的瞬时利率是正值。但这两个模型都是考虑参数是常数。参数是常数可能不适合市场的数据。众所周知, 市场的利率是有期限结构的性质。而常数的参数, 很难刻画出这种性质。Ho和Lee (1986) 提出部分参数是时间的函数, 但是在Ho和Lee模型中, 波动率仍然是一个常数, 也就是在任何时刻波动率是常数。这种想象也有点违背市场的数据。因此, Hull和White (1990) 延拓了Vasicek模型和CIR模型,
在Hull和White模型中, 所有的参数都是时间的函数, 并在他们的文章, 利用这两个模型对于利率互换的产品定价分析。进一步, Rogers (1995) 在理论上分析Hull和White模型。
而在这篇文章中, 我们将根据市场的数据来估计Hull和White (1990) 延拓的Vasicek模型。考虑这个模型是因为这个模型有很好的分析性质, 特别是应用到定价其他的利率衍生品。
二、模型
首先, 假设短期利率r在风险中性测度下满足如下的随机微分方程:
其中, a (t) , b (t) , σ (t) 均为时间函数, W (t) 为标准的Brown运动, 设θ (t) ={a (t) , b (t) , σ (t) }如果是常数时候, (1) 是Vasicek模型。
基于 (1) , 在时间时刻, 到期日支付为1单位的零息债券价格有如下表达式,
其中, E (·∣·) 表示条件数学期望, 表示为在t时刻的观察值。如果是一个常数时, 基于Vasicek模型 (2) 有解析表达式, 这些结果可以参考Hull (2009) 以及姜礼尚等 (2008) 。考虑到参数是时间时候, Rogers (1995) 也证明 (2) 有显示的表达式。
根据Feynman-Kac公式, 在鞅的测度下, (1) 适合下面的微偏分方程,
众所周知 (3) 有一个仿射结构解,
其中, A (T, T) =1, B (T, T) =0。也就是在到期日零息债券价格为一个单位。即使, A (t, T) 和B (t, T) 有解析的表达式, 但他们是多重积分计算问题 (参考Rogers (1995) ) 。因此, 我们把这类问题归结为求方程的问题。
基于 (4) , A (t, T) 和B (t, T) 满足下面的常微分方程
其中, ø (t) =a (t) b (t) 其终端条件 (Model 1) 能够被发现在Hul和White (1990) 。
相应Model1的离散格式, 设Δti=ti+1, i=1, …, n, 其中t1=t, tn=T为了简化说明, 设A (t, ti) =Ai, A, B (t, ti) =Bi, i=1, …n, Model 1的离散格式为:
其中ai=a (ti) , σi=σ (ti) , Øi=a (ti) , (5) 终端条件, An=1, Bn=0
基于 (4) , 债券价格为
其中, Pn=P (ri, t, tn) 。因此, 对于不同的到期日, (6) 能算出其相应的价格。
我们的问题是在给定市场零息债券价格下, 估计Model 1中的参数, 首先, 假设在t时刻观察到零息债券价格。理论上, 当T足够大时候, 我们能估计未来任何时刻的参数值, 事实上, 市场上的零息价格期限是有限, 因此我们将讨论基于离散的市场数据来估计依赖时间函数的参数。
其次, 在我们的计算方法中, 需要计算理论的价格和市场价格的残值, 这里我们将使用常微分方程计算。即使债券价格有显示表达式, 但它是一个三重积分计算问题, 这就可能产生离散的误差。使用偏微分方程计算时有显著的误差 (看图1) 。
图1基于Vasicek模型通过常微方程和偏微分计算零息债券价格数值解与精确解的比较, 这里假设参数是常数, a=0.15, b=0.08, σ=0.1, r∈[0, 1], T=5基于Vasicek模型, 当r比较小时, 常微分方程的数值解比较好, 这归因于计算其偏微分方程时, 其在r大于零的值依赖于r小于零的值。另一原因在于偏微分数值结果依赖对空间变量r的离散误差。而常微分方程数值解仅依赖t, r是一个外在的变量。
三、参数估计计算方法
为了简化问题我们假设当前时刻t=0, 进一步我们假设有N个市场价格PM= (ri, 0, T) , i=l, …, N这就意味着到期日相同的债券有不同价格, 这个根源于初始利率的不同, 这个假设隐含着, 市场上债券价格通过同一个模型而得到的。
进一步我们假设, 市场价格满足如下方程
其中, p (ri, 0, T, θ (t) ) 被称为理论价格, εi表示理论价格和市场价格之间的误差, 直接求解 (7) 可能是不适定, 这个归因于理论价格可能是非线性的形式解而且市场价格报价存在误差。这就导致来自 (7) 的解可能是不适定的, 因此, 我们可以把问题转化为求下列极小值的问题
其中, Ω表示所在领域, 求 (8) 极小解, , 可通过非线性最小二乘法, 如果直接求解 (8) 时, 对于θ (0) 和θ (T) 至少有一个值是不能被计算, 这个依赖于求解理论价格时离散格式。因此我们将通过三次样条差值来计算端点的问题。
四、实证结果
我们使用美国国家债券每天交易市场数据, 这篇文章所有使用的数据都来源于http://www.ustreas.gov.由于其报价是通过收益率曲线来实现, 因此需要通过现值表达式转化为债券的价格。设, 当前时刻t, 在T时刻其相应在收益率曲线上收益率为, 那么其相应的零息债券价格为
如:t=0, y=0.015, T=5, 其债券价格为0.92774.
我们考虑2009年整年的到期为5年债券每天数据, 同时我们限制所有参数取值落在
1. 另一类方法是正则化方法, 这种方法不仅仅可以解决端点的问题, 而且处理参数光滑性的问题, 我们将在今后研究中将使用这类的方法估计参数。
2. 选择美国国家债券数据由于这类债券几乎是无风险债券, 而对于一般债券数据由于隐含违约风险因素, 可能不满足我们所考虑的模型。
3. 该收益率曲线刻画是实际收益率而不是名义收益率。
区间[0, 1]。对于初始利率将适应美国国库券到期为一周的每天收益率为初始利率。
表1:这个表格描述Model 1的参数估计数值结果, 最后一行是 (8) 的残值, 所有的数据都乘以100。
数据, 对于所有的数据变化大约或不超过0.1%, 这就表明参数有足够的光滑性。
为了测试算法的稳定性问题, 我们在收益率和初始利率上添加扰动项, 即:
其中δ1和δ2是待定常数, ε1和ε2是均值为0方差为1的正太分布的随机变量。
在表2中, 我们列出了RMES的值, 通过观察表2中数据, 对于初始利率相对于初始的收益率数据更加稳定, 根均值误差几乎没有什么变化, 虽然对于初始的收益率变动其根均值有微小的变化, 但这样的变化整体来说是可以接受的, 如我们预期的一样, 我们的计算方法是稳定。
五、结论
我们主要考虑基于零息债券价格估计Model 1的参数, 其中参数是依赖时间函数。我们所使用方法是众所周知的最小二乘法。端点的处理通过三次样条插值计算而得到。通过数值结果可知, 基于债券市场数据, 参数估计值是稳定的, 同时参数的估计值的光滑性比较好。
参考文献
[1]O.A.Vasicek, An Equilibrium Characterization of the Term Structure, Journal of Financial Economics, 5 (1977) :177-188.
[2]J.C.Cox, J.E.Ingersoll, and S.A.Ross, A Theory of the Term Structureof Interest Rate, Econometrica, 53 (1985) :385-407.
[3]T.S.Y.Ho and S.B.Lee, Term Structure Movements and Pricing InterestRate Contigent Claims, Journal of Finance, 41, (1986) :1011-1029.7
[4]J.Hull and A.White, Pricing Interest-Rate-Derivative Securities, The Review of Financial Studies, Vol.3 (4) (1990) , 573-592.
[5]L.G.G., Rogers, Which model for term-structure of interest rates shouldon use·Mathematical Finace, 65, 93-116, 1995
[6]J.C.Hull, Options, Futures and Other Derivatives (7th ed.) , PrenticeEducation, Inc., 2009.
平面参数的稳健估计方法探讨 篇9
三维激光扫描技术作为获取测量数据的一种手段,在三维城市、古文物保护、隧道检测等方面获得了广泛的应用。由于受到传感器的限制、遮挡、多回波反射和噪声的影响,会产生很多偏离物体表面的点,为点云数据的处理带来了困难。平面拟合、估计平面参数是点云分析的基础工作。最小二乘(LS)、主成分分析(PCA)和随机采样一致性(RANSAC)算法是平面拟合的常用算法。很多学者对此进行了研究,提出了很多改进算法。官云兰等提出一种平面拟合算法以特征值最小二乘法拟合平面
Rousseeuw等提出的S估计
1 算法介绍
1.1 S估计和Fast S估计
对于一个样本集X={xi∈Rp,i=1,…,n},S估计可定义为:
式(1)中:S为p×p正定矩阵,m∈Rp,ρ为目标函数,当样本服从正态分布时,常令,E代表期望,Φ0代表正态分布。
目标函数满足两个条件:①函数是偶函数,二阶连续可导;②ρ(0)=0,在区间[0,c0]严格单调递增,在区间[c0,+∝]是常数。(c0>0)。通常选用Tukey函数:
式(2)中:c为常数。
Fast S估计利用重抽样算法计算S估计,其主要思想如下
式(1)中的S改写为σ2Γ(|Γ|=1,σ=Σ1/(2p)),因此式(1)可以改写为:
式(3)中:G为p×p正定矩阵,|G|=1,s为正实数。
首先从观测值中抽取N个容量为p+1的子样本,计算每个子样本的均值向量和协方差矩阵获得N个初始估计
然后对这些初始估计执行K次I-step,其中第j次I-step如下:
Step1计算尺度和权重wi(j)
Step2计算加权均值向量和加权协方差矩阵,更新
接着根据式(4)保持不变,迭代计算直至收敛,
选择t个尺度最小的估计量,进行I-step直至收敛,t个估计量中尺度最小的为最终解最后协方差矩阵为
1.2 Det S估计
根据Det MCD
1.3 稳健距离与异常值检测
通过计算样本到样本集中心之间的距离可以来判别异常值。但是对于三维点云,这种方法并不有效,需要顾及协方差矩阵。马氏距离顾及了样本集的均值向量和协方差矩阵。对于m维多元样本集P={pi∈Rm,i=1,...,n}马氏距离定义为:
式(7)中cMD为均值向量,Σ为协方差矩阵。
Rousseeuw等指出,通过马氏距离可以很容易检测出一个异常值,但是当样本中存在多个异常值时,因为异常值之间的相互影响马氏距离不再有效
因此本文利用Det S估计获得的均值向量和协方差用于公式(7),得到稳健的距离:
Rousseeuw等指出稳健距离服从m(变量维度)自由度的χ2分布,稳健距离大于可以判定为异常值
1.4 本文方法
本文首先利用Det S估计获得数据的均值向量和协方差矩阵,再利用改进的马氏距离检测异常值,得到含有较少异常值的数据,最后利用选权迭代法进一步减弱异常值的影响,得到平面参数的最佳估值。
2 实验
为了验证本文方法的有效性和稳健性,分别通过仿真数据和真实场景数据对本文方法进行验证并与经典算法进行对比。
2.1 仿真数据实验
在仿真试验中采用两个指标来评价方法的有效性:
(1)拟合的平面单位法向量与真实的平面单位法向量n0的夹角
(2)平面的拟合参数和平面的真实参数p=(a,b,c,d)之间的相关系数
根据平面方程0.41519x+0.80525y-0.42332z+0.52627=0生成100个点,分别用整体最小二乘(TLS)、主成分分析(PCA)和本文方法拟合平面参数,结果如表1所示,三种方法都可以获得较精确的参数估值。
表1 不含异常值的拟合结果Table 1 Fitting results of data without outlier
为了探索不同异常值含量对拟合方法的影响,按照以下步骤进行试验:
(1)给定异常值含量百分率(0%~60%);
(2)随机生成1000个平面参数(平面参数设计值);
(3)根据平面参数和异常值含量百分率生成含异常值的点集(100个点);
(4)根据不同的算法计算以上1000个平面参数估值,统计每种算法的θ和r和通过稳健距离探测出的异常值个数。
含有25%的异常值的统计结果如表2、表3和表4所示。从表2和表3可以看出,TLS和PCA算法含有较大的偏差。本文方法几乎不受异常值影响,具有很高的稳健性。表4展示了根据稳健马氏距离所探测到的异常值个数的分布,从中可以看出平均97.2%的异常值被探测出。
表2 θ统计表Table 2 Statistical table ofθ
表3 r统计表Table 3 Statistical table of r
表4 探测出的异常值个数统计表Table 4 Statistical table of number of outliers detected
不同含量异常值的计算结果如图1、图2、图3所示。从图1(a)和图2(a)可以看出,当含有较少的异常值时,TLS和PCA算法失效,本文方法具有较好的稳健性。从图1(b)和图2(b)可以看出,当异常值的含量(<50%)很高时,本文方法依然有效。但是当异常值含量大于50%,本文方法不再有效。从图3可以看出,绝大多数异常值被探测到。
2.2 真实场景实验
分别选择不同的数据进行实验,由于篇幅限制选择2个有代表性数据结果进行展示。数据一选择的是一建筑物立面上的一块平面,该平面被周围的灌木遮挡,含有一定量的异常值。数据二选择的是一道路面,路面上有落叶,碎屑,坑洼等,产生一定含量的异常值。结果如表5、表6所示,本文方法可以获得较为精确的平面参数估值。
表5 数据1三种方法的拟合结果Table 5 Fitting results of data1 using three methods
图1 异常值百分率—θ平均值Fig.1 Line diagrams for outlier percentage vs averageθ
图2 异常值百分率—r的平均值Fig.2 Line diagrams for outlier percentage vs average r
图3 异常值百分率—探测到的异常值个数的平均值Fig.3 Line diagrams for outlier percentage vs average number of outliers detected
表6 数据2三种方法的拟合结果Table 6 Fitting results of data2 using three methods
3 结语
参数估计软件 篇10
L 1+ae-bx
中参数L，a，b的估计往往是基于一定的经济意义或生物统计的背景，但从纯粹的数学方程角度来看，当L，a，b的取值范围不加任何限制时，是否也能估计出L，a，b呢？文章证明了当（L，a，b）
∈R3
时，y=
L 1+ae-bx
中参数L，a，b不能用最小二乘法估计，并给出了其参数在一定的限制条件下，可以用最小二乘法估计.
[关键词] 逻辑曲线 参数 最小二乘估计
[中图分类号] G633.6 [文献标识码] A [文章编号] 1674 6058（2016）17 0049
一、引言
基于残差平方和最小的思想，最小二乘法已在社会生产和实践中得到了广泛的应用.在解决非线性回归模型的参数估计问题时，我们往往将其转化为线性回归模型，进而用最小二乘法估计出其参数，但有的模型并不能通过直观的观察和简单的计算就能转化为线性回归模型，也就谈不上用最小二乘法估计其参数了.这就引发了我们对非线性回归模型的参数是否能用最小二乘法估计的沉思，如常见的逻辑曲线y=
L 1+ae-bx
（（L，a，b∈R3）），其参数L，a，b能否用最小二乘法的思路估计出或至少能被确定属于一定的范围呢？这就是本文所要研究的重点.
二、猜想
任意取n个不同的样本（xi，yi），i=1，2，…，n.只要样本中不含用（0，y（0））；那么逻辑曲线y= L 1+ae-bx （（L，a，b）∈R3）的参数L，a，b就不能用最小二乘法估计出且不能估出（L，a，b）属于R3中某一真子空间，即无法确定（L，a，b）∈H，HR3.
三、证明
假设用最小二乘法可以估计出（L，a，b），尽管（L，a，b）并不一定唯一，即（L，a，b）也可属于R3中一真子空间.
第一步，观察y= L 1+ae-bx ，
∵（xi，yi）可取n个不同的样本，
∴可推断出L的估计，L≠0.
第二步，（Ⅰ）y= L 1+ae-bx （*）两边同时对x求导，得
dy dx =by（1-
y L
）
，其中，y（0）=
L 1+a
（**），
∴参数（L，a，b）的估计若满足（*）式，则也满足（**）.
（Ⅱ）考查微分方程y= L 1+ae-bx ，y（0）= L 1+a （**）.
令f（x，y）=by（
1- y L ）
，取R3平面上一区域D：|x| ≤h1|y|≤h2.
不难发现：f（x，y）在D上连续，且
|f（x1，y1）-f（x2，y2）|= 1 L |bL
（y1-y2）
-b（y21-y22）|≤
b|y1-y2|+2
bh2 L
|y1-y2|=（b+2 bh2 L ）|y1-y2|.
故f（x，y）在D上满足Lipschitz条件.
∴由常系数微分方程的解存在唯一性定理知，（**） 式满足y（0）= L 1+a 的解在|x|≤h上存在且唯一.
这里的h=min{h1，h2}，M=max|f（x，y）|，（x，y∈D）.
（Ⅲ）现在解（**）式：
1）观察f（x，y），得：y=0或y=L为（**）的解.
2）在f（x，y）中，令z= 1 y ，则
dz dx =
dz dy × dy dx =
- 1 y2 · b L y（L-y）
=-bz+
b L .
（1）
再令z=c（x）e-bx，则
dz dx =
dc（x） dx e-bx-
be-bxc（x）
= dc（x） dx e-bx
-bz. （2）
比较（1）（2）得： dc（x） dx = b L ebx.
利用变量分离法，求得：c（x）= 1 L ebx
+C，C为任一常数.
∴y= 1 z =
1 c（x）
ebx=
L ebx+CL ebx=
L 1+CLe-bx .
（3）
将y（0）=
L 1+a
代入（3），得：C= a L .
∴y= L 1+ae-bx .
综上，（**）式的解为y=0或L或 L 1+ae-bx .
但由（Ⅰ）知，样本（xi，yi）对（**）式中参数L，a，b的估计问题同样适用.故实际操作中，由于n个样本（xi，yi）是不同的，所以y=0，L要排除.
综合（Ⅱ）（Ⅲ）知，在整个实数域 R 上，（**）式的解都是存在参数（L，a，b）作最小二乘估计.
第三步，多元线性回归方程的一般形式是y=β0+β1x1+…+βpxp，不失一般性，总可以设β0=0.因为可以形式上引入自变量x0=1，所以在理论研究时，可以假设多元线性回归模型的形式为y=β1x1+…+βpxp.
由假设知，y= L 1+ae-bx 中参数a，b可作最小二乘估计，故y= L 1+ae-bx 必可划为线性回归形式：g（x，y）=p（L）u（x，y）+q（a）u（x，y）+h（b）w（x，y）.（4）
由于我们是先估计出p（L），q（a），h（b）后，利用L=p（L），
q（a）=q（a），h（b）=h（b），进而估计出，L，a，b的，所以（4）式与y= L 1+ae-bx 实质上对L，a，b的估计并没有任何区别.
故不妨设y= L 1+ae-bx 可划为y=Lx1+ax2+bx3（5）.
（5）式中，残差θi=
（yi-Lxi1-axi2-bxi3）=（L，a，b）xx′（Lab）-2y′x（Lab）+y′y.
其中，x=x11x22…xnn，y=y1y2…yn，
∴2xx′（Lab）-2x′y=0.
由假设知，L，a，b可用最小二乘估计出：若（L，a，b）唯一，则r（xx′）=3；若（L，a，b）不唯一，则（L，a，b）属于R3中一真子空间，此时r（xx′）=1或2.
综上，由假设得到r（xx′）>0.
第四步，由第二步知，对y= L 1+ae-bx 中参数L，a，b的估计等价于对 dy dx
=by（1- y L ），y（0）= L 1+a 中参数L，a，b的估计.
分析后者，我们发现通过样本（xi，yi）至多只能直接估出b，L，而a与样本（xi，yi）并没有直接关系，它只与y（0）有直接联系，a= L y（0） -1（∵由第一步知，L≠0，∴y（0）= L 1+a ≠0）.
但实际的样本数据中，由于（0，y（0））并未给出，故y（0）的取值范围为R＼{0}.
现任意取定L，b，一方面残差θ退化为只含a的二次函数，即r，s，t∈R1→θ=θ（a）=ra2+sa+t，
则a的估计值必须满足θ′（a）=0.即2ra+s=0. （6）
另一方面，a=
L y（0）
-1，由L≠0，y（0）的取值范围为R＼{0}，得到a的取值范围为R＼{-1}，故R＼{-1}中任一点都可作为a，显然也要满足（6）式.
∴r=s=0，
∴θ（a）=t为常数，故残差θ与a无关.
同理，任意取定a，b，也可得到残差θ与L无关.
故残差θ至多与b有关.
如果θ与b有关，则θ=θ（b），即y= L 1+ae-bx 可转化为线性回归形式，g（x，y）=h（b）w（x，y）；但前式含有3个参数L，a，b，后式只含有1个参数b，显然，这种转化不成立.
故残差θ=t为常数.
反馈到第三步，可知此时的xx′=03×3，即r（xx′）=0，与假设得到的r（xx′）>0矛盾.
∴原假设不成立，故猜想成立.
证毕.
四、总结与归纳
尽管L，a，b在非限制条件下，y= L 1+ae-bx 中的参数L，a，b不能用最小二乘法估计，但在实际经济模型或生物统计模型中，L，a，b往往具有特殊的经济含义和现实意义.我们可以用其他方法先估计出L或b等，然后再利用最小二乘法估计剩余的两个参数.比如，著名的Logistic人口发展模型y= L 1+ae-bx 中的参数L代表自然资源和环境条件所能容纳的最大人口数量，我们可以先根据生态学的知识预测出L，然后将y= L 1+ae-bx 转化为线性回归形式，即
-bx+lna=ln（ L y -1）
，再利用不同年份的人口统计数据和最小二乘法估计出a，b.
[ 参 考 文 献 ]
[1] 王高雄.周之铭等.常微分方程[M].北京：高等教育出版社，2007.
[2] 北大数学系几何与代数教研室前代数小组.高等代数[M].北京：高等教育出版社，2003.
[3] 陈家鼎，孙山泽等.数量统计学讲义[M].北京：高等教育出版社，2011.
[4] 吴梅村.数理统计学基本原理和方法[M].成都：西南财经大学出版社，2006.
【参数估计软件】相关文章：
多普勒参数估计06-10
两参数Weibull分布联合置信区间估计04-11
适合软件无线电的载波频偏估计技术研究09-11
ghost 参数详细说明软件教程04-10
空间估计05-10
系数估计06-02
函数估计06-04
估计预测06-26
交通状态估计05-29
密度核估计05-30