這篇文章給了我一個啟發，我們可以自己用已知分布的密度函數進行組合，然后構建一個新的密度函數啦，然后用極大似然估計MLE進行估計。
代碼和結果演示
代碼：
#取出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)
調試的說明
nlminb函數:
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))
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 <- !ivpars if (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 <- NULL if (!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) 0L else 1L), 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))}
par最先獲得了初始值
obj <- quote(objective(.par, ...))
讓obj表示一個名為
objective，接受形參.par和...的函數，即我們傳入的對數似然函數
rho <- new.env(parent = environment())
創建新環境
assign(".par", par, envir = rho)
將
par賦給
rho環境中的
.par變量
.Call(C_port_nlminb,
obj, grad, hess,
rho,
low,
upp, d = rep_len(as.double(scale),
length(par)), iv, v)
C_port_nlminb在名為stats.dll的包里，我用reflector和VS都沒能看到什么東西,唉，最關鍵的東西是缺失的。
此時rho包含了參數的初始值，而obj接受.par和多出來的data，即參數初始值和數據。
在
nlminb
函數的幫助文檔里提到：
...
Further arguments to be supplied to objective.
即在這里是傳遞給對數似然函數的更多參數，這里是data，之所以LL的第一個參數不需要傳入是因為我們在源碼看到其是.par，也可以在幫助文檔看到:
objective
Function to be minimized. Must return a scalar value.
The first argument to objective is the vector of parameters to be optimized, whose initial values are supplied through start(start 參數).
Further arguments (fixed during the course of the optimization) to objective may be specified as well (see ...(參數)).
---------------------------------------
參考：
原文檔：
附件列表