补充资料——自己实现极大似然估计(最大似然估计)MLE
2017-01-27 22:59
399 查看
这篇文章给了我一个启发,我们可以自己用已知分布的密度函数进行组合,然后构建一个新的密度函数啦,然后用极大似然估计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)[/code]调试的说明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), v[iv1]) else sprintf("V[IV[1]] = V[%d] = %g is out of range (see PORT docu.)",
iv1, v[iv1])
} else port_msg(iv1))
}[/code]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.[b] 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 ...(参数)).---------------------------------------参考:http://cos.name/2009/07/maximum-likelihood-estimation-in-r/原文档:null
附件列表
相关文章推荐
- 参数估计(Parameter Estimation):频率学派(最大似然估计MLE、最大后验估计MAP)与贝叶斯学派(贝叶斯估计BPE)
- 最大似然估计 (MLE)与 最大后验概率(MAP)在机器学习中的应用
- 最大似然估计 (MLE) 最大后验概率(MAP)
- 最大似然估计 (MLE)与 最大后验概率(MAP)在机器学习中的应用
- 最大似然估计 (MLE) 最大后验概率(MAP)
- 最大似然估计(MLE)和最大后验概率(MAP)
- 最大似然估计 (MLE) 最大后验概率(MAP)
- 详解最大似然估计(MLE)、最大后验概率估计(MAP),以及贝叶斯公式的理解
- 最大似然估计MLE_和_最大后验概率MAP 的区别与联系
- 最大似然估计 (MLE) 最大后验概率(MAP)
- 最大似然估计(MLE)和最大后验概率(MAP)
- 最大似然估计(MLE)和最大后验概率(MAP)
- 最大似然估计(MLE)与最小二乘估计(LSE)的区别
- 【MLE】最大似然估计Maximum Likelihood Estimation
- 参数估计(Parameter Estimation):频率学派(最大似然估计MLE、最大后验估计MAP)与贝叶斯学派(贝叶斯估计BPE)
- 参数估计:最大似然估计MLE
- 最大似然估计 (MLE)与 最大后验概率(MAP)在机器学习中的应用
- 最大似然估计法(MLE)
- 参数估计:最大似然估计MLE
- [最大似然估计 MLE] Codeforces 802DEF Helvetic Coding Contest 2017 D. E. F. Marmots