非线性最小二乘-Toy模板网

这篇具有很好参考价值的文章主要介绍了非线性最小二乘。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

非线性最小二乘

1 非线性最小二乘估计

在经典最小二乘法估计中，假定被解释变量的条件期望是关于参数的线性函数，例如
$E (y ∣ x) = a + b x$
其中 $a, b$ 为待估参数， $E (y ∣ x)$ 是关于参数 $a, b$ 的线性函数。但 $E (y ∣ x)$ 是关于参数的非线性函数，则利用ols求出的正规方程组没有解析解。只能通过相关数值计算。考虑一个简单的非线性模型
$Y_{i}=\beta X_{1 i}+\beta^{2} X_{2 i}+\varepsilon_{i}$
其中扰动项 $\varepsilon_i$ 满足 $\mathrm{E}\left(\varepsilon_{i}\right)=0,\operatorname{var}\left(\varepsilon_{i}\right)=\sigma^{2}$ ,且为独立同分布。其残差平方和为
$\begin{aligned} S(\beta) &=\sum_{i=1}^{n} \varepsilon_{i}^{2}=\sum_{i=1}^{n}\left[Y_{i}-f\left(X_{i}, \beta\right)\right]^{2} \\ &=\sum_{i=1}^{n}\left[Y_{i}-\beta X_{1 i}-\beta^{2} X_{2 i}\right]^{2} \end{aligned}$
为了使回归线尽可能接近观测值，要求残差平方和最小。根据微积分的知识
$\begin{aligned} \frac{\mathrm{d} S}{\mathrm{~d} \beta} &=2 \sum_{i=1}^{n}\left[Y_{i}-f\left(X_{i}, \beta\right)\right]\left(-\frac{\mathrm{d} f\left(X_{i}, \beta\right)}{\mathrm{d} \beta}\right) \\ &=2 \sum_{i=1}^{n}\left[Y_{i}-\beta X_{1 i}-\beta^{2} X_{2 i}\right]\left[-X_{1 i}-2 \beta X_{2 i}\right]=0 \end{aligned}$
整理得：
$\beta^{3} \sum_{i=1}^{n} X_{2 i}^{2}+3 \beta^{2} \sum_{i=1}^{n} X_{1 i} X_{2 i}+\beta\left(\sum_{i=1}^{n} X_{1 i}^{2}-2 \sum_{i=1}^{n} X_{2 i} Y_{i}\right)-\sum_{i=1}^{n} X_{1 i} Y_{i}=0$
这是关于参数 $\beta$ 的三次函数。尽管三次函数存在解析解（利用卡丹或盛金公式），其结果极为复杂。若上述三次方程存在实根 $\beta_i(i=1,2,3)$ （最多三个），则将 $\beta_i$ 代入残差平方和，取 $S(\beta)$ 最小所对应的 $\beta_i$ 。上述例子中，被解释变量条件期望是关于参数的二次函数，如果将这种函数形式改为指数、对数或三角函数形式，则一般不存在解析解。

