第1章 非线性最小二乘法与数值最优化
变量之间的关系更多地表现为非线性特征。线性模型作为基础模型是非线性的近似，即任何非线性模型都可以通过线性模型来近似表达。比如，模型01x y e u ββ=++通过泰勒级数展开表述为
0000100101**01|()x x x x x y e x x u
e x e x u x u
βββββββ=≈+-+ =-++ =++
模型201y x u ββ=++的线性近似表达式为
0010201010**01(2)|()22x x y x x x u
x x x u
x u
βββββββ=≈+-+ =-++ =++
例 1.1 利用Monte Carlo 模拟的方法观察线性模型对非线性模型的近似。
设DGP 为：y=10+0.2*exp(x)+u ，x 为[1，3]区间的均匀分布。利用线性模型与指数模型分别回归模型，并计算x 对y 的（平均）边际影响与（平均）弹性。（数据文件：nonlin ）
但线性模型对非线性模型的近似程度取决于高阶部分是否充分小。即使在样本内线性模型能够较好地拟合数据，也不能准确地体现变量的结构关系。非线性模型中，x 对y 的边际影响（或弹性）是变化的；而线性模型中，x 对y 的边际影响（或弹性）是常数。很多情况下，线性模型与非线性模型对边际影响或弹性的估计存在非常大的差异。另外，利用线性模型拟合非线性数据存在潜在的危险，即区间外预测会存在越来越大的误差。因此，正确设定模型的形式是进行准确推断和预测的重要环节。
对于一般的回归模型，如以下形式的模型，
(,)f =+y X βu 1.1
OLS 一般不能得到其解析解。比如，运用OLS 方法估计模型（1.1），令S(β)表示残差平方和，即 2
211()[(;)]n n
i
i i i i S u y f ====-∑∑βX β 1.2 最小化S(β)，即根据一阶条件可以得到
1
(;)()2[(;)]n i i i i f S y f =??=--=??∑X ββX β0ββ 以模型y x u γαβ=++为例，其一阶条件为
2011
0()1[]02i n
x i i S y e ββββ=?=---=?∑β
22011
1()1[]02i i n
x x i i S y e e βββββ=?=---=?∑β 220121
2()1[]02i i n
x x i i i S y e x e ββββββ=?=---=?∑β 上述方程组没有解析解，需要一般的最优化方法。很多数值最优化算法都可以完成这一类任务，这些方法的总体思路是一样的。即，从初始值出发，按照一定的方向搜寻更好的估计量，并反复迭代直至收敛。各种不同的最优化算法的差异主要体现在三个方面：搜寻的方向、估计量变化的幅度和迭代停止法则。本章主要介绍非线性最小二乘法、梯度最优化方法及其其在Stata 中的实现。
1.1 非线性最小二乘法
非线性最小二乘法的思路是，通过泰勒级数将均值函数展开为线性模型。即，只包括一阶展开式，高阶展开式都归入误差项。然后再进行OLS 回归，将得到的估计量作为新的展开点，再对线性部分进行估计。如此往复，直至收敛。
设模型中共存在（k +1）个参数β =（β0, β1 , …, βk ）。首先为参数选择一组初始值00,01,0,0(,,)k βββ=β 。其中下标零表示初始值。然后将f (X , β) 按泰勒级数在β0点展开。
0(0)(0)(,)(,)'()f f R =+-+X βX βg ββ 1.3 其中g 0表示一阶导数在 β0=( β0, 0 , β1, 0 , …, βk , 0)时的取值，R 为高阶部分。上式中只保留β的线性部分，将高阶部分归入误差项，可得
0(0)(0)(0)0(0)(0)1
(,)(,)'()'(,)'f f R f =+ =+-++ =+-+y X βu
X βg ββu g βX βg βu 1.4
其中，随机扰动项u 1包含u 和泰勒级数展开式中高阶成分。得到新的回归模型
0(0)(0)(0)1(,)''f -+=+y X βg βg βu 1.5
新的目标函数为*11??'Q =u
u 。模型的OLS 估计量为， ()()()()1
(1)(0)(0)(0)(0)(0)(0)1(0)(0)(0)(0)(0)'(,)''(,)f f --=-+ =+-βg g g y X βg ββg g g y X β 1.6
因此，NLS 的迭代估计式为：
()()1
(1)()()()()()'(,)j j j j j j f -+=+-ββg g g y X β 1.7 上述方法也被称作高斯牛顿方法，即利用迭代方法估计非线性最小二乘模型。在（1.7）式中，第二项恰好是模型
()()(,)j j f -=+y X βg γv 1.8 γ的OLS 估计量。把1.8式称作高斯-牛顿回归。β的迭代公式又可以写作
(1)()()j j j +=+ββγ 1.9 这种迭代估计方法必须设定初始值和停止法则。初始值的选择对于迅速找到最优解非常重要。如果目标函数不是严格的凹函数或凸函数，或者存在多个局部最优值，可以设定多个初始值，观察最优解。如果不同的初始值得到相同的最优解，则结论是比较稳健的。
停止法则用以设定满足一定的标准后终止迭代过程，否则迭代过程会无限继续下去。可用的停止法则包括：目标函数Q (β(j+1))- Q (β(j))没有明显的变化，或者g (j)的每个元素都非常小，或者β(j+1) - β(j)没有明显变化等。迭代法则也可以同时设定最高迭代次数n ，如果经过n 次迭代仍然没有能够达到收敛，则停止迭代。
例 1.2 利用NLS 方法估计非线性消费函数。（数据文件：usmacro ）
cs inc u γαβ=++
利用NLS 方法估计模型。
. nl (realcons = {alpha} + {beta}*realgdp^{gamma})
参数初始值可以通过两种方法进行设定。其一，直接在参数表达式中设定。比如，上述模型中设定alpha=0，beta=0.5，gamma=1。Stata 命令可以表述为：
. nl (realcons = {alpha=0} + {beta=0.5}*realgdp^{gamma=1})
另外一种方法是通过命令选项initial 进行设定。比如上述命令也可以表述为：
. nl (realcons = {alpha} + {beta}*realgdp^{gamma}), initial(alpha 0 beta 0.5 gamma 1)
注意，Stata 的nl 命令用于估计不存在缺失值的区间。如果某些变量存在缺失值，则必须通过variables 选项设定模型中的解释变量。比如，
. nl (realcons = {alpha} + {beta}*x^{gamma=1}), variables(realgdp)
在一般的非线性模型中，解释变量对被解释变量的结构影响（边际影响或弹性）都不是常数，甚至影响的方向与参数估计量的符号也不同。比如，在本模型中，收入对消费的边际影响为βγ?inc γ-1。Stata 提供了一个方便的计算边际影响的命令mfx ，用于计算不同变量在不同点上的边际影响。但如果在命令栏中使用nl 命令，必须设定variables 选项，才能使用mfx 命令。比如，计算平均边际消费倾向
. nl (realcons = {alpha} + {beta}*x^{gamma=1}), variables(realgdp)
. mfx, at(mean)
计算收入inc=3000美元时的边际消费倾向，
. mfx at(realgdp=3000)
例 1.3 利用NLS 方法估计如下生产函数模型。文件包含产出、资本、劳动力等数据。（数据文件：production ）
（1） 广义CD （GCD ）生产函数
012ln ln ln y y K L u θβββ+=+++
实际上，GCD 模型为线性模型，仍然可以采用nl 命令进行估计，与regress 命令得到完全相同的结果。
. nl (lnout=-{theta}*out + {beta0} + {beta1}*lnk+{beta2}*lnl)
（2） 不变替代弹性（CES ）生产函数
()0ln 1/ln (1)y K L u ρρβρδδ--=-+-+
. nl ( lnout={beta0}-1/{rho=1}*ln( {delta=0.5}*capital^(-{rho})
+(1-{delta})*labor^(-{rho}) ) )
资本劳动力的要素替代弹性为1/(1+rho)。
. nlcom (1/(1+_b[rho]))
（3） 广义CES 生产函数
()0ln 1/ln (1)y y K L u ρρθβρδδ--+=-+-+
. nl ( lnout={beta0}-1/{rho=1}*ln( {delta=0.5}*capital^(-{rho})
+(1-{delta})*labor^(-{rho}) ) )
1.2 最优化方法概述
非线性最小二乘法是利用泰勒级数展开的一种迭代估计方法，它是一般数值最优化的特例。对于一般的非线性的模型的其它估计方法，如极大似然法或广义矩方法等，都涉及到数值最优化的问题。接下来，我们介绍数值最优化方法的一般原理及其在Stata 中的实现。
设参数为β =（β0, β1 , …, βk ），目标函数为Q (β)。估计参数β使目标函数最优化（极大化或极小化）。
1.2.1 格点搜寻法
格点搜寻法即令参数在其取值范围内取不同的数值，从中选择使目标函数最优的数值作为参数估计量。比如，模型中只有一个参数β，已知β介于[0，1]。我们可以令β从0逐渐变化到1，每次递增0.01，即0，0.01，0.02，…, 1。对于每一个取值计算目标函数，进而得到使目标函数最优的β估计量。格点搜寻法只适用于参数个数比较少，而且取值范围具有一定的信息的模型。如果参数个数较多，而且取值范围较大的话，那么格点搜寻会耗费很长时间。比如，5个参数，每个参数取10个数值，那么需要计算105次。
1.2.2 迭代方法
另外一类最优化方法是迭代方法，其中一类最常见的方法即是梯度方法。以极大化问题为例，迭代方法的基本公式为：
(1)()()()j j j j +=+ββA g 1.10 其中，g (j)为一阶导数梯度矩阵（gradient ），元素为/i Q β??。
根据泰勒级数的线性展开式(0)(0)(0)()()'()Q Q =+-ββg ββ，第j 次估计量表示为β(j+1)，则目标函数Q (β(j+1))在β(j)的展开式为：
(1)()()(1)()()()'()j j j j j Q Q R ++=+-+ββg ββ 1.11 其中，R 表示被忽略的高阶项。将迭代公式代入可得
(1)()()()()()()'j j j j j Q Q R +-=+ββg A g 1.12 对于极大化问题，如果A 为正定矩阵且R 比较小的话，则Q (β(j+1))> Q (β(j))。对于极小化问题，如果A 为负定矩阵且R 比较小的话，则Q (β(j+1)) < Q (β(j))。
1.2.3 其它方法
梯度方法要求目标函数充分平滑，可以计算其梯度向量。如果目标函数不存在梯度向量，则需要利用其它最优化方法。比如，最小绝对值方法的目标函数是，
11(,),()|||(;)|n n
i i i i i f S u y f ===+==-∑∑y X βu βX β
可以利用线性规划方法。其它的更复杂的模型需要用到模拟退火法（annealling ）或者遗传算法（genetic ）。
1.3 梯度方法
1.3.1 牛顿（拉弗森）方法
初始值为(0)0,01,0,0(,,)k βββ=β 。利用泰勒级数展开得到目标函数的二阶近似表达式：
*(0)(0)(0)(0)(0)(0)1()()'()()'()2
Q Q =+-+--ββg ββββH ββ 其中，g (0)为一阶导数梯度矩阵（gradient ），元素为/i Q β??；H (0)为二阶导数海塞矩阵（Hessian ），元素为2/i j Q ββ???。如果函数为二次（或一次）线性函数，则泰勒级数展开式为精确展开。
根据一阶条件
*(0)(0)(0)()/()Q ??=+-=ββg H ββ0
可得β的第一次估计量
1(1)(0)(0)(0)-=-ββH g
如果Q *为凹函数，即H 为正定矩阵，则β(1)是最小值。如果Q *可以作为Q 的较好的近
似的话，β(1)也是?β
（目标函数Q 的估计量）的较好的近似。Newton 方法即是利用公式 1(1)()()()j j j j -+=-ββH g
进行迭代计算。
根据泰勒级数的线性展开式(0)(0)(0)()()'()Q Q =+-ββg ββ，目标函数Q (β(j+1))在β(j)的展开式为：
(1)()()(1)()()()'()j j j j j Q Q R ++=+-+ββg ββ 1.13 其中，R 表示被忽略的高阶项。将迭代公式代入可得
1(1)()()()()()()'j j j j j Q Q R -+-=-+ββg H g 1.14 对于极小化问题，如果H 为正定矩阵且R 比较小的话，则Q (β(j+1)) < Q (β(j))。对于极大化问题，如果H 为负定矩阵且R 比较小的话，则Q (β(j+1)) > Q (β(j))。
实践中，梯度方法对应两种修正方法。其一，根据（1.5）式，如果A 比较小，则目标函数的变化也比较小，每次迭代的进展速度会比较慢。反之，如果A 比较大，则有可能导致高阶项R 比较大，将其忽略会产生较大的误差。因此，对梯度方法的一个修正是加上一个调整因子α。这种方法称之为修正的牛顿方法。
1(1)()()()()j j j j j α-+=-ββH g 1.15 α使得下式最大化，()
?1()()()
()j j j Q Q αα-=-βH g 。 其二，如果出现H 为奇异矩阵的情况，则上式无法进行计算。这时，需要对奇异矩阵进行调整。比如，令H 加上某个矩阵H ，得到D =(H +C )-1。利用下式进行迭代估计。 1(1)()()()()j j j j j α-+=-ββD g 1.16 这种方法称之为拟牛顿方法。其中，D (j)作为H (j)的近似，总是正定矩阵。如果D (j)=H (j)，则得到了修正的Newton 方法。如果D (j)=H (j)，α(j)=1，则得到了普通的Newton 方法。
多情况下，比如Q *不是严格的凹函数，Newton 方法可能不能正确地找到β的最优值。 实际上，对比NR 公式可以发现，高斯牛顿方法是牛顿拉弗森方法的一个特例。 拟牛顿方法的基本步骤如下。
第一步，计算g (j)和D (j)，用于决定下一步的搜寻方向1()()j j -H g 。
第二步，计算α(j)，常用计算方法是最小化如下函数
()
?1()()()()j j j Q Q αα-=-βH g 第三步，根据停止法则，判断继续迭代还是停止迭代。
1.3.2 得分法
与NR 方法相比，得分法将其中的海塞矩阵H 替换为H 的期望矩阵的估计量，即 2?E 'MS Q ???=??????
H βββ 1.3.3 BHHH 法
Berndt, Hall, Hall, and Hausman (1974)方法将海塞矩阵H 替换为
?'BHHH Q Q
??=-??∑βH ββ
BHHH 估计只需要计算目标函数的一阶导数，而无需计算二阶导数，因此计算比较方便。
1.3.4 陡峭爬山法
陡峭爬山法直接令H=I 。
1.3.5 DFP 、BFGS 方法
Davidson, Fletcher, and Powell 是通过下面的迭代公式计算H 矩阵，
11'
'''',j j j j j j j j j j j j j j j j j j j j
++=++=+ =-δδA γγA A A g δγγA γδA g γg g
Boyden, Fletcher, Goldfarb, and Shannon 对其进行了进行了进一步修正。
1'
''''''(/')(/')
j j j j j j j j j j j j j j j j j j j j j j j j j j j j +=++-= -δδA γγA A A g γA γηηδγγA γηδδγA γγA γ
1.4 最优化方法在Stata 中的实现
Stata 中通过Mata 的最优化函数optimize 对目标函数进行最优化。基本步骤包括三步。设定目标函数，设置最优过程，进行最优化并观察估计结果。
1.4.1 设定目标函数
目标函数的一般格式为：
void evaluator (todo, p , X 1[,X 2, ... , X k ],v , g , H )
其中，todo 为0、1、2。0表示只设定目标函数（默认选项），1表示设定目标函数和梯度向量，2表示设定目标函数、梯度向量和海塞矩阵。P 为待估参数，v 为目标函数。g 、H 表示梯度向量和海塞矩阵。X 1 [,X 2, ... , X k ]表示模型中的变量，最多可以设置9个。
Stata 设定了三种类型的估值器：v0、v1、v2。如果采用v1，则需要用户设定梯度向量的表达式；如果采用v2，则需要用户设定海塞矩阵的表达式。
模型中的变量必须通过下面的命令进行设置。如果有两组变量X 、y ，则进行如下设置。 void evaluator (todo, p , X , y , v , g , H )
optimize_init_argument(S, 1, X )
optimize_init_argument(S, 2, y )
1.4.2设置最优化过程
最优化过程的控制主要包括最优化类型（最大化、最小化）、估值器类型、最优化方法、初始值、迭代停止法则等。
1．初始化最优设置
命令格式为：
S = optimize_init()
这一命令将最优化结果保存到S中，并将optimize_init_*()的设置恢复到默认状态。如果需要更改设置，通过后面的各种optimize_init_*()进行单独设定。
optimize_init_evaluator(S, & evaluator())
这一命令将最优化设定赋值到目标函数
2．设置最优化类型
命令格式为：
optimize_init_which(S, {"max" | "min"} )
默认选项为max。
3．设置估值器类型
命令格式为：
optimize_init_evaluatortype(S, evaluatortype)
evaluatortype包括："v0", "v1", "v2", "v1debug", "v2debug"。默认选项
为d0。
注：对于数学函数的最优化，evaluatortype包括"d0", "d1", "d2", "d1debug", "d2debug"。
4．最优化方法
命令格式为：
optimize_init_technique(S [,technique] )
technique包括：
"nr" 修正的Newton-Raphson
"dfp" Davidon-Fletcher-Powell
"bfgs" Broyden-Fletcher-Goldfarb-Shanno
"bhhh" Berndt-Hall-Hall-Hausman
"nm" Nelder-Mead
5．奇异海塞矩阵情况的处理方法
命令格式为：
optimize_init_singularHmethod(S [, singularHmethod] ) singularHmethod包括：
"m-marquardt" 修正的Marquardt算法
"hybrid" 垂直爬升法和牛顿方法的混合
6．初始值的设定
命令格式为：
optimize_init_params(S, real rowvector initialvalues)
7．停止法则的设定
设置最高迭代次数：
optimize_init_conv_maxiter(S, real scalar max)
设置收敛标准：
对于nr、dfp、bfgs、bhhh方法，收敛标准为：
（1）根据参数估计量的变化。即，mreldif(p, p_prior) < ptol
optimize_init_conv_ptol(S, real scalar ptol)
（2）根据目标函数的变化。即，reldif(v, v_prior) < vtol。
optimize_init_conv_vtol(S, real scalar vtol)
（3）根据一阶导数的变化。即，g*invsym(-H)*g' < nrtol。
optimize_init_conv_nrtol(S, real scalar vtol)
（4）对于极大化问题，要求-H为正半定矩阵。对于极小化问题，可以通过-f转为极大化问题。
当满足条件（1）、（3）、（4）时，或者满足（2）、（3）、（4）时，Stata将停止迭代。
默认的收敛标准为ptol=1e-6，vtol=1e-7，nrtol=1e-5。
8．设置跟踪程序
命令格式为：
optimize_init_tracelevel(S, string scalar tracelevel)
其中，tracelevel包括：
------------------------------------------------------
"none" 不进行记录
"value" 目标函数值（默认选项）
"tolerance" 收敛值
"params" 参数估计值
"step" 变化幅度
"gradient" 梯度向量
"hessian" 海塞矩阵
-----------------------------------------------------
1.4.3最优化并观察估计结果
进行最优化的命令格式：
p = optimize(S)
Stata通过optimize_result_*观察估计结果，包括参数估计量及其协方差矩阵、目标函数值、梯度向量、得分向量、海塞矩阵、迭代次数、是否收敛等。
参数估计量：real rowvector optimize_result_params(S)
目标函数最优值：real scalar optimize_result_value(S)
目标函数初始值：real scalar r optimize_result_value0(S)
梯度向量：real rowvector optimize_result_gradient(S)
得分向量：real matrix optimize_result_scores(S)
海塞矩阵：real matrix optimize_result_Hessian(S)
协方差矩阵：real matrix optimize_result_V(S)
协方差矩阵（仅对ML方法有意义）：
real matrix optimize_result_V_oim(S)
real matrix optimize_result_V_opg(S)
real matrix optimize_result_V_robust(S)
其中，oim 表示Cov=invsym(-H)（最大化）（对于最小化问题，Cov=invsym(H)），即直接利用观测到的信息矩阵（observed informatin matrix ）来计算。Opg 表示Cov=invsym(S'S)，其中S 表示得分向量，即利用梯度向量的外积公式（outer product of gradient ）。Robust 表示Cov=H*invsym(S'S)*H ，其中S 表示得分向量，即稳健协方差矩阵。
迭代次数：real scalar optimize_result_iterations(S )
是否收敛：real scalar optimize_result_converged(S )
如果达到收敛标准，则optimize_result_converged() = 1。
观察初始设定：optimize_query(S )
观察当前的收敛标准的设定值：
real scalar optimize_init_conv_ptol(S )
real scalar optimize_init_conv_vtol(S )
real scalar optimize_init_conv_nrtol(S )
1.5 最优化案例
例 1.4 最优化如下函数。
2min 3y x x =-+-
一阶导数为g=-2x +1，二阶导数为H=-2。根据NR 迭代公式1(1)()()()j j j j -+=-ββH g ，
估计量的迭代计算公式为()(1)()0.5(21)|j j j x x x x +=-?-+。设定初始值为x (0) = 0，则
(0)0(1)0.5(1)(0)(2)(1)0.5(21)|00.50.5
0.5(21)|0.500.5x x x x x x x x ===+?-+=+==+?-+=+=
目标函数不再变化，停止迭代。事实上，由于目标函数为严格凸函数，初始点无论选择在哪里，都可以经过一次迭代达到收敛。
Stata 的最优化命令如下。（程序文件：optim_func1.do ）。
mata:
mata clear
void myeval(todo, x, y, g, H)
{
y = -x^2 + x - 3
if (todo>=1) {
g = -2*x+1
if (todo==2) {
H = -2
}
} }
S = optimize_init()
optimize_init_evaluator(S, &myeval())
optimize_init_params(S, 10)
optimize_init_tracelevel(S,"step")
x = optimize(S)
end
例 1.5 最优化如下函数。
22121212exp(3)Q x x x x x x =++-++
梯度向量为
112221/(21)/(21)Q x x x Q Q x x x Q ??+-????==??????++????
g 海塞矩阵为
2221
11212122122221122121222(21)(21)(21)(21)(21)2(21)Q Q x x x x Q x x Q Q x x x x Q Q x x x x Q Q x x Q Q Q x x x x ????????????++-++-++??==????++++-+++??????????????
H 根据NR 迭代公式，1(1)()()()j j j j -+=-ββH g ，估计量的迭代计算公式为
1()1(1)1()12(1)2()|j j j x j j x x x x +-+????=-?????????????
H g 设定初始值为x 1 (0) = 0，x 2 (0) = 0。Stata 的命令函数如下（程序文件：optim_func2.do ）。 mata:
mata clear
void myeval(todo, x, y, g, H)
{
y = x[1]^2 + x[2]^2 +x[1]*x[2]-x[1]+x[2]+3
if (todo>=1) {
g[1]=(2*x[1]+x[2]-1)*y
g[2]=(2*x[2]+x[1]+1)*y
if (todo==2) {
H[1,1] = 2*y + (2*x[1]+x[2]-1)*g[1]
H[1,2] = y + (2*x[2]+x[1]+1)*g[1]
H[2,2] = 2*y + (2*x[2]+x[1]+1)*g[2]
_makesymmetric(H)
}
} }
S = optimize_init()
optimize_init_evaluator(S, &myeval())
optimize_init_which(S, "min")
optimize_init_params(S, (0,0))
optimize_init_tracelevel(S,"value")
x = optimize(S)
end
例1.6利用Stata的最优化方法估计例1.1中的消费函数模型。（数据文件：usmacro；程序文件：usmacro.do）
/* usmacro.dta */
mata:
mata clear
Data=st_data(., ("realcons", "realgdp"))
y=Data[.,1]
e=J(rows(Data),1,1)
X=e,Data[.,2::cols(Data)]
void eval(todo, p, X, y, v, g, H)
{
beta0=p[1]
beta1=p[2]
beta2=p[3]
mu=beta0*X[.,1] :+ beta1*X[.,2]:^beta2
v = (y:-mu):^2
}
S = optimize_init()
optimize_init_evaluator(S, &eval())
optimize_init_which(S, "min")
optimize_init_evaluatortype(S, "v0")
param=(100,0.7,1)
optimize_init_params(S, param)
optimize_init_argument(S, 1, X)
optimize_init_argument(S, 2, y)
p=optimize(S)
cov=optimize_result_V_oim(S)
st_matrix("Coef",p)
st_matrix("Cov",cov)
mata clear
end
matrix list Coef
matrix list Cov
例 1.7 利用Stata的最优化方法估计例 1.2中的各种生产函数模型。（数据文件：production）。
以CES函数为例。Stata的程序文件如下。（程序文件：ces.do）
mata:
Data=st_data(., ("lnout", "capital", "labor"))
param=(0,1,0.5)
e=J(rows(Data),1,1)
y=Data[.,1]
X=e,Data[.,2::cols(Data)]
void eval(todo, p, X, y, v, g, H)
{
beta0=p[1]
rho=p[2]
delta=p[3]
mu=beta0:-1/rho*ln(delta*X[.,2]:^(-rho)+(1-delta)*X[.,3]:^(-rho))
v = (y:-mu):^2
}
S = optimize_init()
optimize_init_evaluator(S, &eval())
optimize_init_which(S, "min")
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, param)
optimize_init_argument(S, 1, X)
optimize_init_argument(S, 2, y)
p=optimize(S)
cov=optimize_result_V_oim(S)
st_matrix("Coef",p)
st_matrix("Cov",cov)
mata clear
end
matrix list Coef
matrix list Cov
练习题
1．设DGP为：y=10+0.2*x^2+u，x为[1，3]区间的均匀分布。利用线性模型与二次多项式模型分别回归模型，并计算x对y的（平均）边际影响与（平均）弹性。（数据文件：nonlin2）
2．对于模型y=β0+β1*x+u，目标函数定义为。利用optimize函数估计模型。
非线性最小二乘法 编辑词条分享 ?新知社新浪微博腾讯微博人人网QQ空间网易微博开心001天涯飞信空间MSN移动说客 非线性最小二乘法 非线性最小二乘法是以误差的平方和最小为准则来估计非线性静态模型参数的一种参数估 计方法。 编辑摘要 目录 1 简介 2 推导 3 配图 4 相关连接 非线性最小二乘法 - 简介 以误差的平方和最小为准则来估计非线性静态模型参数的一种参数估计方法。设非线性系统的模型为y＝f(x，θ) 式中y是系统的输出，x是输入，θ是参数（它们可以是向量）。这里的非线性是指对参数θ的非线性模型，不包括输入输出变量随时间的变化关系。在估 计参数时模型的形式f是已知的，经过N次实验取得数据(x1,y1),(x2,y1)， ,(xn，yn)。估计参数的准则（或称目标函数）选为模型的误差平方和非线性最小二乘法就是求使Q达到极小的参数估计值孌。 推导 非线性最小二乘法 - 推导 以误差的平方和最小为准则来估计非线性静态模型参数的一种参数估计方法。设非线 性系统的模型为 y＝f(x，θ) 式中y是系统的输出，x是输入，θ是参数（它们可以是向量）。这里的非线性是指对参数θ的非线性模型，不包括输入输出变量随时间的变化关系。在估计参数时模型的形式f是已知的，经过N次实验取得数据(x1,y1),(x2,y1)， ,(x n，y n)。估计参数的准则（或称目标函数）选为模型的误差平方和
非线性最小二乘法就是求使Q达到极小的参数估计值孌。 由于f的非线性，所以不能象线性最小二乘法那样用求多元函数极值的办法来得到参 数估计值，而需要采用复杂的优化算法来求解。常用的算法有两类，一类是搜索算法，另 一类是迭代算法。 搜索算法的思路是：按一定的规则选择若干组参数值，分别计算它们的目标函数值并 比较大小；选出使目标函数值最小的参数值，同时舍弃其他的参数值；然后按规则补充新 的参数值，再与原来留下的参数值进行比较，选出使目标函数达到最小的参数值。如此继 续进行，直到选不出更好的参数值为止。以不同的规则选择参数值，即可构成不同的搜索 算法。常用的方法有单纯形搜索法、复合形搜索法、随机搜索法等。 迭代算法是从参数的某一初始猜测值θ(0)出发，然后产生一系列的参数点θ(1)、θ(2) ，如果这个参数序列收敛到使目标函数极小的参数点孌，那么对充分大的N就可用θ(N)作为孌。迭代算法的一般步骤是： ① 给出初始猜测值θ(0)，并置迭代步数i＝1。 ② 确定一个向量v(i)作为第i步的迭代方向。 ③ 用寻优的方法决定一个标量步长ρ(i)，使得 Q(θ(i))=Q(θ(i))，其中θ(i)＝θi-1+ρ(i)v(i)。 ④ 检查停机规则是否满足，如果不满足,则将i加1再从②开始重复；如果满足，则 取θ(i)为孌。 典型的迭代算法有牛顿－拉夫森法、高斯迭代算法、麦夸特算法、变尺度法等。 非线性最小二乘法除可直接用于估计静态非线性模型的参数外，在时间序列建模、连 续动态模型的参数估计中，也往往遇到求解非线性最小二乘问题。 非线性最小二乘法 - 配图 非线性最小二乘法
题目：用外点法求解非线性约束最优化问题 学院信息管理学院 学生姓名余楠学号 81320442 专业数量经济学届别 2013 指导教师易伟明职称教授 二O一三年十二月
用外点法求解非线性约束最优化问题 摘要 约束最优化问题是一类重要的数学规划问题。本文主要研究了用外点罚函数法对约束非线性规划问题进行求解。对于一个约束非线性规划用罚函数求解的基本思路是通过目标函数加上惩罚项，将原约束非线性规划问题转化为求解一系列无约束的极值问题。本文最后利用c语言编程得到满足允许误差内的最优解。 本文主要对一个约束非线性规划问题的实例，首先利用上述迭代的方法，计算出各迭代点的无约束极值问题的最优解。然后应用c语言编程，得到精确地最优解，需迭代四次次才使得ε≤0.001，得到的最优解为* X=(333 .0)T。 .0, 666 关键词：外点罚函数法非线性规划约束最优化迭代最优解
一、背景描述 线性规划的目标函数和约束条件都是决策变量的线性函数，但如果目标函数或约束条件中含有决策变量的非线性函数，就称为非线性规划。非线性规划与线性规划一样，也是运筹学的一个极为重要的分支，它在经济、管理、计划、统计以及军事、系统控制等方面得到越来越广泛的应用。 非线性规划模型的建立与线性规划模型的建立类似，但是非线性规划问题的求解却是至今为止的一个研究难题。虽然开发了很多求解非线性规划的算法，但是目前还没有适用于求解所有非线性规划问题的一般算法，每个方法都有自己特定的适用范围。 罚函数法是应用最广泛的一种求解非线性规划问题的数值解法，它的基本思路是通过目标函数加上惩罚项，将原约束非线性规划问题转化为求解一系列无约束的极值问题。这种惩罚体现在求解过程中，对于企图违反约束的那些迭代点，给予很大的目标函数值，迫使这一系列无约束问题的极小值点无限的向可行集（域）逼近，或者一直保持在可行集（域）内移动，直到收敛于原来约束问题的极小值点。 外点法的经济学解释如下：若把目标函数看成“价格”，约束条件看成某种“规定”，采购人员在规定的范围内采购时价格最便宜，但若违反了规定，就要按规定加收罚款。采购人员付出的总代价应是价格和罚款的综合。采购人员的目标是使总代价最小，当罚款规定的很苛刻时，违反规定支付的罚款很高。这就迫使采购人员在规定的范围内采购。数学上表现为罚因子足够大时，无约束极值问题的最优解应满足约束条件，而成为约束问题的最优解。 二、基础知识 2.1 约束非线性优化问题的最优条件 该问题是一个约束非线性优化问题，利用外点罚函数法求解该问题，约束非线性优化问题的最优解所要满足的必要条件和充分条件是我们设计算法的依据，为此有以下几个定理。
非线性曲线拟合最小二乘法 一、问题提出 设数据（i i y x ,）,(i=0,1,2,3,4).由表3-1给出，表中第四行为i i y y =ln ,可以看出数学模型为bx ae y =,用最小二乘法确定a 及b 。 i 0 1 2 3 4 i x 1.00 1.25 1.50 1.75 2.00 i y 5.10 5.79 6.53 7.45 8.46 i y 1.629 1.756 1.876 2.008 2.135 二、理论基础 根据最小二乘拟合的定义:在函数的最佳平方逼近中],[)(b a C x f ∈，如果f(x)只在一组离散点集{i x ,i=0,1,…,m}，上给定，这就是科学实验中经常见到的实验数据{（i i y x ,）, i=0,1,…,m}的曲线拟合，这里)(i i x f y =，i=0,1,…,m,要求一个函数)(*x S y =与所给数据{（i i y x ,）,i=0,1,…,m}拟合，若记误差 i i i y x S -=)(*δ,i=0,1,…,m,T m ),,(10δδδδ， =,设)(,),(),(10x x x n ??? 是] ,[b a C 上线性无关函数族，在)}(,),(),({10x x x span n ???? =中找一函数)(*x S ，使误差平方和 ∑∑∑===∈ -=-==m i m i m i i i x S i i i y x S y x S 0 2 )(2 * 2 22 ])([])([min ? δδ , 这里 )()()()(1100x a x a x a x S n n ???+++= (n 线性和非线性最优化理论、方法、软件及应用 最优化在航空航天、生命科学、水利科学、地球科学、工程技术等自然科学领域和经济金融等社会科学领域有着广泛和重要的应用, 它的研究和发展一直得到广泛的关注. 最优化的研究包含理论、方法和应用.最优化理论主要研究问题解的最优性条件、灵敏度分析、解的存在性和一般复杂性等.而最优化方法研究包括构造新算法、证明解的收敛性、算法的比较和复杂性等.最优化的应用研究则包括算法的实现、算法的程序、软件包及商业化、在实际问题的应用. 这里简介一下线性和非线性最优化理论、方法及应用研究的发展状况. 1. 线性最优化 线性最优化, 又称线性规划, 是运筹学中应用最广泛的一个分支.这是因为自然科学和社会科学中许多问题都可以近似地化成线性规划问题. 线性规划理论和算法的研究及发展共经历了三个高潮, 每个高潮都引起了社会的极大关注. 线性规划研究的第一高潮是著名的单纯形法的研究. 这一方法是Dantzig在1947年提出的,它以成熟的算法理论和完善的算法及软件统治线性规划达三十多年. 随着60年代发展起来的计算复杂性理论的研究, 单纯形法在七十年代末受到了挑战. 1979年前苏联数学家Khachiyan提出了第一个理论上优于单纯形法的所谓多项式时间算法--椭球法, 曾成为轰动一时的新闻, 并掀起了研究线性规划的第二个高潮. 但遗憾的是广泛的数值试验表明, 椭球算法的计算比单纯形方法差. 1984年Karmarkar提出了求解线性规划的另一个多项式时间算法. 这个算法从理论和数值上都优于椭球法,因而引起学术界的极大关注, 并由此掀起了研究线性规划的第三个高潮. 从那以后, 许多学者致力于改进和完善这一算法,得到了许多改进算法.这些算法运用不同的思想方法均获得通过可行区域内部的迭代点列,因此统称为解线性规划问题的内点算法. 目前内点算法正以不可抗拒的趋势将超越和替代单纯形法. 线性规划的软件, 特别是由单纯形法所形成的软件比较成熟和完善.这些软件不仅可以解一般线性规划问题, 而且可以解整数线性规划问题、进行灵敏度分析, 同时可以解具有稀疏结构的大规模问题.CPLEX是Bi xby基于单纯形法研制的解线性和整数规划的软件, CPLEX的网址是https://www.doczj.com/doc/924129047.html,/. 此外,这个软件也可以用来解凸二次规划问题, 且特别适合解大规模问题. PROC LP是SAS软件公司研制的SAS商业软件中OR模块的一个程序. 这个程序是根据两阶段单纯形法研制的,可以用来解线性和整数规划问题并可进行灵敏度分析, 是一个比较完善的程序.用户可以根据需要选择不同的参数来满足不同的要求。关于内点法的软件也在研制之中.BPMP D是Cs.Mzos基于原始对偶内点法研制的解线性和整数规划的软件，其FTP地址是ftp://ftp.sztaki.hu/pub /oplab/SOFTWARE/BPMPD/ ，可以自由下载.此外,在互联网上能访问到的解线性和整数规划问题的软件还有:EQPS（线性，整数和非线性规划）,FMP（线性和混合整数规划），HS/LPLO（线性规划），KORBX（线性规划），LAMPS（线性和整数规划），LPBLP（线性规划），MILP（混合整数规划），MINTO（混合整数规划），MPSIII（线性和混合整数规划），OML（线性和混合整数规划），OSL（线性，二次和混合整数规划），PROCLP（线性和整数规划），WB（线性和混合整数规划），WHIZARD（线性和混合整数规划），XPRESSMP（线性和混合整数规划）等. matlab最小二乘法的非线性参数拟合 首先说一下匿名函数：在创建匿名函数时，Matlab记录了关于函数的信息，当使用句柄调用该函数的时候，Matlab不再进行搜索，而是立即执行该函数，极大提高了效率。所以首选匿名函数。具体拟合时可以使用的方法如下： 1 曲线拟合工具箱提供了很多拟合函数，使用简单 非线性拟合nlinfit函数 clear all; x1=[0.4292 0.4269 0.381 0.4015 0.4117 0.3017]'; x2=[0.00014 0.00059 0.0126 0.0061 0.00425 0.0443]'; x=[x1 x2]; y=[0.517 0.509 0.44 0.466 0.479 0.309]'; f=@(p,x) 2.350176*p(1)*(1-1/p(2))*(1-(1-x(:,1).^(1/p(2))).^p(2)).^2.*(x(:,1).^ (-1/p(2))-1).^(-p(2)).*x(:,1).^(-1/p(2)-0.5).*x(:,2); p0=[8 0.5]'; opt=optimset('TolFun',1e-3,'TolX',1e-3);% [p R]=nlinfit(x,y,f,p0,opt) 2 最小二乘法在曲线拟合中比较普遍。拟合的模型主要有 1.直线型 2.多项式型 3.分数函数型 4.指数函数型 5.对数线性型 6.高斯函数型 一般对于LS问题，通常利用反斜杠运算“\”、fminsearch或优化工具箱提供的极小化函数求解。在Matlab中，曲线拟合工具箱也提供了曲线拟合的图形界面操作。在命令提示符后键入：cftool，即可根据数据，选择适当的拟合模型。 “\”命令 1.假设要拟合的多项式是：y=a+b*x+c*x^ 2.首先建立设计矩阵X： X=[ones(size(x)) x x^2]; 执行： para=X\y para中包含了三个参数：para(1)=a;para(2)=b;para(3)=c; 这种方法对于系数是线性的模型也适应。 2.假设要拟合：y=a+b*exp(x)+cx*exp(x^2) 设计矩阵X为 X=[ones(size(x)) exp(x) x.*exp(x.^2)]; para=X\y 3.多重回归(乘积回归) 设要拟合：y=a+b*x+c*t，其中x和t是预测变量，y是响应变量。设计矩阵为X=[ones(size(x)) x t] %注意x,t大小相等！ para=X\y 习题二包括题目：P36页5（1）（4） 5（4） 习题三 包括题目：P61页1(1)(2); 3; 5; 6; 14；15(1) 1(1)(2)的解如下 3题的解如下 5,6题 14题解如下 14. 设22121212()(6)(233)f x x x x x x x =+++---, 求点在(4,6)T -处的牛顿方向。 解：已知 (1) (4,6)T x =-，由题意得 121212212121212(6)2(233)(3)()2(6)2(233)(3)x x x x x x x f x x x x x x x x +++-----?? ?= ?+++-----?? ∴ (1)1344()56g f x -?? =?= ??? 21212122211212122(3)22(3)(3)2(233)()22(3)(3)2(233)22(3)x x x x x x x f x x x x x x x x +--+--------? ??= ? +--------+--?? ∴ (1)2(1)1656()()564G x f x --?? =?= ?-?? (1)1 1/8007/400()7/4001/200G x --?? = ?--?? ∴ (1)(1)11141/100()574/100d G x g -?? =-= ?-?? 15（1）解如下 15. 用DFP 方法求下列问题的极小点 （1）22 121212min 353x x x x x x ++++ 解：取 (0) (1,1)T x =，0H I =时，DFP 法的第一步与最速下降法相同 2112352()156x x f x x x ++???= ?++??， (0)(1,1)T x =，(0) 10()12f x ???= ??? (1)0.07800.2936x -??= ?-??， (1) 1.3760() 1.1516f x ???= ?-?? 以下作第二次迭代 (1)(0) 1 1.07801.2936x x δ-??=-= ?-??， (1)(0) 18.6240()()13.1516f x f x γ-??=?-?= ?-?? 0110 111011101 T T T T H H H H H γγδδδγγγ=+- 非线性最小二乘平差 6-1问题的提出 经典平差是基于线性模型的平差方法。然而在现实世界中，严格的线性模型并不多见。测量上大量的数学模型也是非线性模型。传统的线性模型平差中的很多理论在非线性模型平差中就不一定适用；线性模型平差中的很多结论在非线性模型平差中就不一定成立；线性模型平差中的很多优良统计性质在非线性模型平差中就不一定存在。例如，在线性模型平差中，当随机误差服从正态分布时，未知参数X 的最小二乘估计具有一致无偏性和方差最小性。但在非线性模型平差中，即使随机误差严格服从正态分布，未知参数X的非线性最小二乘估计也是有偏的。其方差一般都不能达到最小值。 对于测量中大量的非线性模型，在经典平差中总是进行线性近似（经典的测量平差中称之为线性化），即将其展开为台劳级数，并取至一次项，略去二次以上各项。如此线性近似，必然会引起模型误差。过去由于测量精度不高，线性近似所引起的模型误差往往小于观测误差，故可忽略不计。随着科学技术的不断发展，现在的观测精度已大大提高，致使因线性近似所产生的模型误差与观测误差相当，有些甚至还会大于观测误差。例如，GPS载波相位观测值的精度很高，往往小于因线性近似所产生的模型误差。因此，用近似的理论、模型、方法去处理具有很高精度的观测结果，从而导致精度的损失，这显然是不合理的。现代科学技术要求估计结果的精度尽可能高。这样，传统线性近似的方法就不一定能满足当今科学技术的要求。另外，有些非线性模型对参数的近似值十分敏感，若近似值精度较差，则线性化会产生较大的模型误差。由于线性近似后，没有顾及因线性近似所引起的模型误差，而用线性模型的精度评定理论去评定估计结果的精度，从而得到一些虚假的优良统计性质，人为地拔高了估计结果的精度。 鉴于上述各种原因，对非线性模型平差进行深入的研究是很有必要的。非线性模型的平差和精度估计以及相应的误差理论研究也是当前国内外测绘界研究的前沿课题之一。 电子教材 > 第六章非线性模型平差 > 6-2 非线性模型平差原理 Advances in Education教育进展, 2019, 9(4), 450-453 Published Online July 2019 in Hans. https://www.doczj.com/doc/924129047.html,/journal/ae https://https://www.doczj.com/doc/924129047.html,/10.12677/ae.2019.94076 A Good Textbook for Optimization —Review of Theory and Algorithms on Nonlinear Programming Naiyang Deng School of Science, China Agricultural University, Beijing Received: Jun. 30th, 2019; accepted: Jul. 12th, 2019; published: Jul. 22nd, 2019 Abstract This paper provides a review for the textbook written by Professors Yiju Wang and Naihua Xiu, named Theory and Method of Nonlinear Optimization published by Science Publishing House in 2012 by pointing out some highlights and features of the book, and giving some suggestions for improvement. Keywords Book Review, Highlights and Features, Suggestions 最优化课程的一本好教材 ——《非线性最优化理论与方法》书评 邓乃扬 中国农业大学理学院，北京 收稿日期：2019年6月30日；录用日期：2019年7月12日；发布日期：2019年7月22日 摘要 本文对王宜举教授和修乃华教授2012年在科学出版社出版的《非线性最优化理论与方法》进行了评述，指出了该书的一些亮点和特色，同时也给出了进一步提升的建议。 关键词 书评，亮点与特色，建议 非线性最小二乘lsqnonlin 数学规划模型的matlab求解 数学规划模型的matlab求解 var OsObject = ""; if(https://www.doczj.com/doc/924129047.html,erAgent.indexOf("MSIE")>0) { document.write(""); } if (isFirefox=https://www.doczj.com/doc/924129047.html,erAgent.indexOf("Firefox")>0){ document.write(" "); } if(isSafari=https://www.doczj.com/doc/924129047.html,erAgent.indexOf("Safari")>0) { //return "Safari"; } if(isCamin o=https://www.doczj.com/doc/924129047.html,erAgent.indexOf("Camino")>0){ //return "Camino"; } if(isMozilla=navigato https://www.doczj.com/doc/924129047.html,erAgent.indexOf("Gecko/")>0){ //return "Gecko"; } 今天胡老师给我们讲了数学规划模型，数学规划模型是优化模型的一种，包括线性规划模型(目标函数和约束条件都是线性函数的优化问题); 非线性规划模型(目标函数或者约束条件是非线性的函数); 整数规划(决策变量是整数值得规划问题); 多目标规划(具有多个目标函数的规划问题) ；目标规划(具有不同优先级的目标和偏差的规划问题) 动态规划(求解多阶段决策问题的最优化方法) 。数学规划模型相对比较好理解，关键是要能熟练地求出模型的解。 以下是解线性规划模型的方法： 1.线性规划问题 线性规划问题的标准形式为： min f ' *x sub.to：A*x 五种最优化方法 1.最优化方法概述 最优化问题的分类 1）无约束和有约束条件； 2）确定性和随机性最优问题（变量是否确定）； 3）线性优化与非线性优化（目标函数和约束条件是否线性）； 4）静态规划和动态规划（解是否随时间变化）。 最优化问题的一般形式（有约束条件）： 式中f（X）称为目标函数（或求它的极小，或求它的极大），si（X）称为不等式约束，hj（X）称为等式约束。化过程就是优选X，使目标函数达到最优值。 2.牛顿法 简介 1）解决的是无约束非线性规划问题； 2）是求解函数极值的一种方法： 3）是一种函数逼近法。 原理和步骤 3.最速下降法（梯度法） 最速下降法简介 1）解决的是无约束非线性规划问题； 2）是求解函数极值的一种方法； 3）沿函数在该点处目标函数下降最快的方向作为搜索方向； 最速下降法算法原理和步骤 4?模式搜索法（步长加速法） 简介 1）解决的是无约束非线性规划问题； 2）不需要求目标函数的导数，所以在解决不可导的函数或者求导异常麻烦的函数的优化问题时非常有效。 3）模式搜索法每一次迭代都是交替进行轴向移动和模式移动。轴向移动的目的 是探测有利的下降方向，而模式移动的目的则是沿着有利方向加速移动。 模式搜索法步骤 5.评价函数法 简介 评价函数法是求解多目标优化问题中的一种主要方法。在许多实际问题中，衡量一个方案的好坏标准往往不止一个，多目标最优化的数学描述如下： min (f_1(x),f_2(x),…,f_k(x)) .g(x)<=o 传统的多目标优化方法本质是将多目标优化中的各分目标函数, 经处理或数学变换，转变成一个单目标函数，然后采用单目标优化技术求解。常用的方法有“线性加权和法”、“极大极小法”、“理想点法”。选取其中一种线性加权 求合法介绍。 线性加权求合法 6.遗传算法 智能优化方法是通过计算机学习和存贮大量的输入-输出模式映射关系，进而达到优化的一种方法，主要有人工神经网络法，遗传算法和模拟退火法等。遗传算法基本概念 1.个体与种群 个体就是模拟生物个体而对问题中的对象 (一般就是问题的解)的一种称呼。种群就是模拟生物种群而由若干个体组成的群体，它一般是整个搜索空间的一个很小的子集。 2.适应度与适应度函数 适应度就是借鉴生物个体对环境的适应程度，而对问题中的个体对象所设计的表征其优劣的一种测度。 适应度函数就是问题中的全体个体与其适应度之间的一个对应关系。该函数就是遗传算法中指导搜索的评价函数。 遗传算法基本流程 的就是对一定数量个体组成的生物种群进行选择、交叉、变异等遗传操作，最终求得最优解或近似最优解。 遗传算法步骤 Levenberg-Marquardt Method(麦夸尔特法) Levenberg-Marquardt is a popular alternative to the Gauss-Newton method of finding the minimum of a function that is a sum of squares of nonlinear functions, Let the Jacobian of be denoted , then the Levenberg-Marquardt method searches in the direction given by the solution to the equations where are nonnegative scalars and is the identity matrix. The method has the nice property that, for some scalar related to , the vector is the solution of the constrained subproblem of minimizing subject to (Gill et al. 1981, p. 136). The method is used by the command FindMinimum[f, x, x0] when given the Method -> Levenberg Marquardt option. SEE A LSO:Minimum, Optimization REFERENCES: Bates, D. M. and Watts, D. G. N onlinear Regr ession and Its Applications. New York: Wiley, 1988. Gill, P. R.; Murray, W.; and Wright, M. H. "The Levenberg-Marquardt Method." §4.7.3 in Practical Optim ization. London: Academic Press, pp. 136-137, 1981. Levenberg, K. "A Method for the Solution of Certain Problems in Least Squares." Quart. Appl. Math.2, 164-168, 1944. Marquardt, D. "An Algor ithm for Least-Squares Estimation of Nonlinear Parameters." SIAM J. Appl. Math.11, 431-441, 1963. 统计学梯度下降法（SGDs）易于实现，然而它有两个主要的缺陷。第一个缺陷是它需要手动调谐大量的参数，比如学习速率和收敛准则。第二个缺陷是它本质上是序列方法，不利于并行计算或分布式计算。（然而，在计算资源如RAM受限的情况下，序列方法倒是一个不错的选择。） 这里介绍一些非线性优化算法：牛顿算法，伪牛顿算法和共轭梯度法。其中，伪牛顿算法包括DFP、BFGS和L-BFGS算法。 考虑如下的无约束最小化问题： min x f(x)（1） 其中x=(x1,…,x N)T∈?N. 为简便起见，这里假设f是凸函数，且二阶连续可导。记（1）的解为x?. 牛顿算法（Newton‘s Method） 基本思想：在现有的极小点估计值的附近对f(x)做二阶泰勒展开，进而找到下 其中g(k)=?f(x)| x(k)是梯度矩阵，H(k)=?2f(x)| x(k) 是海森矩阵。 牛顿算法是一种具有二次收敛性的算法。对于非二次函数，若函数的二次性态较强，或迭代点已进入极小点的领域，则其收敛速度也是很快的，这是牛顿算法的主要优点。但牛顿算法由于迭代公式中没有步长因子，而是定步长迭代，所以对于非二次函数，有时会出现f(x(k+1))>f(x(k))的情况，这表明牛顿算法不能保证函数值稳定地下降。由此，人们提出了阻尼牛顿算法，在原始牛顿算法的第4步中，采用一维搜索（line search）算法给d(k)加一个步长因子λ(k)，其中： λ(k)=arg minλ∈?f(x(k)+λd(k))（2）一维搜索算法将另作介绍。 拟牛顿算法（Quasi-Newton Methods） 基本思想：不直接计算二阶偏导数，而是构造出近似海森矩阵（或海森矩阵的逆）的正定对称阵，在拟牛顿条件下优化目标函数。 下文中，用B表示对H的近似，用D表示对H?1的近似，并令s(k)=x(k+1)?x(k)，y(k)=g(k+1)?g(k). 第四章最小二乘法与组合测量 §1概述 最小二乘法是用于数据处理和误差估计中的一个很得力的数学工具。对于从事精密科学实验的人们来说，应用最小乘法来解决一些实际问题，仍是目前必不可少的手段。例如，取重复测量数据的算术平均值作为测量的结果，就是依据了使残差的平方和为最小的原则，又如，在本章将要用最小二乘法来解决一类组合测量的问题。另外，常遇到用实验方法来拟合经验公式，这是后面一章回归分析方法的内容，它也是以最小二乘法原理为基础。 最小二乘法的发展已经经历了200多年的历史，它最先起源于天文和大地测量的需要，其后在许多科学领域里获得了广泛应用，特别是近代矩阵理论与电子计算机相结合，使最小二乘法不断地发展而久盛不衰。 本章只介绍经典的最小二乘法及其在组合测量中的一些简单的应用，一些深入的内容可参阅专门的书籍和文献。 §2最小二乘法原理 最小二乘法的产生是为了解决从一组测量值中寻求最可信赖值的问题。对某量x 测量一组数据n x x x ,,,21 ，假设数据中不存在系统误差和粗大误差，相互独立，服从正态分布，它们的标准偏差依次为：n ,,21记最可信赖值为x ，相应的残差x x v i i 。测值落入),(dx x x i i 的概率。 根据概率乘法定理，测量n x x x ,,,21 同时出现的概率为 显然，最可信赖值应使出现的概率P 为最大，即使上式中页指数中的因子达最小，即 权因子：2 2o i i w 即权因子i w ∝21i ，则 再用微分法，得最可信赖值x 1 1 n i i i n i i w x x w 即加权算术平均值 这里为了与概率符号区别，以i 表示权因子。 特别是等权测量条件下，有： 以上最可信赖值是在残差平方和或加权残差平方和为最小的意义下求得的，称之为最小二乘法原理。它是以最小二乘方而得名。 为从一组测量数据中求得最佳结果，还可使用其它原理。 例如 （1）最小绝对残差和法：Min v i （2）最小最大残差法：Min v i max （3）最小广义权差法：Min v v i i m in m ax 以上方法随着电子计算机的应用才逐渐引起注意，但最小二乘法便于解析，至今仍用得最广泛。 §3.线性参数最小二乘法 先举一个实际遇到的测量问题，为精密测定三个电容值:321,,x x x 采用的测量方案是，分别等权、独立测得323121,,,x x x x x x ,列出待解的数学模型。 1x =0.3 2x =-0.4 1x +3x =0.5 第1章 非线性最小二乘法与数值最优化 变量之间的关系更多地表现为非线性特征。线性模型作为基础模型是非线性的近似，即任何非线性模型都可以通过线性模型来近似表达。比如，模型01x y e u ββ=++通过泰勒级数展开表述为 0000100101**01|()x x x x x y e x x u e x e x u x u βββββββ=≈+-+ =-++ =++ 模型201y x u ββ=++的线性近似表达式为 0010201010**01(2)|()22x x y x x x u x x x u x u βββββββ=≈+-+ =-++ =++ 例 1.1 利用Monte Carlo 模拟的方法观察线性模型对非线性模型的近似。 设DGP 为：y=10+0.2*exp(x)+u ，x 为[1，3]区间的均匀分布。利用线性模型与指数模型分别回归模型，并计算x 对y 的（平均）边际影响与（平均）弹性。（数据文件：nonlin ） 但线性模型对非线性模型的近似程度取决于高阶部分是否充分小。即使在样本内线性模型能够较好地拟合数据，也不能准确地体现变量的结构关系。非线性模型中，x 对y 的边际影响（或弹性）是变化的；而线性模型中，x 对y 的边际影响（或弹性）是常数。很多情况下，线性模型与非线性模型对边际影响或弹性的估计存在非常大的差异。另外，利用线性模型拟合非线性数据存在潜在的危险，即区间外预测会存在越来越大的误差。因此，正确设定模型的形式是进行准确推断和预测的重要环节。 对于一般的回归模型，如以下形式的模型， (,)f =+y X βu 1.1 OLS 一般不能得到其解析解。比如，运用OLS 方法估计模型（1.1），令S(β)表示残差平方和，即 2 211()[(;)]n n i i i i i S u y f ====-∑∑βX β 1.2 最小化S(β)，即根据一阶条件可以得到 1 (;)()2[(;)]n i i i i f S y f =??=--=??∑X ββX β0ββ 以模型y x u γαβ=++为例，其一阶条件为 2011 0()1[]02i n x i i S y e ββββ=?=---=?∑β 精心整理 五种最优化方法 1.最优化方法概述 1.1最优化问题的分类 1）无约束和有约束条件； 2）确定性和随机性最优问题（变量是否确定）； 3 4 1.2 2. 2.1 1 2 3 2.2 3. 3.1 1 2 3 3.2 4.模式搜索法(步长加速法) 4.1简介 1）解决的是无约束非线性规划问题； 2）不需要求目标函数的导数，所以在解决不可导的函数或者求导异常麻烦的函数的优化问题时非常有效。 3）模式搜索法每一次迭代都是交替进行轴向移动和模式移动。轴向移动的目的是探测有利的下降 方向，而模式移动的目的则是沿着有利方向加速移动。 4.2模式搜索法步骤 5.评价函数法 5.1简介 评价函数法是求解多目标优化问题中的一种主要方法。在许多实际问题中，衡量一个方案的好坏标准往往不止一个，多目标最优化的数学描述如下： min(f_1(x),f_2(x),...,f_k(x)) s.t.g(x)<=0 传统的多目标优化方法本质是将多目标优化中的各分目标函数，经处理或数学变换，转变成一个单目标函数，然后采用单目标优化技术求解。常用的方法有“线性加权和法”、“极大极小法”、“理想点法”。选取其中一种线性加权求合法介绍。 5.2线性加权求合法 6.遗传算法 智能优化方法是通过计算机学习和存贮大量的输入-输出模式映射关系，进而达到优化的一种方法，主要有人工神经网络法，遗传算法和模拟退火法等。 6.1遗传算法基本概念 1.个体与种群 个体就是模拟生物个体而对问题中的对象（一般就是问题的解）的一种称呼。 种群就是模拟生物种群而由若干个体组成的群体,它一般是整个搜索空间的一个很小的子集。 2.适应度与适应度函数 适应度就是借鉴生物个体对环境的适应程度,而对问题中的个体对象所设计的表征其优劣的一种测度。 适应度函数就是问题中的全体个体与其适应度之间的一个对应关系。该函数就是遗传算法中指导搜索的评价函数。 6.2遗传算法基本流程 遗传算法的中心思想就是对一定数量个体组成的生物种群进行选择、交叉、变异等遗传操作,最终求得最优解或近似最优解。 遗传算法步骤 步1在搜索空间U上定义一个适应度函数f(x)，给定种群规模N，交叉率Pc和变异率Pm，代数T； 非线性约束最优化方法 ——乘子法 1.1 研究背景 最优化理论与方法是一门应用性相当广泛的学科，它的最优性主要表现在讨论决策问题的最佳选择性，讨论计算方法的理论性质，构造寻求最佳解的计算方法，以及实际计算能力。伴随着计算数学理论的发展、优化计算方法的进步以及计算机性能的迅速提高，规模越来越大的优化问题得到解决。因为最优化问题广泛见于经济计划、工程设计、生产管理、交通运输、国防等重要领域，它已受到政府部门、科研机构和产业部门的高度重视。然而，随着人们对模型精度和最优性的要求所得到的优化命题往往具有方程数多、变量维数高、非线性性强等特点，使得相关变量的存储、计算及命题的求解都相当困难，从而导致大规模非线性优化很难实现。因此，寻求高效、可靠的大规模非线性优化算法成为近年来研究的热点。 本文讨论的问题属于非线性约束规划的范畴，讨论了其中的非线性等式约束最优化问题方面的一些问题。 1.2非线性约束规划问题的研究方法 非线性约束规划问题的一般形式为 （NPL ） {}{} m in (),, s.t. ()0,1,2,...,, ()0,1,2,...,n i i f x x R c x i E l c x i I l l l m ∈=∈=≤∈=+++ 其中，(),()i f x c x 是连续可微的. 在求解线性约束优化问题时，可以利用约束问题本身的性质， 但是对于非线性约束规划问题，由于约束的非线性使得求解这类问题比较复杂、困难。因此，我们将约束问题转化为一系列无约束优化问题，通过求解一系列无约束优化问题，来得到约束优化问题的最优解。我们用到的几类主要的方法有:罚函数法、乘子法以及变尺度法。 传统上我们所提出的非线性约束最优化方法一般都遵循下列三个基本思路之一 1 借助反复的线性逼近把线性方法扩展到非线性优化问题中来 2 采用罚函数把约束非线性问题变换到一系列无约束问题 3 采用可变容差法以便同时容纳可行的和不可行的X 矢量 其中源于思路2 的乘子罚函数法具有适合于等式及不等式约束不要求初始点为严格内点，甚至不要求其为可行点对自由度的大小无任何要求等特点。 1.3乘子法 罚函数法的主要缺点在于需要惩罚因子趋于无穷大，才能得到约束问题的极小点，这会使罚函数的Hesse矩阵变得病态，给无约束问题的数值求解带来很大问题，为克服这一缺点，Hestenes和Powell 于1964年各自独立地提出乘子法。所谓乘子法是：由问题的Lagrange 函数出发，考虑它的精确惩罚，从而将约束优化问题化为单个函数的无约束优化问题，它同精确罚函数法一样，具有较好的收敛速度和数值稳定性，且避免了寻求精确罚函数法中关于罚参数阈值的困难，它们一直是求解约束优化问题的主要而有效的算法。 考虑如下非线性等式约束优化问题： %后面的为注释,红色部分代码需要根据实际情况更改 %最小二乘法线性拟合y=ax+b x=[0.5,1,1.5,2,2.5,3,3.5,4,4.5,5];%自变量 y=[191,321,442,565,686,819,930,1032,1153,1252];%因变量 xmean=mean(x);ymean=mean(y); sumx2=(x-xmean)*(x-xmean)'; sumxy=(y-ymean)*(x-xmean)'; a=sumxy/sumx2;%解出直线斜率a(即传感器灵敏度) b=ymean-a*xmean;%解出直线截距b z=((a*(x(1,10))+b-(y(1,10)))/(y(1,10)));%“10”是自变量的个数，z为非线性误差(即线性度) a b z %作图,先把原始数据点用蓝色"十"字描出来 figure plot(x,y,'+'); hold on % 用红色绘制拟合出的直线 px=linspace(0,6,50);%(linspace语法(从横坐标负轴起点0画到横坐标正轴终点6，50等分精度)) py=a*px+b; plot(px,py,'r'); 运行结果： a =236.9818 b =87.4000 另一种简单一点的方法： %最小二乘法线性拟合y=ax+b x=[0.5,1,1.5,2,2.5,3,3.5,4,4.5,5];%自变量 y=[191,321,442,565,686,819,930,1032,1153,1252];%因变量p=polyfit(x,y,1); p 运行结果： p = 236.9818 87.4000 五种最优化方法 1. 最优化方法概述 1.1最优化问题的分类 1）无约束和有约束条件； 2）确定性和随机性最优问题（变量是否确定）； 3）线性优化与非线性优化（目标函数和约束条件是否线性）； 4）静态规划和动态规划（解是否随时间变化）。 1.2最优化问题的一般形式（有约束条件）： 式中f(X)称为目标函数(或求它的极小，或求它的极大)，si(X)称为不等式约束，hj(X)称为等式约束。化过程就是优选X，使目标函数达到最优值。 2.牛顿法 2.1简介 1）解决的是无约束非线性规划问题； 2）是求解函数极值的一种方法； 3）是一种函数逼近法。 2.2 原理和步骤 3. 最速下降法（梯度法） 3.1最速下降法简介 1）解决的是无约束非线性规划问题； 2）是求解函数极值的一种方法； 3）沿函数在该点处目标函数下降最快的方向作为搜索方向； 3.2 最速下降法算法原理和步骤 4. 模式搜索法(步长加速法) 4.1 简介 1）解决的是无约束非线性规划问题； 2）不需要求目标函数的导数，所以在解决不可导的函数或者求导异常麻烦的函数的优化问题时非常有效。 3）模式搜索法每一次迭代都是交替进行轴向移动和模式移动。轴向移动的目的是探测有利的下降方向，而模式移动的目的则是沿着有利方向加速移动。 4.2模式搜索法步骤 5.评价函数法 5.1 简介 评价函数法是求解多目标优化问题中的一种主要方法。在许多实际问题中，衡量一个方案的好坏标准往往不止一个，多目标最优化的数学描述如下：min (f_1(x),f_2(x),...,f_k(x)) s.t. g(x)<=0 传统的多目标优化方法本质是将多目标优化中的各分目标函数，经处理或数学变换，转变成一个单目标函数，然后采用单目标优化技术求解。常用的方法有线性和非线性最优化理论、方法、软件及应用
matlab最小二乘法的非线性参数拟合
最优化计算方法课后习题答案----高等教育出版社。施光燕
非线性最小二乘平差
最优化课程的一本好教材 ——《非线性最优化理论与方法》书评
非线性最小二乘lsqnonlin
五种最优化方法
非线性最小二乘法Levenberg-Marquardt method
非线性优化算法-牛顿法_DFP_BFGS_L-BFGS_共轭梯度算法
参数的最小二乘法估计
Chapter2 非线性最小二乘法与数值最优化
五种最优化方法
关于非线性约束最优化方法-乘子法
用Matlab进行最小二乘法线性拟合求传感器非线性误差灵敏度
五种最优化方法
相关主题
文本预览