上文用模擬數(shù)據(jù)演示了雙邊隨機(jī)前沿分析（SFA2T）的建模和參數(shù)估計(jì)過(guò)程，并給出了簡(jiǎn)單的極大似然估計(jì)程序加矛。本文我們用一組實(shí)際數(shù)據(jù)來(lái)繼續(xù)研究SFA2T模型的MLE參數(shù)估計(jì)。
實(shí)驗(yàn)數(shù)據(jù)包含13個(gè)自變量，除3個(gè)二值型自變量外绍弟，均為數(shù)值型锐锣，因變量為某種農(nóng)產(chǎn)品的種植產(chǎn)出率腌闯。首先需對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行歸一化，防止對(duì)數(shù)似然函數(shù)發(fā)生異车胥荆或出現(xiàn)溢出绑嘹。然后將歸一化的實(shí)驗(yàn)數(shù)據(jù)代入上文的程序，得到結(jié)果如下橘茉。從結(jié)果來(lái)看，參數(shù)估計(jì)的結(jié)果比較湊合姨丈，應(yīng)該還有改進(jìn)空間畅卓。
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本質(zhì)上是求解一個(gè)非線性最優(yōu)化問(wèn)題，我們考慮改造似然函數(shù)蟋恬，令其包含更多的信息翁潘。
眾所周知，對(duì)數(shù)似然函數(shù)在某個(gè)參數(shù)組合上的hessian矩陣的逆矩陣就是該參數(shù)組合的協(xié)方差矩陣cov歼争，因此拜马，如果cov對(duì)角線元素為負(fù)，或者cov是奇異矩陣沐绒，表明該參數(shù)組合不可行俩莽，反之如果參數(shù)組合可行，則可以從cov算出p值乔遮，從而知道該參數(shù)的顯著度扮超。
因此除了對(duì)數(shù)似然值外，還可以將hessian矩陣的形態(tài)蹋肮、顯著度等信息體現(xiàn)在優(yōu)化過(guò)程中出刷。改造后的函數(shù)如下，其中nas是hessian矩陣的逆矩陣對(duì)角線元素中負(fù)值的個(gè)數(shù)坯辩，sig指出hessian矩陣的逆矩陣是否近似為奇異矩陣馁龟，spv則可以看作參數(shù)估計(jì)結(jié)果的顯著度得分（小星星的數(shù)目）。
當(dāng)然漆魔，似然函數(shù)經(jīng)過(guò)改造后坷檩，已經(jīng)不能稱為“似然函數(shù)”了，我稱其為適應(yīng)度函數(shù)有送。
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)
}
由于計(jì)算hessian矩陣和矩陣求逆時(shí)間開(kāi)銷(xiāo)都不小淌喻，優(yōu)化過(guò)程需時(shí)較長(zhǎng)。現(xiàn)在對(duì)數(shù)似然值提高到了-65.9雀摘，同時(shí)新增兩個(gè)顯著的參數(shù)裸删。
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
--------------------------------------------
需要注意的是，對(duì)數(shù)似然函數(shù)經(jīng)歷如此改造后阵赠，空間結(jié)構(gòu)已經(jīng)發(fā)生了變化涯塔，所以實(shí)際工作中使用上述方法必須謹(jǐn)慎肌稻，當(dāng)出現(xiàn)對(duì)數(shù)似然值更差，適應(yīng)度函數(shù)的值更優(yōu)的情況時(shí)匕荸，應(yīng)果斷舍棄看起來(lái)更好的擬合結(jié)果爹谭，否則就有勉強(qiáng)提高參數(shù)顯著性的嫌疑了。
