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 |
> num=c(rep(0:5,c(1532,581,179,41,10,4)))#用rep()函數生成樣本，樣本值有。一5的數字構成，函數中的第二個向量對應表示每個數字的重復次數 > lambda=mean(num) > lambda [1] 0.4780571
畫圖比較損失次數的估計值和樣本值之間的差別
> k=0:5
> ppois=dpois(k,lambda)
> poisnum=ppois*length(num)#由poisson分布生成的損失次數
> plot(k,poisnum,ylim=c(0,1600))#畫圖比較，為圖形效果更好，用參數ylim設置縱軸的范圍，最小值為0,最大值要大於樣本的最值，選取1600
> samplenum=as.vector(table(num))#樣本的損失次數
> points(k,samplenum,type="p",col=2)
> legend(4,1000,legend=c("num","poisson"),col=1:2,pch="0")
rootSolve包的函數multiroot()用於解方程組：
> x=c(4,5,4,3,9,9,5,7,9,8,0,3,8,0,8,7,2,1,1,2)
> m1=mean(x)
> m2=var(x)
> model=function(x,m1,m2){}
> model=function(x,m1,m2){
+ c(f1=x[1]+x[2]-2*m1,
+ f2=(x[2]-x[1])^2/12-m2)
+ }
> library(rootSolve)
> multiroot(f=model,start=c(0,10),m1=m1,m2=m2)
$root
[1] -0.7523918 10.2523918
#均勻分布的兩個參數值[0.75, 10.25]
$f.root
f1 f2
-5.153211e-12 1.121688e-09
$iter
[1] 4
$estim.precis
[1] 5.634204e-10
驗證一下：
> m1-sqrt(3*m2);m1+sqrt(3*m2) [1] -0.7523918 [1] 10.25239
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設置參數的下限和上限，如果未指定，則假設所有參數都不受約束。
例：
> library(MASS) > head(geyser,5) waiting duration 1 80 4.016667 2 71 2.150000 3 57 4.000000 4 80 4.000000 5 75 4.000000 > attach(geyser) > hist(waiting,freq=FALSE)#通過直方圖了解數據分布的形態
猜測分布是兩個正態分布的混合，需要估計出函數中的5個參數：p、μ1、μ2、σ1、σ2。
在R中編寫對數似然函數時，5個參數都存放在向量para中，由於nlminb()是計算極小值的，因此函數function中最后返回的是對數似然函數的相反數。
> l1=function(para)
+ {
+ f1=dnorm(waiting,para[2],para[3])
+ f2=dnorm(waiting,para[4],para[5])
+ f=para[1]*f1+(1-para[1])*f2
+ l1=sum(log(f))
+ return(-11)
+ }
做參數估計，使用nlminb()之前最大的要點是確定初始值，初始值越接近真實值，計算的結果才能越精確。我們猜想數據的分布是兩個正態的混合，概率P直接用0.5做初值即可。通過直方圖中兩個峰對應的x軸數值(大概為50和80>，就可以將初值設定為μ1和μ2。而概率P處於((0,1)區間內，參數σ1，σ2是正態分布的標准差，必須大於0，所以通過lower和upper兩個參數進行一定的約束。
> geyser.est=nlminb(c(0.5,50,10,80,10),l1,lower=c(0.0001,-Inf,0.0001,-Inf,0.0001),upper=c(0.9999,Inf,Inf,Inf,Inf)) > options(digits=3) > geyser.est$par [1] 0.308 54.203 4.952 80.360 7.508 > p=geyser.est$par[1] > mu1=geyser.est$par[2];sigma1=geyser.est$par[3] > mu2=geyser.est$par[4];sigma2=geyser.est$par[5] > x=seq(40,120) >#將估計的參凌丈函數代入原密度函數 > f=p*dnorm(x,mu1,sigma1)+(1-p)*dnorm(x,mu2,sigma2) > hist(waiting,freq=F) > lines(x,f)
(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。如果要估計β的值，應當轉換一下形式。
> num=c(rep(0:5,c(1532,581,179,41,10,4)))
> loglik=function(para)
+ {
+ f=dnbinom(num,para[1],1/(1+para[2]))
+ l1=sum(log(f))
+ return(l1)
+ }
> library(maxLik)
> para=maxLik(loglik,start=c(0.5,0.4))$estimate
> r=para[1];beta=para[2]
通過圖形來觀察估計的效果，比較損失次數的樣本值和估計值：
> l=length(num)
> nbinomnum=dnbinom(0:5,r,1/(1+beta))*l;nbinomnum
[1] 1530.12 588.08 170.66 44.17 10.74 2.51
> plot(0:5,nbinomnum,ylim=c(0,1600))
> points(0:5,nbinomnum,type="p",col=2)
> legend(3,1000,legend=c("num","poisson"),col=1:2,lty=1)
可以看出，負二項分布的極大似然估計效果非常好，估計值與樣木值幾乎完全重合，可以得出結論，損失次數服從負二項分布。
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()：
> 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))
+ }
> set.seed(111)
> x=rnorm(20,10,2)
> conf.int(x,2,0.05)
[1] 8.42 10.17
用函數z.test（）也可以直接得到這一結果：
> library(BSDA) > z.test(x,sigma.x=2)$conf.int [1] 8.42 10.17 attr(,"conf.level") [1] 0.95
simple.z.test()，可以直接得到區間估計結果：
> library(UsingR) > simple.z.test(x,2) [1] 8.42 10.17
三種方法的結果均顯示，該樣本的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()計算置信區間后:
> t.test(x)
One Sample t-test
data: x
t = 22.6, df = 19, p-value = 3.407e-15
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
8.43 10.15
sample estimates:
mean of x
9.29
如果只要區間估計的結果，則用符號“$”選取conf.int的內容:
> t.test(x)$conf.int [1] 8.43 10.15 attr(,"conf.level") [1] 0.95
6.2.2方差σ2的區間估計
(1)μ已知
(2) μ未知
在R中沒有直接計算方差的置信區間的函數，我們可以把上面兩種情況寫在一個函數里，通過一個if語句進行判斷，只要是方差的區間估計，都調用這個函數即可。在R中寫函數時，參數可以事先設定一個初值，例如設mu=Inf，代表均值未知的情況，調用函數時如果沒有特殊說明mu的值，將按照均值未知的方法計算;如果均值己知，在調用函數時應該對mu重新賦值。
> var.conf.int=function(x,mu=Inf,alpha){
+ n=length(x)
+ if(mu<Inf){
+ s2=sum((x-mu)^2)/n
+ df=n
+ }
+ else{
+ s2=var(x)
+ df=n-1
+ }
+ c(df*s2/qchisq(1-alpha/2,df),df*s2)/qchisq(alpha/2,df)
+ }
> var.conf.int(x,alpha=0.05)
[1] 5.35 39.50
計算得到總體方差的置信區間為【5.35,39.5]，置信水平是95%