因此，数值分析自然成为解决上述问题的有力武器。考虑一般化的非线性回归问题，设总体回归模型满足
$\beta)+\varepsilon$
对应的残差平方和为
$S(\beta)=\sum_{i=1}^{n}\left[Y_{i}-f\left(X_{i}, \beta\right)\right]^{2}$
要使其最小化，需要满足一阶条件
$\frac{\mathrm{d} S}{\mathrm{~d} \beta}=-2\left[\sum_{i=1}^{n}\left[Y_{i}-f\left(X_{i}, \beta\right)\right]\left(-\frac{\mathrm{d} f\left(X_{i}, \boldsymbol{\beta}\right)}{\mathrm{d} \beta}\right)\right]=0$
显然，上述问题不存在解析解，因此考虑对 $f(X_i, \beta)$ 进行一阶泰勒展开。设参数向量 $\beta$ 的初始值为 $\beta_1$ ,则可以在 $\beta_1$ 附近找到函数 $f(X_i, \beta)$ 使得
$f\left(X_{i}, \beta\right) \approx f\left(X_{i}, \beta_{1}\right)+\frac{\mathrm{d} f\left(X_{i}, \beta\right)}{\mathrm{d} \beta} \mid_{\beta = \beta_{1}}\left(\beta-\beta_{1}\right)$
记 $\left.\frac{\mathrm{d} f\left(X_{i}, \beta\right)}{\mathrm{d} \beta}\right|_{\beta_{1}} \approx \frac{f\left(X_{i}, \beta\right)-f(X, \beta)}{\beta-\beta_{1}}$ ，简记 $\widetilde{X}_{i}\left(\beta_{1}\right)=\left.\frac{\mathrm{d} f\left(X_{i}, \beta\right)}{\mathrm{d} \beta}\right|_{\beta_{1}}$ ，则
$\begin{aligned} S(\beta) &=\sum_{i=1}^{n}\left[Y_{i}-f\left(X_{i}, \beta_{1}\right)-\widetilde{X}_{i}\left(\beta_{1}\right)\left(\beta-\beta_{1}\right)\right]^{2} \\ &=\sum_{i=1}^{n}\left[\widetilde{Y}_{i}\left(\beta_{1}\right)-X_{i}\left(\beta_{1}\right) \beta\right]^{2} \end{aligned}$
其中
$\widetilde{Y}_{i}\left(\beta_{1}\right)=Y_{i}-f\left(X_{i}, \beta_{1}\right)+\widetilde{X}_{i}\left(\beta_{1}\right) \beta_{1}$
给定初始值向量 $\beta_i$ ，则 $\widetilde{Y}_{i}\left(\beta_{1}\right)$ 与 $\widetilde{X}_{i}\left(\beta_{1}\right)$ 可计算，从而求出最小残差平方和。 $S(\beta)$ 对应的回归方程为
$\widetilde{Y}_{i}\left(\beta_{1}\right)=\widetilde{X}_{i}(\beta) \beta+\varepsilon_{i}$
最小二乘估计量为
$\beta_{2}=\left[\widetilde{X}\left(\beta_{1}\right)^{\prime} \widetilde{X}\left(\beta_{1}\right)\right]^{-1} \widetilde{X}\left(\beta_{1}\right)^{\prime} \widetilde{Y}\left(\beta_{1}\right)$
其中
$\widetilde{X}\left(\beta_{1}\right)=\left[\begin{array}{c} \widetilde{X}_{1}\left(\beta_{1}\right) \\ \vdots \\ \widetilde{X}_{n}\left(\beta_{1}\right) \end{array}\right], \quad \hat{Y}\left(\beta_{1}\right)=\left[\begin{array}{c} \widetilde{Y}_{1}\left(\beta_{1}\right) \\ \vdots \\ \widetilde{Y}_{n}\left(\beta_{1}\right) \end{array}\right]$
此时我们求出 $\beta_2$ ,再将 $\beta_2$ 作为初始值依次迭代计算，得到关于向量参数 $\beta_i$ 的一个序列，当且仅当
$||\beta^{(k+1)}-\beta^{(k)}||<\delta$
其中 $\delta>0$ 为事先预定的绝对误差。不难得到，参数 $\beta$ 满足递推关系
$\begin{aligned} \boldsymbol{\beta}_{n+1} &=\left[\widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)\right]^{-1} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime} \widetilde{Y}\left(\boldsymbol{\beta}_{n}\right) \\ &=\left[\widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)\right]^{-1} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime}\left[\boldsymbol{Y}-f\left(\widetilde{X}, \boldsymbol{\beta}_{n}\right)+\widetilde{X}\left(\boldsymbol{\beta}_{n}\right) \boldsymbol{\beta}_{n}\right] \\ &=\boldsymbol{\beta}_{n}+\left[\widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)\right]^{-1} \widetilde{X}\left(\boldsymbol{\beta}_{n}\right)^{\prime}\left[Y-f\left(X, \boldsymbol{\beta}_{n}\right)\right] \end{aligned}$
通过证明，随着样本容量 $n\to\infty$ ,参数 $\beta$ 估计量服从渐进正态分布，即
$\widetilde{\beta} \sim N\left(\beta, \hat{\sigma}^{2}\left[\widetilde{X}(\beta)^{\prime} \widetilde{X}(\beta)\right]^{-1}\right), \hat{\sigma}^{2}=\frac{S(\widetilde{\beta})}{n-1}$

