Как вычислить лог-прогнозирующие баллы в R

Question

Я использую Байесовскую модель усреднения и Байесовскую регрессию Лассо для прогнозирования и хочу оценить точность прогнозов плотности, используя прогностические логарифмические баллы.

Я использую пакет bms для байесовского модельного усреднения и пакет monomvn для байесовского LASSO. В пакете bms функция для вычисления прогнозирующих баллов уже реализована, но для monomvn {bayesian lasso} она не реализована.

Я смог вычислить прогнозируемые плотности из объекта лассо, умножив каждый задний рисунок на соответствующую объясняющую переменную для каждого отдельного наблюдения, поэтому теперь у меня есть прогнозные плотности для каждого наблюдения.

Как я могу оценить прогностический логарифмический балл в R, учитывая прогностическую плотность и реализованные значения?

Лучший

Обновление (решено)

После обращения к одному из авторов пакета BMS, вот моя реализация:

TrainingIdx <- 1:900
TestIdx <- 901:1000

# d = draws
# n = length of test data
# SigmaSq = Error variance draws # has dimension of (1 X d) 
# PredictiveDensity = X.beta draws # has dimension of (n x d)


scores <- matrix(0, nrow = length(TestIdx), ncol = dim(PredictiveDensity)[2]) # create an empty matrix for log-predictive scores

for(obs in 1:length(TestIdx)){ # for each observation
  for(draw in 1:dim(PredictiveDensity)[2]){    # for each draw
    scores[obs,draw] <- dnorm(y[TestIdx,1][obs], mean = PredictiveDensity[obs,draw], sd = sqrt(SigmaSq[draw]))
  }  
}
lps <- -sum(log(rowMeans(scores)))/ length(TestIdx)

Одна из проблем здесь может заключаться в том, что эта реализация не учитывает неопределенность параметров, поэтому вполне вероятно, что она будет предпочтительнее для более крупных моделей, поскольку SigmaSq в выборке имеет тенденцию быть меньше для больших моделей.