（一）极大似然估计

在R中作极大似然估计，主要就是定义似然后函数，然后再用nlminb函数对参数进行估计。
对于r语言也可以直接使用RLRsim包做计算

1. 数据与模型
我们要使用的数据来自于“MASS”包中的geyser数据。先把数据调出来，看看它长什么样子。
 geyser
    waiting duration
1 80 4.0166667
2 71 2.1500000
3 57 4.0000000
4 80 4.0000000
5 75 4.0000000
……
该数据采集自美国黄石公园内的一个名叫Old Faithful 的喷泉。“waiting”就是喷泉两次喷发的间隔时间，“duration”当然就是指每次喷发的持续时间。在这里，我们只用到“waiting”数据，为了简单一点，可以使用attach()函数。
attach(geyser)
2. 绘制出数据的频率分布直方图：
 hist(waiting)
从图中可以看出，其分布是两个正态分布的混合。

具体的代码实现为：

require(MASS)
attach(geyser)
head(geyser)
hist(waiting)
lines(density(waiting))


第一部分 #定义log-likelihood函数
LL<-function(params,data)
  {#参数"params"是一个向量，依次包含了五个参数：p,mu1,sigma1,
  #mu2,sigma2.
  #参数"data"，是观测数据。
t1<-dnorm(data,params[2],params[3])
t2<-dnorm(data,params[4],params[5])
  #这里的dnorm()函数是用来生成正态密度函数的。
f<-params[1]*t1+(1-params[1])*t2
  #混合密度函数
ll<-sum(log(f))
  #log-likelihood函数
return(-ll)
  #nlminb()函数是最小化一个函数的值，但我们是要最大化log-
  #likeilhood函数，所以需要在“ll”前加个“-”号。
  }


第二部分 参数估计
  #用hist函数找出初始值
hist(waiting,freq=F)
lines(density(waiting))
  #拟合函数####optim####
geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=waiting,lower=c(0.0001,-Inf,0.0001,-Inf,-Inf,0.0001),upper=c(0.9999,Inf,Inf,Inf,Inf))
  #初始值为p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10
  #LL是被最小化的函数。
  #data是拟合用的数据
  #lower和upper分别指定参数的上界和下界。




第三部分： 估计结果
  #查看拟合的参数
geyser.res$par
  #拟合的效果
X<-seq(40,120,length=100)
  #读出估计的参数
p<-geyser.res$par[1]
mu1<-geyser.res$par[2]
sig1<-geyser.res$par[3]
mu2<-geyser.res$par[4]
sig2<-geyser.res$par[5]
  #将估计的参数函数代入原密度函数。
f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2)
  #作出数据的直方图
hist(waiting,probability=T,col=0,ylab="Density",ylim=c(0,0.04),xlab="Eruption waiting times")
  #画出拟合的曲线
lines(X,f)
lines(density(waiting),col='red')
  #补充密度曲线做对比
detach()

（二）置信区间

用R语言求置信区间是很方便的，而且很灵活，至少比spss好多了。
如果你要求的只是95%的置信度的话，那么用一个很简单的命令就可以实现了
首先，输入da=c(你的数据，用英文逗号分割），然后t.test(da)，运行就能得到结果了。
我的数据是newbomb <- c(28,26,33,24,34,-44,27,16,40,-2,29,22,24,21,25,30,23,29,31,19)
t.test(newbomb)得到的结果如下

如果要求任意置信度下的置信区间的话，就需要自己编一个函数了。
当然，有两点要记住的，置信区间的计算在知道方差和不知道方差的情况下，计算公式是不一样的。
下面做一个两种情况下都可以用的函数。
confint<-function(x,sigma=-1,alpha=0.05)
  {
      n<-length(x)
      xb<-mean(x)
      if(sigma>=0)
          {
             tmp<-sigma/sqrt(n)*qnorm(1-alpha/2);df<-n
           }
       else{
           tmp<-sd(x)/sqrt(n)*qt(1-alpha/2,n-1);df<- n-1
           }
       data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp)
   }
这个函数的使用：
   如果不知道方差，则confint（x,alpha） 知道方差，则confint（x,sigma,alpha)
这样就能计算出结果了

转自：https://www.douban.com/note/420549640/?type=like