вычисление логарифмической вероятности для многомерной линейной регрессии с использованием R - PullRequest
Новая
       81

вычисление логарифмической вероятности для многомерной линейной регрессии с использованием R

1 голос
/ 11 июля 2019

Я хочу вычислить логарифмическую вероятность для многомерной линейной регрессии. Я не уверен, является ли этот код верным или нет.

Я рассчитал вероятность записи в журнал, используя функцию dmvnorm в пакете mvtnorm r. 

sdmvn_mle <- function(obj){
  sdmvn_mle_1 <- apply(obj$residuals^2,2,mean)
  sdmvn_mle_2 <- mean(residuals(obj)[,1] * residuals(obj)[,2])
  return(matrix(c(sdmvn_mle_1[1], sdmvn_mle_2, sdmvn_mle_2, sdmvn_mle_1[2]), nrow = 2))
} 


llmvn <- function(obj, sd){
  lr <- c()
  for( i in 1: nrow(obj$fitted.values)){
    lr <- c(lr, mvtnorm::dmvnorm(model.response(model.frame(obj))[i,], mean=fitted(obj)[i,], sigma=sd, log=TRUE))
  }
  return(sum(lr))
}

Y <- as.matrix(mtcars[,c("mpg","disp")])
(mvmod <- lm(Y ~ hp + drat + wt, data=mtcars))
# Call:
# lm(formula = Y ~ hp + drat + wt, data = mtcars)

# Coefficients:
#             mpg        disp     
# (Intercept)   29.39493   64.52984
# hp            -0.03223    0.66919
# drat           1.61505  -40.10238
# wt            -3.22795   65.97577

llmvn(mvmod, sdmvn_mle(mvmod))
# [1] -238.7386

Я не уверен, что результат правильный или нет.

Кроме того, пожалуйста, дайте мне знать, если есть другие стратегии для вычисления логарифмической вероятности для многомерной линейной регрессии.

...