How to compute log-predictive scores in R

Question

I am using Bayesian Model Averaging and Bayesian Lasso regression for prediction and I want to evaluate the accuracy of the density forecasts using predictive log-scores.

I am using the bms package for Bayesian Model Averaging and monomvn package for Bayesian LASSO. In the bms package, a function to compute the predictive log-scores is already implemented but for monomvn {bayesian lasso} it is not implemented.

I was able to compute the predictive densities from the lasso object by multiplying each posterior draw with the corresponding explanatory variable for each single observation, so I have the predictive densities for each observation now.

How can I estimate the predictive log-score in R, given the predictive densities and the realized values?

Best

Update (SOLVED)

After contacting one of the authors of BMS package, here is my implementation:

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)

One problem here could be that, this implementation doesn't account for parameter uncertainty, so it is likely that it could favor larger models, as in-sample SigmaSq tends to be smaller for large models.

How to compute log-predictive scores in R

0 Answers0