Re: KDE user-defined distribution as hmc prior??

[Perry de Valpine]

Hi Alex,

Thanks for the question. In nimble model code (but not nimbleFunction run code), one needs to include square brackets for non-scalars, even if passing the whole object. So you would need 

beta[j] ~ dkde(j,kde_array) to be changed to 

beta[j] ~ dkde(j,kde_array[,,]), if the sizes of kde_array can be easily determined by nimbleModel (which should be the case, and you can use the dimensions argument if needed) or to make it explicit,

beta[j] ~ dkde(j,kde_array[1:n, 1:m, 1:r]), where you would change n, m, and r to whatever is appropriate.

If you haven't had a chance to check out the nimble user manual, it is at r-nimble.org. Also on the documentation page is a link to a short guide to differences in how to code models between jags and nimble, in case you are a jags user.

I am not sure if there aren't other issues. I just went to the most proximate issue with kde_array based on the error message.

HTH and good luck.

Perry

On Fri, Jul 18, 2025 at 11:20 PM Alex <sasha.b...@gmail.com> wrote:

Hi,

This my attempt at a user-defined distribution of a kde prior:

library(nimble)
dkde <- nimbleFunction(
  run = function(x = double(0), j = integer(0), kde_array = double(3), log = integer(0, default = 0)) {
    returnType(double(0))
   
    kde_x <- kde_array[, j, 1]
    kde_y <- kde_array[, j, 2]
   
    #out of bounds
    if (x < min(kde_x) | x > max(kde_x)) return(1e-12)
   
    #find index
    i <- max(which(kde_x <= x))
    if (i == length(kde_x)) return(y[i])
   
    #approxfun
    x0 <- kde_x[i]
    x1 <- kde_x[i+1]
    y0 <- kde_y[i]
    y1 <- kde_y[i+1]
   
    out <- y0 + (y1 - y0)*(x - x0)/(x1 - x0)
   
    if (out <= 0) out <- 1e-12
    if (log) return(log(out))
    else return(out)
  }
)

rkde <- nimbleFunction(
  run = function(n = integer(0), samples = double(1)) {
    returnType(double(0))

    len <- length(samples)
    idx <- as.integer(1 + floor(runif(1, 0, len)))
    out <- samples[idx]

    return(out)
  }
)

registerDistributions(list(
  dkde = list(
    BUGSdist = "dkde(j, kde_array)",
    discrete = FALSE
  )
))

and I wrote this hmc function:

hmc <- function(X, y, niter = 5000, nburnin = 1000, nchains = 1, inits = NULL, KDE = FALSE, C = NA) {
  #0) Define constants
  #---------------
  n = nrow(X)
  M = ncol(X)
 
  #1) Define model
  #---------------
  library(nimble)
  if (KDE == FALSE) {
    log_reg_code <- nimbleCode({
      #Priors
      for (j in 1:M) {
        beta[j] ~ dnorm(mean = 0, sd = 100)
      }
     
      #Likelihood
      for (i in 1:n) {
        logit(p[i]) <- inprod(X[i, 1:M], beta[1:M])
        y[i] ~ dbern(p[i])
      }
    })
  } else {
    n_grid <- 512
    M <- ncol(C)
   
    kde_x <- list()
    kde_y <- list()
    for (j in 1:M) {
      kde <- density(C[,j], n = n_grid)
      kde_x[[j]] <- kde$x
      kde_y[[j]] <- kde$y
    }
    kde_array <- array(0, dim = c(n_grid, M, 2))
    for (j in 1:M) {
      kde_array[, j, 1] <- kde_x[[j]]
      kde_array[, j, 2] <- kde_y[[j]]
    }

    log_reg_code <- nimbleCode({
      #Priors
      for (j in 1:M) {
        beta[j] ~ dkde(j,kde_array)
      }
     
      #Likelihood
      for (i in 1:n) {
        logit(p[i]) <- inprod(X[i, 1:M], beta[1:M])
        y[i] ~ dbern(p[i])
      }
    })
  }
 
  #2) Prepare data, constants, init values
  #---------------
  if (KDE == FALSE) {
    hmc_consts <- list(n = n, M = M)
    hmc_data <- list(y = y, X = X)
    hmc_dim <- list(y = n, X = c(n, M))
  } else {
    hmc_consts <- list(n = n, M = M, kde_array = kde_array)
    hmc_data <- list(y = y, X = X)
    hmc_dim <- list(y = n, X = c(n,M), kde_array = c(n_grid,M,2))
  }
 
  if(is.null(inits)){
    hmc_inits <- list(beta = rep(0, M))
  } else {
    hmc_inits <- list(beta = inits)
  }
 
  #3) Create model
  #---------------
  model <- nimbleModel(
    code = log_reg_code,
    constants = hmc_consts,
    data = hmc_data,
    inits = hmc_inits,
    calculate = FALSE,
    buildDerivs = TRUE,
    dimensions = hmc_dim
  )
 
  #4) Compile model and mcmc to C++ for speed
  #---------------
  compiled_model <- compileNimble(model)
 
  #5) Configure model
  #---------------
  library(nimbleHMC)
  model_config <- configureHMC(model)
  model_config$removeSamplers() #remove all
  model_config$addSampler(target = "beta")
 
  #6) Build HMC algo
  #---------------
  HMC <- buildHMC(model = compiled_model, hmcConf = model_config)
 
  #7) Compile HMC algo
  #---------------
  compiled_HMC <- compileNimble(HMC, project = model)
 
  #8) Run HMC sampler
  #---------------
  samples <- runMCMC(
    compiled_HMC,
    niter = niter,
    nburnin = nburnin,
    nchains = nchains,
    summary = TRUE
  )
 
  return(samples)
}

but when I run: 

hmc(parts_X[[j]], unlist(list(parts_y[[j]])), niter = niter, nburnin = nburnin, nchains = 1, inits = modes, KDE = TRUE, C = C)$samples

I get the following warning:

Warning in replaceConstantsRecurse(x, constEnv, constNames) : Code kde_array was given as known but evaluates to a non-scalar. This is probably not what you want.

and error:

Error in replaceConstantsRecurse(x, constEnv, constNames) : Error, hit end

Essentially, my goal is to create a user-defined kde distribution dkde() created from a previous sample of betas as the new prior for betas in hmc.

Sure, 

      #Priors
      for (j in 1:M) {
        beta[j] ~ dnorm(mean = mean(C[1:n,j]), sd = sd(C[1:n,j]))
      }

works in hmc if I assume the priors are approximately normal, but I don't want to make such assumptions in general, so I'd like to use a custom dkde of the samples C instead. 

Any help would be greatly appreciated. I've been really stuck on this.Thanks!

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/03312d64-705f-4b4a-b36c-59c3dc289609n%40googlegroups.com.