做参数估计，使用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()，它可以对基于正态分布的单样本和双样本进行假设检验、区间估计

程序包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.17attr(,"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-testdata:  xt = 22.6, df = 19, p-value = 3.407e-15alternative hypothesis: true mean is not equal to 095 percent confidence interval:  8.43 10.15sample estimates:mean of x      9.29

如果只要区间估计的结果，则用符号“$”选取conf.int的内容:

> t.test(x)$conf.int[1]  8.43 10.15attr(,"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%