我的模型中有以下對數似然，我試圖將其撰寫為 R 中的函式。
我的問題來了，因為我不知道如何根據函式撰寫 theta。我已經對此進行了幾次嘗試，如下所示，如果這些提示/建議接近正確，我們將不勝感激。
theta 寫為 theta 的函式
#my likelihood function
mylikelihood = function(beta) {
#log-likelihood
result = sum(log(dengue$cases theta 1 / dengue$cases))
sum(theta*log(theta / theta exp(beta[1] beta[2]*dengue$time)))
sum(theta * log(exp(beta[1] beta[2]*dengue$time / dengue$cases exp(beta[1] beta[2]*dengue$time))))
#return negative log-likelihood
return(-result)
}
我的下一次嘗試將 thetas 替換為我的資料集中的 Xi，這里是 (dengue$time)
#my likelihood function attempt 2
mylikelihood = function(beta) {
#log-likelihood
result = sum((log(dengue$Cases dengue$Time 1 / dengue$Cases)))
sum(dengue$Time*log(dengue$time / dengue$Time exp(beta[1] beta[2]*dengue$Time)))
sum(dengue$Cases * log(exp(beta[1] beta[2]*dengue$Time / dengue$Cases
exp(beta[1] beta[2]*dengue$Time))))
#return negative log-likelihood
return(-result)
}
資料
head(dengue)
Cases Week Time
1 148 36 1
2 275 37 2
3 205 38 3
4 133 39 4
5 123 40 5
6 138 41 6
這些中的任何一個都接近正確，如果不是，我哪里出錯了？
更新了對數似然的來源；
該模型;
具有均值 μ 和色散引數 θ 的負二項分布具有 pmf；
uj5u.com熱心網友回復：
基本問題是您必須同時傳遞beta（問題的一個組成部分的截距和斜率）以及theta作為單個引數向量的一部分。我修復了括號放置的其他問題，我稍微重新組織了運算式。
您的代碼中有幾個更重要的錯誤。
- 第一項不是分數；它是一個二項式系數。（即，您應該使用
lchoose()，如下所示） - 您在第一學期將 1 更改為 -1。
nll <- function(pars) {
beta <- pars[1:2]
theta <- pars[3]
##log-likelihood
yi <- dengue$Cases
xi <- dengue$Time
ri <- exp(beta[1] beta[2]*xi)
result <- sum(lchoose(yi theta - 1,yi))
sum(theta*log(theta / (theta ri)))
sum(yi * log(ri/(theta ri)))
##return negative log-likelihood
return(-result)
}
讀取資料
dengue <- read.table(row.names = 1, header = TRUE, text = "
Cases Week Time
1 148 36 1
2 275 37 2
3 205 38 3
4 133 39 4
5 123 40 5
6 138 41 6
")
配件
猜測 (1,1,1) 的起始引數有點危險——了解引數的含義并猜測生物學上合理的值會更有意義——但這似乎沒問題。
nll(c(1,1,1))
optim(par = c(1,1,1), nll)
由于我們沒有將 theta 限制為正數，因此我們會收到一些關于記錄負數日志的警告，但這些可能是無害的（例如，請參見此處）
備擇方案
R 有很多內置的機器來擬合負二項式模型（我應該知道你在做什么！）
MASS::glm.nb自動為您設定一切，您只需指定預測變數（它使用對數鏈接并添加截距，因此指定~Time將使平均值等于exp(beta0 beta1*Time)）。
library(MASS)
glm.nb(Cases ~ Time, data = dengue)
bbmle自動化程度較低，但更靈活（這里我適合theta對數刻度以避免嘗試任何負值）
library(bbmle)
mle2(Cases ~ dnbinom(mu = exp(logmu), size = exp(logtheta)),
parameters = list(logmu ~ Time),
data = dengue,
start = list(logmu = 0, logtheta = 0))
所有這三種方法（校正的負對數似然函式 optim, MASS::glm.nb, bbmle::mle2）都給出了相同的結果。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/484902.html
上一篇：基于向量長度的遞回函式
下一篇：如何重用每個函式的條件檢查？