上文用模拟数据演示了双边随机前沿分析（SFA2T）的建模和参数估计过程，并给出了简单的极大似然估计程序。本文我们用一组实际数据来继续研究SFA2T模型的MLE参数估计。
实验数据包含13个自变量，除3个二值型自变量外，均为数值型，因变量为某种农产品的种植产出率。首先需对实验数据进行归一化，防止对数似然函数发生异常或出现溢出。然后将归一化的实验数据代入上文的程序，得到结果如下。从结果来看，参数估计的结果比较凑合，应该还有改进空间。
Maximum Likelihood estimation
Nelder-Mead maximization, 9709 iterations
Return code 0: successful convergence
Log-Likelihood: -79.11712
16 free parameters
Estimates:
Estimate Std. error t value Pr(> t)
par1 0.224194 0.025818 8.683 < 2e-16 ***
par2 0.298188 0.023763 12.548 < 2e-16 ***
par3 0.052398 0.031010 1.690 0.09108 .
par4 0.183463 0.056721 3.234 0.00122 **
par5 -0.008597 0.024388 -0.353 0.72445
par6 -0.045980 0.014175 -3.244 0.00118 **
par7 0.009187 0.019488 0.471 0.63735
par8 -0.009617 0.019852 -0.484 0.62807
par9 0.010799 0.017797 0.607 0.54399
par10 0.005012 0.019052 0.263 0.79249
par11 0.021041 0.019748 1.065 0.28668
par12 0.017241 0.017623 0.978 0.32792
par13 -0.015531 0.015802 -0.983 0.32566
par14 -0.100216 0.045260 -2.214 0.02681 *
par15 -0.012184 0.035846 -0.340 0.73393
par16 -0.142891 0.027261 -5.242 1.59e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------
既然MLE本质上是求解一个非线性最优化问题，我们考虑改造似然函数，令其包含更多的信息。
众所周知，对数似然函数在某个参数组合上的hessian矩阵的逆矩阵就是该参数组合的协方差矩阵cov，因此，如果cov对角线元素为负，或者cov是奇异矩阵，表明该参数组合不可行，反之如果参数组合可行，则可以从cov算出p值，从而知道该参数的显著度。
因此除了对数似然值外，还可以将hessian矩阵的形态、显著度等信息体现在优化过程中。改造后的函数如下，其中nas是hessian矩阵的逆矩阵对角线元素中负值的个数，sig指出hessian矩阵的逆矩阵是否近似为奇异矩阵，spv则可以看作参数估计结果的显著度得分（小星星的数目）。
当然，似然函数经过改造后，已经不能称为“似然函数”了，我称其为适应度函数。
llh=function(p){
hess=numDeriv::hessian(nll,p)
if(sum(is.na(hess))>0)
return(-Inf)
eg=abs(eigen(hess,symmetric = T,only.values = T)$values)
sig=(min(eg) <=(1e-12*max(eg)))*200
v=diag(solve(hess))
nas=sum(v<0)*100
if(is.na(nas)|is.na(sig))
return(-Inf)
if(nas==0){
std=sqrt(v)
t=p/std
pv=2*pnorm(-abs(t))
spv=sum(apply(as.array(pv),MARGIN=1,sigscore))
}else{
spv=0
}
ll=ll(p)
res=ll-nas-sig+spv
print(c(res,ll,nas,sig,spv))
return(res)
}
由于计算hessian矩阵和矩阵求逆时间开销都不小，优化过程需时较长。现在对数似然值提高到了-65.9，同时新增两个显著的参数。
summary(a)
--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 278 iterations
Return code 0: successful convergence
Log-Likelihood: -65.87475
16 free parameters
Estimates:
Estimate Std. error t value Pr(> t)
[1,] 0.187664 0.019789 9.483 < 2e-16 ***
[2,] 0.293360 0.022862 12.832 < 2e-16 ***
[3,] 0.021298 0.054563 0.390 0.696286
[4,] 0.168370 0.059778 2.817 0.004853 **
[5,] 0.001071 0.030012 0.036 0.971535
[6,] -0.045803 0.017623 -2.599 0.009349 **
[7,] 0.007218 0.027785 0.260 0.795037
[8,] -0.025107 0.011717 -2.143 0.032126 *
[9,] 0.016953 0.022023 0.770 0.441435
[10,] -0.014796 0.013567 -1.091 0.275455
[11,] -0.011237 0.015191 -0.740 0.459464
[12,] 0.015300 0.014001 1.093 0.274507
[13,] 0.004043 0.016861 0.240 0.810515
[14,] -0.145904 0.042235 -3.455 0.000551 ***
[15,] -0.092548 0.024308 -3.807 0.000140 ***
[16,] -0.134274 0.026857 -5.000 5.74e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------
需要注意的是，对数似然函数经历如此改造后，空间结构已经发生了变化，所以实际工作中使用上述方法必须谨慎，当出现对数似然值更差，适应度函数的值更优的情况时，应果断舍弃看起来更好的拟合结果，否则就有勉强提高参数显著性的嫌疑了。
网友评论