【資料分析 R語言實戰】學習筆記 第六章 引數估計與R實現（上）
6.1點估計及R實現
6.1.1矩估計
R中的解方程函式:
函式及所在包：功能
uniroot()@stats：求解一元（非線性)方程
multiroot()@rootSolve：給定n個(非線性)方程，求解n個根
uniroot.all()@rootSolve：在一個區問內求解一個方程的多個根
BBsolve()@BB：使用Barzilai-Borwein步長求解非線性方程組
uniroot(f,interval, ...,lower = min(interval), upper = max(interval),f.lower = f(lower,...), f.upper = f(upper, ...),extendInt = c("no", "yes","downX", "upX"), check.conv = FALSE,tol =.Machine$double.eps^0.25, maxiter = 1000, trace = 0)
其中f指定所要求解方程的函式:interval是一個數值向量，指定要求解的根的區間範圍:或者用lower和upper分別指定區間的兩個端點;tol表示所需的精度(收斂容忍度):maxiter為最人迭代次數。
如果遇到多元方程的求解，就需要利用rootSolve包的函式multiroot()來解方程組。multiroot()用於對n個非線性方程求解n個根，其要求完整的雅可比矩陣，採用Newton-Raphson方法。其呼叫格式為:
multiroot(f, start, maxiter = 100,
rtol = 1e-6, atol = 1e-8, ctol = 1e-8,
useFortran = TRUE, positive = FALSE,
jacfunc = NULL, jactype = "fullint",
verbose = FALSE, bandup = 1, banddown = 1,
parms = NULL, ...)
f指定所要求解的函式;由於使用的是牛頓迭代法，因而必須通過start給定根的初始值，其中的name屬性還可以標記輸出變數的名稱;maxiter是允許的最大迭代次數;rtol和atol分別為相對誤差和絕對誤差，一般保持預設值即可;ctol也是一個用於控制迭代次數的標量，如果兩次迭代的最大變化值小於ctol,那麼迭代停止，得到方程組的根。
例如，己知某種保險產品在一個保單年度內的損失情況如下所示，其中給出了不同損失次數下的保單數，我們對損失次數的分佈進行估計。已知分佈型別是泊松(Poisson ) ,其樣本均值即為引數λ的矩估計。
|
損失次數 |
0 |
1 |
2 |
3 |
4 |
5 |
|
保單數 |
1532 |
581 |
179 |
41 |
10 |
4 |
|
1 2 3 4 |
|
畫圖比較損失次數的估計值和樣本值之間的差別
|
1 2 3 4 5 6 7 |
|
rootSolve包的函式multiroot()用於解方程組：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
|
驗證一下：
|
1 2 3 |
|
6.1.2極大似然估計
R中計算極值的函式（stats包）
optimize( ) 計算單引數分佈的極人似然估計值
optim() 計算多個引數分佈的極大似然估計值
nlm() 計算非線性函式的最小值點
nlminb( ) 非線性最小化函式
1.函式optimize()
當分佈只包含一個引數時，我們可以使用R中計算極值的函式optimize()求極大似然估計值。
optimize(f = , interval = , ..., lower = min(interval),upper = max(interval), maximum = FALSE,tol = .Machine$double.eps^0.25)
其中f是似然函式:interval指定引數的取值範圍;lower/upper分別是引數的下界和上界:maximum預設為FALSE，表示求似然函式的極小值，若為TRUE則求極大值:tol表示計算的精度。
2.函式optim()和nlm()
當分佈包含多個引數時，用函式optim()或nlm()計算似然函式的極大值點。
optim(par, fn, gr = NULL, ...,method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent"),lower = -Inf, upper = Inf,control = list(), hessian = FALSE)
par設定引數的初始值；fn為似然函式；method提供了5種計算極值的方法
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)
nlm是非線性最小化函式，僅使用牛頓一拉夫遜演算法，通過迭代計算函式的最小值點。一般只布要對前兩個引數進行設定:f是需要最小化的函式:P設定引數初始值。
3.函式nlminb()
在實際應用中，上面這三個基本函式在遇到資料量較大或分佈較複雜的計算時，就需要使用優化函式nlminb()
nlminb(start, objective, gradient = NULL, hessian = NULL, ...,scale = 1, control = list(), lower = -Inf, upper = Inf)
引數start是數值向量，用於設定引數的初始值;objective指定要優化的函式:gradient和hess用於設定對數似然的梯度，通常採用預設狀態;control是一個控制引數的列表:lower和upper設定引數的下限和上限，如果未指定，則假設所有引數都不受約束。
例：
|
1 2 3 4 5 6 7 8 9 10 |
|
猜測分佈是兩個正態分佈的混合，需要估計出函式中的5個引數：p、μ1、μ2、σ1、σ2。
在R中編寫對數似然函式時，5個引數都存放在向量para中，由於nlminb()是計算極小值的，因此函式function中最後返回的是對數似然函式的相反數。
|
1 2 3 4 5 6 7 8 |
|
做引數估計，使用nlminb()之前最大的要點是確定初始值，初始值越接近真實值，計算的結果才能越精確。我們猜想資料的分佈是兩個正態的混合，概率P直接用0.5做初值即可。通過直方圖中兩個峰對應的x軸數值(大概為50和80>，就可以將初值設定為μ1和μ2。而概率P處於((0,1)區間內，引數σ1，σ2是正態分佈的標準差，必須大於0，所以通過lower和upper兩個引數進行一定的約束。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
|
(2)使用極大似然估計函式maxLik()計算
程式包maxLik中同名的函式maxLik()可以直接計算極大似然估計值，呼叫格式如下:
maxLik(logLik, grad = NULL, hess = NULL, start, method,constraints=NULL, ...)
logLik是對數似然函式，grad和hess用於設定對數似然的梯度，通常不需要進行設定，採用預設值NULL即可;start是一個數值向量，設定引數的初始值;method選擇求解最大化的方法，包括“牛頓-拉夫遜”、"BFGS". "BFGSR", "BHHH","SANK”和“Nelder-Mead"，如果不設定，將自動選擇一個合適的方法;constraints指定對似然估計的約束。
例：
採用兩引數的負二項分佈做極大似然估計，具體說明離散分佈的擬合：
編寫R程式時首先要寫出對數似然函式loglik，用到R中的負二項函式dnbinom()，它的引數是r、p。如果要估計β的值，應當轉換一下形式。
|
1 2 3 4 5 6 7 8 9 10 |
|
通過圖形來觀察估計的效果，比較損失次數的樣本值和估計值：
|
1 2 3 4 5 6 |
|
可以看出，負二項分佈的極大似然估計效果非常好，估計值與樣木值幾乎完全重合，可以得出結論，損失次數服從負二項分佈。
6.2單正態總體的區間估計
6.2.1均值μ的區間估計
(1 )σ2已知
R中沒有計算方差己知時均值置信區間的內建函式，需要自己編寫：
conf.int=function(x,sigma,alpha){
mean=mean(x)
n=length(x)
z=qnorm(1-alpha/2,mean=0,sd=1,lower.tail=TRUE)
c(mean-sigma*z/sqrt(n),mean+sigma*z/sqrt(n))
}
其中x為資料樣本;sigma是已知總體的標準差;alpha表示顯著性水平。通常我們作區間估計時，都會估計出雙側的置信區間，因為它為待估引數提供了上下限兩個參考值。但如果要估計單.側的置信區間，理論上與雙側相同，只需要使用標準正態分佈的α分位點即可，編寫函式時也做同樣變動即可。
現在基本統計和資料分析程式包BSDA (Basic Statisticsand Data Analysis )中己經提供了函式z.test()，它可以對基於正態分佈的單樣本和雙樣本進行假設檢驗、區間估計，其使用方法如下:
z.test(x, y = NULL, alternative = "two.sided", mu = 0, sigma.x = NULL,sigma.y = NULL, conf.level = 0.95)
其中，x和Y為數值向量，預設y=NULL，即進行單樣本的假設檢驗;alternative用於指定所求置信區間的型別，預設為two.sided，表示求雙尾的置信區間，若為less則求置信上限，為greater求置信卜限;mu表示均值，它僅在假設檢驗中起作用，預設為0; sigma.x和sigma.y分別指定兩個樣本總體的標準差:conf.level指定區間估計時的置信水平。
程式包UsingR中的函式simple.z.test()，它專門用於對方差己知的樣本均值進行區間估計，與z.test()的不同點在於它只能進行置信區間估計，而不能實現Z檢驗。simple.z.test()
的使用方法如下:
simple.z.test (x,sigma, conf.level=0.95)
其中，x是資料向量:sigma是己知的總體標準差;conf.level指定區間估計的置信度，預設
為95% 。
例：
從均值為10、標準差為2的總體中抽取20個樣本，因此這是一個方差己知
的正態分佈樣本。計算置信水平為95%時x的置信區間，首先呼叫自行編寫的函式conf.int()：
|
1 2 3 4 5 6 7 8 9 10 |
|
用函式z.test（）也可以直接得到這一結果：
|
1 2 3 4 5 |
|
simple.z.test()，可以直接得到區間估計結果：
|
1 2 3 |
|
三種方法的結果均顯示，該樣本的95%置信區間為[8.42, 10.17]
(2 )σ2未知
總體方差未知時，用t分佈的統計量來替代z，方差也要由樣本方差s2代替
t.test(x, y = NULL,alternative = c("two.sided", "less", "greater"),mu = 0, paired = FALSE, var.equal = FALSE,conf.level = 0.95, ...)
其中，x為樣本資料;若x和Y同時輸入，則做雙樣本t檢驗;alternative用於指定所求置信區間的型別，預設為two.sided，表示求雙尾的置信區間，若為less則求置信上限，為greater求置信下限;mu表示均值，其僅在假設檢驗中起作用，預設為0.
仍使用上例中的向量x，假設總體方差未知時，用函式t.test()計算置信區間後:
|
1 2 3 4 5 6 7 8 9 10 11 |
|
如果只要區間估計的結果，則用符號“$”選取conf.int的內容:
|
1 2 3 4 |
|
6.2.2方差σ2的區間估計
(1)μ已知
(2) μ未知
在R中沒有直接計算方差的置信區間的函式，我們可以把上面兩種情況寫在一個函式裡，通過一個if語句進行判斷，只要是方差的區間估計，都呼叫這個函式即可。在R中寫函式時，引數可以事先設定一個初值，例如設mu=Inf，代表均值未知的情況，呼叫函式時如果沒有特殊說明mu的值，將按照均值未知的方法計算;如果均值己知，在呼叫函式時應該對mu重新賦值。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
|
計算得到總體方差的置信區間為【5.35,39.5]，置信水平是95%