#取出MASS包这中的数据data(geyser,package="MASS")head(geyser)attach(geyser)par(bg=‘lemonchiffon‘)hist(waiting,freq=F,col="lightcoral")#freq=F要加上，否则就无法添加线了lines(density(waiting),lwd=2,col="cadetblue4")#根据图像，我们认为其在前后分别是两个正态分布函数的组合#定义 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])#f是概率密度函数 f<-params[1]*t1+(1-params[1])*t2#混合密度函数 ll<-sum(log(f))#log‐likelihood 函数return(-ll)#nlminb()函数是最小化一个函数的值，#但我们是要最大化 log‐likeilhood 函数#所以需要在“ ll”前加个“ ‐”号。}#估计函数####optim##### debugonce(nlminb)geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=waiting, lower=c(0.0001,-Inf,0.0001,-Inf,0.0001), upper=c(0.9999,Inf,Inf,Inf,Inf))#初始值为 p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10#初始值也会被传递给LL#LL 是被最小化的函数。#data 是估计用的数据(传递给我们的LL)#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=‘lightpink3‘, ylab="Density",ylim=c(0,0.04), xlab="Eruption waiting times")#画出拟合的曲线lines(X,f,col=‘lightskyblue3‘,lwd=2)detach(geyser)function (start,objective, gradient = NULL, hessian = NULL,..., scale =1, control = list(), lower =-Inf, upper =Inf){ par <- setNames(as.double(start), names(start)) n <- length(par) iv <- integer(78+3* n) v <- double(130+(n *(n +27))/2).Call(C_port_ivset,2, iv, v)if(length(control)){ nms <- names(control)if(!is.list(control)||is.null(nms)) stop("‘control‘ argument must be a named list") pos <- pmatch(nms, names(port_cpos))if(any(nap <-is.na(pos))){ warning(sprintf(ngettext(length(nap),"unrecognized control element named %s ignored","unrecognized control elements named %s ignored"), paste(sQuote(nms[nap]), collapse =", ")), domain = NA) pos <- pos[!nap] control <- control[!nap]} ivpars <- pos <=4 vpars <-!ivparsif(any(ivpars)) iv[port_cpos[pos[ivpars]]]<-as.integer(unlist(control[ivpars]))if(any(vpars)) v[port_cpos[pos[vpars]]]<-as.double(unlist(control[vpars]))} obj <- quote(objective(.par,...)) rho <- new.env(parent = environment()) assign(".par", par, envir = rho) grad <- hess <- low <- upp <- NULLif(!is.null(gradient)){ grad <- quote(gradient(.par,...))if(!is.null(hessian)){if(is.logical(hessian)) stop("logical ‘hessian‘ argument not allowed. See documentation.") hess <- quote(hessian(.par,...))}}if(any(lower !=-Inf)|| any(upper !=Inf)){ low <- rep_len(as.double(lower), length(par)) upp <- rep_len(as.double(upper), length(par))}else low <- upp <- numeric().Call(C_port_nlminb, obj, grad, hess, rho, low, upp, d = rep_len(as.double(scale), length(par)), iv, v) iv1 <- iv[1L] list(par = get(".par", envir = rho), objective = v[10L], convergence =(if(iv1 %in%3L:6L)0Lelse1L), iterations = iv[31L], evaluations = c(`function`= iv[6L], gradient = iv[30L]), message =if(19<= iv1 && iv1 <=43){if(any(B <- iv1 == port_cpos)) sprintf("‘control‘ component ‘%s‘ = %g, is out of range", names(port_cpos)[B], v[iv1])else sprintf("V[IV[1]] = V[%d] = %g is out of range (see PORT docu.)", iv1, v[iv1])}else port_msg(iv1))}
原文：http://www.cnblogs.com/xuanlvshu/p/6354033.html