眾所周知，R語言是個不錯的統計軟體。今天分享一下利用R語言做點估計的內容。主要有：矩估計、極大似然估計、EM演算法、最小二乘估計、刀切法（Jackknife）、自助法（Bootstrap）的相關內容。

          點估計是引數估計的一個組成部分。有許多的估計方法與估計理論，具體內容可以參見lehmann的《點估計理論》（推薦第一版，第二版直接從UMVU估計開始的）

一、矩估計

         對於隨機變數來說，矩是其最廣泛，最常用的數字特徵，母體的各階矩一般與的分佈中所含的未知引數有關，有的甚至就等於未知引數。由辛欽大數定律知，簡單隨機子樣的子樣原點矩依概率收斂到相應的母體原點矩

          因為不同的分佈有著不同的引數，所以在R的基本包中並沒有給出現成的函式，我們通常使用人機互動的辦法處理矩估計的問題，當然也可以自己編寫一些函式。

首先，來看看R中給出的一些基本分佈，如下表：

           雖然R中基本包中沒有現成求各階矩的函式，但是對於給出的樣本，R可以求出其平均值（函式：mean），方差（var），標準差（sd），在fBasics包中還提供了計算偏度的函式

          注：在actuar包中，函式emm（）可以計算樣本的任意階原點矩。但在引數估計時需要注意到原點矩的存在性

例如我們來看看正態分佈N（0,1）的矩估計效果。

> x<-rnorm(100)     #產生N（0,1）的100個隨機數

> mu<-mean(x)     #對N(mu,sigma)中的mu做矩估計

> sigma<-var(x)    #這裡的var並不是樣本方差的計算函式，而是修正的樣本方差，其實也就是

> mu

[1] -0.1595923

> sigma

[1] 1.092255

可以看出，矩估計的效果還是可以的。

我們再來看一個矩估計的例子：設總體X服從二項分佈B(k,p)，X1,X2,…,Xn，是總體的一個樣本。K,p為未知引數。那麼k,p的矩估計滿足方程：kp – M₁ = 0, kp(1 − p) −M2 = 0.

我們可以編寫函式：

moment <-function(p){

f<-c(p[1]*p[2]-M1,p[1]*p[2]-p[1]*p[2]^2-M2)

J<-matrix(c(p[2],p[1],p[2]-p[2]^2,p[1]-2*p[1]*p[2]),nrow=2,byrow=T)#jicobi矩陣

list(f=f,J=J)

}# p[2]=p, p[1]=k,

檢驗程式

 x<-rbinom(100, 20, 0.7); n<-length(x)

 M1<-mean(x);M2<-(n-1)/n*var(x)

 p<-c(10,0.5)

Newtons(moment, p)$root #是用牛頓法解方程的程式，見附件1

執行結果為：

 [1]22.973841  0.605036

可以得到k,p的數值解

二、極大似然估計（MLE）

            極大似然估計的基本思想是：基於樣本的資訊引數的合理估計量是產生獲得樣本的最大概率的引數值。值得一提的是：極大似然估計具有不變性，這也為我們求一些奇怪的引數提供了便利。

           在單引數場合，我們可以使用函式optimize（）來求極大似然估計值，函式的介紹如下：

optimize(f = , interval = ,  ..., lower = min(interval),

        upper = max(interval), maximum = FALSE,

        tol = .Machine$double.eps^0.25)

           例如我們來處理Poisson分佈引數lambda的MLE。

         設X1,X2,…,X100為來自P（lambda）的獨立同分布樣本，那麼似然函式為：

L（lambda，x）=lambda^(x1+x2+…+x100)*exp(10*lambda)/(gamma(x1+1)…gamma(x100+1))

          這裡涉及到的就是一個似然函式的選擇問題：是直接使用似然函式還是使用對數似然函式，為了說明這個問題，我們可以看這樣一段R程式：

> x<-rpois(100,2)

> sum(x)

[1] 215

> ga(x)#這是一個求解gamma(x1+1)…gamma(x100+1)的函式，用gamma函式求階乘是為了提高計算效率（原始碼見附1）

[1] 1.580298e+51

> f<-function(lambda)lambda^(215)*exp(10*lambda)/(1.580298*10^51)#這裡有一些magic number + hard code 的嫌疑，其實用ga(x)帶入，在函式引數中多加一個x就好

> optimize(f,c(1,3),maximum=T)

$maximum

[1] 2.999959

$objective

[1] 2.568691e+64

> fun<-function(lambda)-100*lambda+215*log(lambda)-log(1.580298*10^51)

> optimize(fun,c(1,3),maximum=T)

$maximum

[1] 2.149984  #MLE

$objective

[1] -168.3139

              為什麼會有這樣的差別？這個源於函式optimize，這個函式本質上就是求一個函式的最大值以及取最大值時的自變數。但是這裡對函式的穩定性是有要求的，取對數無疑增加了函式的穩定性，求極值才會合理。這也就是當你擴大了MLE存在區間時warning會出現的原因。當然，限定範圍時，MLE會在邊界取到，但是，出現邊界時，我們需要更多的資訊去判斷它。這個例子也說明多數情況下利用對數似然估計要優於似然函式。

            在多元ML估計時，你能用的函式將變為optim,nlm,nlminb它們的呼叫格式分別為：

optim(par, fn, gr = NULL, ...,

      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),

      lower = -Inf, upper = Inf,

      control = list(), hessian = FALSE)

nlm(f, p, ..., hessian = FALSE, typsize = rep(1, length(p)),

    fscale = 1, print.level = 0, ndigit = 12, gradtol = 1e-6,

    stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000),

    steptol = 1e-6, iterlim = 100, check.analyticals = TRUE)

nlminb(start, objective, gradient = NULL, hessian = NULL, ...,

       scale = 1, control = list(), lower = -Inf, upper = Inf)

> x<-rnorm(1000,2,6)  #6是標準差，而我們估計的是方差

> ll<-function(theta,data){

+ n<-length(data)

+ ll<--0.5*n*log(2*pi)-0.5*n*log(theta[2])-1/2/theta[2]*sum((data-theta[1])^2)

+ return(-ll)

+ }

>nlminb(c(0.5,2),ll,data=x,lower=c(-100,0),upper=c(100,100)) $par

  [1]  1.984345  38.926692

看看結果估計的還是不錯的，可以利用函式mean，var驗證對正態分佈而言，矩估計與MLE是一致.

然而這裡還有一些沒有解決的問題，比如使用nlminb初始值的選取。希望閱讀到這的朋友給出些建議。

最後指出在stata4，maxLik等擴充套件包中有更多關於mle的東西，你可以通過檢視幫助文件來學習它。

附1：輔助程式程式碼

Newtons<-function (fun, x, ep=1e-5,it_max=100){

 index<-0; k<-1

 while (k<=it_max){

 x1 <- x; obj <- fun(x);

 x <- x - solve(obj$J, obj$f);

 norm <- sqrt((x-x1) %*% (x-x1))

 if (norm<ep){

 index<-1; break

 }

 k<-k+1

 }

 obj <- fun(x);

 list(root=x, it=k, index=index, FunVal= obj$f)

 }

ga<-function(x){
ga<-1
for(i in 1:length(x)){
ga<-ga*gamma(x[i]+1)
}
ga
}