yifei-philip
commented
6 days ago I used your method to update splitting prob. and alpha. However, my algorithm always gets stuck in alpha with very little value, which causes only a certain variable to split.
Do you have any ideas about it？
##### auxiliar functions #####
## calculate log(sum(exp(v))), avoiding overflow
log.sum.exp = function(v) {
max.v = max(v)
sum.exp.v.modify = sum(exp(v-max.v))
return (max.v + log(sum.exp.v.modify))
}
## draw rho
draw_rho = function(varcount, p, alpha) {
# rho = rdirichlet(1, as.vector((alpha / p ) + varcount))
rho.par = as.vector((alpha / p ) + varcount)
# Sample unnormalized s on the log scale
templog.rho = rep(NA, p)
for(i in 1:p) {
templog.rho[i] = SoftBart:::rlgam(shape = rho.par[i]) # improvement of drawing from gamma distribution
}
# Normalize s on the log scale, then exponentiate
log.rho = templog.rho - log.sum.exp(templog.rho)
rho = exp(log.rho)
return(rho)
}
## draw alpha
draw_alpha = function(rho, alpha.a0, alpha.b0, p) {
x = 1:1000 / 1001 # generates grids of x = alpha/(alpha+p) ~ beta
alpha.grid = x * p / (1 - x)
alpha.log.lik = lgamma(alpha.grid) - p*lgamma(alpha.grid/p) + alpha.grid * mean(log(rho))
beta.log.prior = (alpha.a0-1)*log(x) + (alpha.b0-1)*log(1-x)
alpha.log.post = alpha.log.lik + beta.log.prior
alpha.log.post.sum = log.sum.exp (alpha.log.post)
# calculate the weights
alpha.weight = exp(alpha.log.post - alpha.log.post.sum)
# sample alpha
alpha = sample(alpha.grid, 1, prob = alpha.weight)
return(alpha)
}
##### DART simulation #####
library(dbarts)
library(SoftBart)
## data generation
f = function(x) {
10 * sin(pi * x[,1] * x[,2]) + 20 * (x[,3] - 0.5)^2 +
10 * x[,4] + 5 * x[,5]
}
set.seed(99)
sigma = 1.0
n = 100
x = matrix(runif(n * 10), n, 10)
Ey = f(x)
y = rnorm(n, Ey, sigma)
data = data.frame(y, x)
## bart sampler initialization
control = dbartsControl(n.samples = 1L,
n.chains = 1L,
n.threads = 1L,keepTrees = TRUE)
rho = rep(1/10, 10)
sampler = dbarts(y ~ ., data = data, test = NULL, control = control,
tree.prior = dbarts:::cgm(power = 2, base = 0.95, split.probs = rho))
## MCMC process
alpha = 1.0
nburn = 1000L
nsave = 1000L
ntotal = nburn + nsave
# storage
sampler.store = list()
sigma.store = array(NA,nsave)
rho.store = array(NA,c(nsave,10))
alpha.store = array(NA,nsave)
for(nrep in 1:ntotal){
# update split prob. after warmup iterations
if(nrep > floor(nburn/2)){
tempmodel = sampler$model
tempmodel@tree.prior@splitProbabilities = rho
sampler$setModel(newModel = tempmodel)
}
# sample from bart/dart
sample = sampler$run()
# update rho and alpha after warmup iterations
if(nrep > floor(nburn/2)){
# draw rho
# tempcounts = table(sampler$getTrees()$var)[-1]
tempcounts = as.vector(sample$varcount)
rho = draw_rho(tempcounts, p = 10, alpha)
# draw alpha
alpha = draw_alpha(rho, alpha.a0 = 0.5, alpha.b0 = 1.0, p = 10)
}
# store samples
if(nrep > nburn){
n = nrep - nburn
sampler.store[[n]] = sampler
sigma.store[n] = sample$sigma
rho.store[n,] = rho
alpha.store[n] = alpha
}
}
EoghanONeill
To implement BART models with hyperpriors on splitting probabilities, for example as in Linero (2018), it is necessary to update the splitting probability in each MCMC iteration. It would be more convenient to do this with a function like
setSplitProbabilities( )for an object of the classdbartsSampler.It is also possible to update splitting probabilities with
setModel( ). I include an example below for anyone who would like to know how to update splitting probabilities. Perhaps an example with the Linero (2018) Dirichlet hyperprior could be added to a vignette.Linero, A. R. (2018). Bayesian regression trees for high-dimensional prediction and variable selection. Journal of the American Statistical Association, 113(522), 626-636.