如何在R中編寫最大似然程序？

Question

我想在我創建的數據樣本上運行一些最大可能性代碼。這是我到目前爲止有：

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的每一個樣品，我修改了應用功能，以儘量滿足這一點，但什麼也還沒有工作。有任何想法嗎？

Answer 1

創建（在這個例子​​）一個函數，數據的矢量和基於它計算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) 

這會給你一個2×20矩陣，你的估計：

> 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 

Answer 2

我真的覺得有沒有需要編寫任何功能，以獲得最大似然（ML下同）估計的平均值和標準差。如果X是一個正態隨機變量，那麼總體均值和sd的ML估計量就是樣本均值和樣本sd，我們知道樣本均值是總體均值的無偏估計量，但方差的ML估計量是有偏差的（向下），因爲方差的分母是n而不是n-1。

因此R計算樣本準方差（針對自由度進行校正），這是無偏估計量，所以它不是ML估計量，但我們可以從R估計中獲得ML估計量，必須乘以（n-1）（1/n），結果將是方差的ML估計，然後應用平方根和voilá，你將有sd的ML估計，但我喜歡容易的東西所以，只需乘以（n-1）（1/n）的sd，這就是你的答案。有關詳細說明，請參見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 

所以，你可以刪除功能（通過@大衛羅賓遜寫了一個很不錯的功能）如果你只是使用一個簡單的產品。這是一個簡單的理論統計觀點。