非线性最小二乘

非线性最小二乘

目录

1 非线性最小二乘估计

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

---

因此，数值分析自然成为解决上述问题的有力武器。考虑一般化的非线性回归问题，设总体回归模型满足
$$Y=f(X,\beta )+\epsilon$$
对应的残差平方和为
$$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,\mathbf{\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)$$
记${\frac{\mathrm{d}f\left(X_i,\beta \right)}{\mathrm{d}\beta }|}_{{\beta }_1}\approx \frac{f\left(X_i,\beta \right)-f(X,\beta )}{\beta -{\beta }_1}$，简记${\tilde{X}}_i\left({\beta }_1\right)={\frac{\mathrm{d}f\left(X_i,\beta \right)}{\mathrm{d}\beta }|}_{{\beta }_1}$，则
$$\begin{array}{rl}S(\beta )& =\sum _{i=1}^n{\left[Y_i-f\left(X_i,{\beta }_1\right)-{\tilde{X}}_i\left({\beta }_1\right)\left(\beta -{\beta }_1\right)\right]}^2\\ & =\sum _{i=1}^n{\left[{\tilde{Y}}_i\left({\beta }_1\right)-X_i\left({\beta }_1\right)\beta \right]}^2\end{array}$$
其中
$${\tilde{Y}}_i\left({\beta }_1\right)=Y_i-f\left(X_i,{\beta }_1\right)+{\tilde{X}}_i\left({\beta }_1\right){\beta }_1$$
给定初始值向量${\beta }_i$，则${\tilde{Y}}_i\left({\beta }_1\right)$与${\tilde{X}}_i\left({\beta }_1\right)$可计算，从而求出最小残差平方和。$S(\beta )$对应的回归方程为
$${\tilde{Y}}_i\left({\beta }_1\right)={\tilde{X}}_i(\beta )\beta +{\epsilon }_i$$
最小二乘估计量为
$${\beta }_2={\left[\tilde{X}{\left({\beta }_1\right)}^{\prime }\tilde{X}\left({\beta }_1\right)\right]}^{-1}\tilde{X}{\left({\beta }_1\right)}^{\prime }\tilde{Y}\left({\beta }_1\right)$$
其中
$$\tilde{X}\left({\beta }_1\right)=\left[\begin{array}{c}{\tilde{X}}_1\left({\beta }_1\right)\\ ⋮\\ {\tilde{X}}_n\left({\beta }_1\right)\end{array}\right],\hat{Y}\left({\beta }_1\right)=\left[\begin{array}{c}{\tilde{Y}}_1\left({\beta }_1\right)\\ ⋮\\ {\tilde{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{array}{rl}{\mathbf{\beta }}_{n+1}& ={\left[\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\tilde{X}\left({\mathbf{\beta }}_n\right)\right]}^{-1}\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\tilde{Y}\left({\mathbf{\beta }}_n\right)\\ & ={\left[\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\tilde{X}\left({\mathbf{\beta }}_n\right)\right]}^{-1}\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\left[\mathbf{Y}-f\left(\tilde{X},{\mathbf{\beta }}_n\right)+\tilde{X}\left({\mathbf{\beta }}_n\right){\mathbf{\beta }}_n\right]\\ & ={\mathbf{\beta }}_n+{\left[\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\tilde{X}\left({\mathbf{\beta }}_n\right)\right]}^{-1}\tilde{X}{\left({\mathbf{\beta }}_n\right)}^{\prime }\left[Y-f\left(X,{\mathbf{\beta }}_n\right)\right]\end{array}$$
通过证明，随着样本容量$n\to \infty$,参数$\beta$估计量服从渐进正态分布，即
$$\tilde{\beta }\sim N\left(\beta ,{\hat{\sigma }}^2{\left[\tilde{X}(\beta {)}^{\prime }\tilde{X}(\beta )\right]}^{-1}\right),{\hat{\sigma }}^2=\frac{S(\tilde{\beta })}{n-1}$$

---

3 非线性最小二乘的实现

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

设一国产出取决于资本、劳动与全要素的投入，即
$$Y=A{K}^{\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,即估计
$$lnY=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=2sin(x)+4cos(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$即为接近

---

-END-

参考文献

王斌会(2015).计量经济学建模及R语言应用[M].北京大学出版社