如何在R中编写最大似然例程？
如何在R中编写最大似然例程？
提问于 2012-07-11 01:19:56
我想在我创建的数据样本上运行一些最大似然代码。这就是我到目前为止所知道的：
library("maxLik")
data <- replicate(20, rnorm(100))
logLikFun <- function(param) {
mu <- param[1]
sigma <- param[2]
sum(dnorm(data, mean = mu, sd = sigma, log = TRUE))
}
mle <- maxLik(logLik = logLikFun, start = c(mu = 0, sigma = 1))
summary(mle)我在提取20个样本中每个样本的平均值和标准差时遇到了一些问题，我修改了应用函数以尝试适应这一点，但还没有起到任何作用。有什么想法吗？
回答 2
Stack Overflow用户
回答已采纳
发布于 2012-07-11 01:30:47
创建一个函数(在本例中为find.mle)，该函数接受一个数据矢量并基于该矢量计算最大似然估计，然后使用apply将其应用于data的列
library("maxLik")
data <- replicate(20, rnorm(100))
find.mle = function(d) {
logLikFun <- function(param) {
mu <- param[1]
sigma <- param[2]
sum(dnorm(d, mean = mu, sd = sigma, log = TRUE))
}
maxLik(logLik = logLikFun, start = c(mu = 0, sigma = 1))$estimate
}
mles = apply(data, 2, find.mle)这将为您提供一个2x20矩阵，其中包含您的估计：
> mles
[,1] [,2] [,3] [,4] [,5] [,6]
mu 0.03675611 0.1129927 -0.06499549 0.04651673 0.06593217 -0.08753828
sigma 0.93497523 0.9817961 0.84734600 0.93139761 1.01083924 1.04114752
[,7] [,8] [,9] [,10] [,11] [,12]
mu 0.1629807 0.01665411 0.2306688 -0.02147982 0.07723695 0.009476477
sigma 1.0428713 1.01658241 1.0073277 0.99781761 0.99327722 0.983356049
[,13] [,14] [,15] [,16] [,17] [,18]
mu 0.06524147 0.02442983 -0.1305258 -0.1050299 0.1449996 0.1172218
sigma 1.04004799 0.89963009 0.9979824 1.0227063 0.9319562 0.9916734
[,19] [,20]
mu -0.1288296 -0.05769467
sigma 0.9975368 0.89506586Stack Overflow用户
发布于 2012-07-12 07:41:39
我真的认为不需要编写任何函数来获得均值和sd的最大似然(以下简称ML)估计值。如果X是正态随机变量，则总体均值和sd的ML估计器是样本均值和样本sd，我们知道样本均值是总体均值的无偏估计器，但方差的ML估计器是有偏的(向下)，因为方差的分母是n而不是n-1。
所以R计算样本准方差(修正了自由度)，这是无偏估计器，所以它不是最大似然估计器，但我们可以从R估计器获得最大似然估计器，我们只需将其乘以(n-1)(1/n)，结果将是方差的最大似然估计值，然后应用平方根，然后得到sd的最大似然估计值，但我喜欢简单的东西，只要将sd乘以(n-1)(1/n)，这就是你的答案。有关详细说明，请参阅http://en.wikipedia.org/wiki/Variance上的总体方差和样本方差
现在，您只需在R中执行以下操作：
## Reproducing @ David Robinson code
install.packages('maxLik')
library("maxLik")
set.seed(007) ## making it reproducible
data <- replicate(20, rnorm(100))
find.mle = function(d) {
logLikFun <- function(param) {
mu <- param[1]
sigma <- param[2]
sum(dnorm(d, mean = mu, sd = sigma, log = TRUE))
}
maxLik(logLik = logLikFun, start = c(mu = 0, sigma = 1))$estimate
}
mles = apply(data, 2, find.mle)
apply(data, 2, function(x) c(Mean=mean(x), SD=(n-1)*(1/n)*sd(x))) # my simple answer.
# Comparing results:
> mles
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
mu 0.1386966 0.1304418 -0.03515036 -0.05065659 0.04170382 0.0007424064 -0.07625412
sigma 0.9540009 0.9442371 1.07218240 1.03162817 0.96140925 1.0274500157 0.87450358
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
mu 0.02024026 -0.1732926 0.03401213 -0.1254751 0.05263887 -0.01258275 -0.02843866
sigma 0.98456202 0.9628233 0.95087131 0.9912367 1.01347266 0.99542339 1.03761674
[,15] [,16] [,17] [,18] [,19] [,20]
mu 0.02441331 -0.03021781 0.2170172 0.02271656 -0.04946737 0.115728
sigma 1.03889635 1.02796932 1.0457951 1.07906578 0.93627993 1.009641
> apply(data, 2, function(x) c(Mean=mean(x), SD=(n-1)*(1/n)*sd(x)))
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
Mean 0.1386966 0.1304418 -0.03515036 -0.05065659 0.04170382 0.0007424064 -0.07625412
SD 0.9492189 0.9395041 1.06680802 1.02645707 0.95659012 1.0222998579 0.87012008
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
Mean 0.02024026 -0.1732926 0.03401213 -0.1254751 0.05263887 -0.01258275 -0.02843866
SD 0.97962684 0.9579971 0.94610501 0.9862680 1.00839257 0.99043377 1.03241563
[,15] [,16] [,17] [,18] [,19] [,20]
Mean 0.02441331 -0.03021781 0.2170172 0.02271656 -0.04946737 0.115728
SD 1.03368881 1.02281656 1.0405530 1.07365689 0.93158677 1.004580所以如果你只使用一个简单的产品，你可以删除这个函数( @David Robinson写的一个非常好的函数)。这是一个简单的理论统计观点。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/11418534
复制相关文章
相似问题
腾讯云开发者
