getSamplesDPMeasure returns replacement length error

[Ethan A]

Hello,

I trying to use the CRP sampler to estimate an ANOVA-DDP type model with a multivariate normal kernel. The sampler runs fine, but when I call getSamplesDPmeasure, I get the following error message:

Error in mvIndexes[l, (tildeVarsColsSum[j] + 1):tildeVarsColsSum[j + 1]] <- which(aux !=  :
  number of items to replace is not a multiple of replacement length

Below is a reproducible example that runs relatively quickly:

library(nimble)
library(posterior)
library(bayesplot)

set.seed(123)

## Nimble model code
code_bnp_mvn_reg = nimbleCode({
  alpha  ~ dgamma(alpha_shape, alpha_rate)
  z[1:n] ~ dCRP(alpha, size = n)
  for ( k in 1:Kmax ) {
    beta[1:pJ, k]      ~ dmnorm(beta_mean[1:pJ], cov = beta_cov[1:pJ, 1:pJ])
    Sigma[1:J, 1:J, k] ~ dinvwish(S = Sigma_scale[1:J, 1:J], df = Sigma_df)
  }
  for ( i in 1:n ) {
    mu[i, 1:J] <- X[id_start[i]:id_end[i], 1:pJ] %*% beta[1:pJ, z[i]]
    y[id_start[i]:id_end[i]] ~ dmnorm(mu[i, 1:J], cov = Sigma[1:J, 1:J, z[i]])
  }
})

##-------------------------------
## Generate simulated data
##-------------------------------


J    = 3       ## Number of time points
n    = 50      ## Sample size
rho  = 0.6     ## Correlation between time points (assumes AR(1) structure)
ptrt = 0.5     ## probability of receiving treatment
pstr = 0.75    ## probability to belong in stratum

## Create AR(1) correlation matrix
Rho = matrix(nrow = J, ncol = J)
for ( s in 1:J )
  for ( J in 1:J )
    Rho[s,J] = rho^(abs(s-J))
sds     = rep(0.15, J)
Sigma   = diag(sds) %*% Rho %*% diag(sds)
epsilon = mvtnorm::rmvnorm(n, sigma = Sigma)
a       = rbinom(n, 1, ptrt)    ## treatment indicator
x       = rbinom(n, 1, pstr)    ## stratum indicator
Mu      = matrix(nrow = n, ncol = J)
for ( i in 1:n ) {
  Mu[i,1] = -0.40 - 0.10 * a[i] + 0.20 * x[i] - 0.10 * a[i] * x[i]
  for ( J in 2:J ) {
    Mu[i,J] = Mu[i,J-1] - 0.05
  }
}
Y = Mu + epsilon

##---------------------------------
## Use MLE to elicit base measure
##---------------------------------
fit.mle       = lm(Y ~ a + x)
summ.mle      = summary(fit.mle)
beta.mle      = coef(fit.mle)
beta.vcov     = vcov(fit.mle)
beta.mean     = as.vector( beta.mle )
beta.cov      = n / 5 * beta.vcov
eta           = model.matrix(fit.mle) %*% beta.mle
Sigma.mle     = crossprod(Y - eta) / n
Sigma.mean    = Sigma.mle
Sigma.df      = J + 4
Sigma.scale   = Sigma.mean * (Sigma.df - J - 1)


##-----------------------------------
## Create variables to easily index
##-----------------------------------
p             = nrow(beta.mle)
beta.start    = integer(J)
beta.end      = integer(J)
beta.start[1] = 1
beta.end[1]   = p
for ( J in 2:J ) {
  beta.end[J]   = beta.end[J-1] + p
  beta.start[J] = beta.end[J] - (p-1)
}
id.start    = integer(n)
id.end      = integer(n)
id.start[1] = 1
id.end[1]   = J
for ( i in 2:n ) {
  id.end[i]   = id.end[i-1] + J
  id.start[i] = id.end[i] - (J-1)
}


## Construct design matrix as block diagonal
Xone  = model.matrix(fit.mle)
X     = matrix(nrow = 0, ncol = p*J)
for ( i in 1:n ) {
  Xi = lapply(1:J, function(t) Xone[i, , drop = FALSE])
  Xi = as.matrix( Matrix::bdiag(Xi) )
  X  = rbind(X, Xi)
}


##-----------------------------------
## Construct nimble model
##-----------------------------------

constants = list(
  Kmax = floor(n/3)      ## maximum number of possible clusters
  , n = n                ## number of observations
  , J = J                ## number of time points
  , pJ = p*J             ## number of regression coefficients
  , id_start = id.start  ## n-dimensional vector giving start index of stacked y's
  , id_end   = id.end    ## n-dimensional vector giving end index of stacked y's
)

data = list(
  y = c(Y), X = X
  , beta_mean = beta.mean, beta_cov = beta.cov
  , Sigma_scale = Sigma.scale, Sigma_df = Sigma.df
  , alpha_shape = 8, alpha_rate = 2
)
inits = list(
  alpha  = 2
  , z     = sample(1:min(constants$Kmax - 1, 5), n, replace = TRUE)  ## initial cluster membership
  , beta  = t(mvtnorm::rmvnorm(constants$Kmax, data$beta_mean, data$beta_cov))
  , Sigma = replicate(constants$Kmax, Sigma.mle)
)

## Build and compile the model and MCMC sampler
model    = nimbleModel(
  code_bnp_mvn_reg, constants, data, inits
)
mcmcConf = configureMCMC(model, monitors = c('alpha', 'beta', 'Sigma', 'z') )
cmodel   = compileNimble(model)
mcmc     = buildMCMC(mcmcConf)
cmcmc    = compileNimble(mcmc, project = cmodel)

## Obtain MCMC samples
fit = runMCMC(cmcmc, niter = 500, nburnin = 0)

## Approximate the Bayesian nonparametric result with a mixture model
fit.stick = getSamplesDPmeasure(cmcmc)

[Chris Paciorek]

Hi Ethan,

Thanks for the report. I'll take a look and see if we can fix this for the next release, which is coming soon-ish (hopefully in the next month or two).

-chris

--
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/38a39b7a-4a90-4509-9995-9e45f488f733n%40googlegroups.com.

[Chris Paciorek]

Hi Ethan,

It looks like the problem is that for a multivariate case like yours, we were expecting `beta` to have as many elements as `mu` (in your case expecting `J` rather than `pJ` elements).

I've made a fix on a GitHub branch. If you'd like to try it out now (or anytime before our next release) you can do:

remotes::install_github("nimble-dev/nimble", ref = "dpmeas-dim-fix", subdir = "packages/nimble")

-chris

[Ethan A]

Hey Chris,

Looks like you pushed it to to the main dev version already. It's working now. Thanks a lot for fixing it so quickly!

Best,

Ethan