3 非线性最小二乘的实现

在R语言中，可以适用nls函数实现非线性最小二乘法。以C-D函数为例，

设一国产出取决于资本、劳动与全要素的投入，即
$AK^{\alpha}L^{\beta}\mu$
下面通过R代码运行实现对参数 $\alpha,\beta$ 的估计

t = 1:12 #时间设定
Y=c(26.74, 34.81, 44.72, 57.46, 73.84, 88.45, 105.82,126.16, 150.9, 181.6, 204.3, 222.8) #产出序列
K=c(23.66,30.55,38.12,46.77,56.45,67.15,78.92,91.67,105.5, 121.3, 128.6, 132.5) #资本序列
L=c(26, 28, 32, 36, 41, 45, 48, 52, 56, 60, 66, 70) #劳动投入序列
Cdnls <- nls(Y~A*K^a*L^b,start = list(A = 0.1,a = 0.5,b = 0.5)) #非线性最小二乘,start为参数初始值向量
summary(Cdnls)
#-------------------运行结果---------------------------
#Formula: Y ~ A * K^a * L^b
Parameters:
  Estimate Std. Error t value Pr(>|t|)    
A   0.1129     0.0159    7.12  5.6e-05 ***
a   0.6568     0.0652   10.07  3.4e-06 ***
b   1.0298     0.1044    9.86  4.0e-06 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.7 on 9 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 7.55e-06

结果显示，参数 $\alpha = 0.6568$ , $\beta = 1.0298$ 。对比直接取对数的OLS,即估计
$lnA+\alpha lnK+\beta lnL+e$

CDlm <- lm(log(Y)~log(K)+log(L))  #对数形式
summary(CDlm)
#--------------------运行结果--------------
Call:
lm(formula = log(Y) ~ log(K) + log(L))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.02714 -0.00595 -0.00118  0.00764  0.02557 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.0737     0.2355   -8.80  1.0e-05 ***
log(K)        0.6258     0.0916    6.83  7.6e-05 ***
log(L)        1.0379     0.1621    6.40  0.00012 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0173 on 9 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:  0.999 
F-statistic: 9.16e+03 on 2 and 9 DF,  p-value: 1.29e-15

结果显示，参数 $\alpha = 0.6268$ , $\beta = 1.0379$ 。因此，CD函数对数化的结果回归与非线性最小二乘回归的参数基本一致。但一些不能对数化的方程，非线性最小二乘的作用更为明显。考虑真实模型
$y = 2 s in (x) + 4 cos (x)$
接下来我们进行仿真模拟

set.seed(123) #随机种子
x <- seq(1,100,by = 0.1) #1-100，步长为0.1
e <- rnorm(length(x),0,1) #长度为序列x的长度，服从标准正态分布的误差
y <- 2*sin(x)+4*cos(x)+e #实际观测的被解释变量
plot(x,y,type = "o") #打印散点图

nls1 <- nls(y~a*sin(x)+b*cos(x),
            start = list(a = 0,b =0)) #非线性最小二乘，初始值设定为0，0
nls1
#-------------运行结果------------------
Nonlinear regression model
  model: y ~ a * sin(x) + b * cos(x)
   data: parent.frame()
   a    b 
1.92 4.03 
 residual sum-of-squares: 974

Number of iterations to convergence: 1 
Achieved convergence tolerance: 6.73e-10

结果显示估计量 $a = 1.92$ , $b = 4.03$ ，与总体参数 $a = 2, b = 4$ 即为接近