社区首页 >问答首页 >如何在R中编写最大似然例程？

问如何在R中编写最大似然例程？
EN

Stack Overflow用户

提问于 2012-07-10 17:19:56

回答 2查看 1.7K关注 0票数 3

我想在我创建的数据样本上运行一些最大似然代码。这就是我到目前为止所知道的：

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个样本中每个样本的平均值和标准差时遇到了一些问题，我修改了应用函数以尝试适应这一点，但还没有起到任何作用。有什么想法吗？

CDN&音视频通信出海专场

提供游戏出海、社交娱乐等方案，助力企业抢占出海市场

回答 2

Stack Overflow用户

回答已采纳

发布于 2012-07-10 17: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.89506586

票数 5

Stack Overflow用户

发布于 2012-07-11 23: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