dimension issues with user-defined distribution

[Liz Lawler]

Hi there,

I'm trying to use a custom distribution in a model, and I keep getting the following error:

Error in model$checkBasics() :
  Dimension of distribution argument(s) 'sigma,xi,kappa' does not match required dimension(s) for the distribution 'dEGPDpower'. Necessary dimension(s) are 1,1,1. You may need to ensure that you have explicit vectors and not one-row or one-column matrices.

I've tried a few different ways of doing the matrix multiplication in the likelihood section of my model, but keep getting the same error message. My distribution function does work outside of the model, so I think the issue might be somewhere in my model but I'm not sure where.

Below, I've included my code with some toy data. Thank you so much for your help with this!

Best,

Liz

library(nimble)
library(coda)
library(tidyverse)

g1_cdf_inv <- function(u, sigma = sigma, xi = xi, kappa = kappa) {
  (sigma/xi) * ((1-u^(1/kappa))^-xi - 1)
}
g1_random <- function(n = n, sigma = 1, xi = 0.5, kappa = 1) {
  u <- runif(n)
  return(g1_cdf_inv(u, sigma, xi, kappa))
}

p <- 15
n <- 100
betas_kappa <- rmnorm(1, mean = rep(0,p), varcov = diag(p))
betas_nu <- rmnorm(1, mean = rep(0,p), varcov = diag(p))
betas_xi <- rmnorm(1, mean = rep(0,p), varcov = diag(p))
X <- matrix(rnorm(p*n), nrow = n, ncol = p)
kappa_true <- exp(X %*% t(betas_kappa))
nu_true <- exp(X %*% t(betas_nu))
xi_true <- exp(X %*% t(betas_xi))
sigma_true <- nu_true/(1 + xi_true)
y <- g1_random(n = 100, sigma = sigma_true, xi = xi_true, kappa = kappa_true)

dEGPDpower <- nimbleFunction(
  run = function(x = double(1), sigma = double(1), xi = double(1), kappa = double(1),
                 log = logical(0, default = 1)) {
    returnType(double())
    logProb <- log(kappa) - log(sigma) - (1/xi + 1) * log(1 + xi * (x/sigma)) + (kappa-1) * log(1 - (1 + xi * (x/sigma))^(-1/xi))
    return(logProb)
  }
)


toy_egpd <- nimbleCode({
  # priors
  iden[1:p, 1:p] <- diag(p)
  betask[1:p] ~ dmnorm(mu[1:p], iden[1:p, 1:p])
  betasnu[1:p] ~ dmnorm(mu[1:p], iden[1:p, 1:p])
  betasxi[1:p] ~ dmnorm(mu[1:p], iden[1:p, 1:p])
 
  # likelihood
  for(i in 1:n) {
    log(kappa[i]) <- inprod(X[i, 1:p], betask[1:p])
    log(nu[i]) <- inprod(X[i, 1:p], betasnu[1:p])
    log(xi[i]) <- inprod(X[i, 1:p], betasxi[1:p])
    sigma[i] <- nu[i] / (1 + xi[i])
    y[i] ~ dEGPDpower(sigma[i], xi[i], kappa[i])
  }
})

X <- sweep(X, 2, colMeans(X))
constants <- list(X = X, p = p, n = n, mu = rep(0,p))
data <- list(y = c(y))
inits <- list(betask = rnorm(p), betasnu = rnorm(p), betasxi = rnorm(p))
toy_model <- nimbleModel(toy_egpd, constants = constants, data = data, inits = inits)

[Daniel Turek]

Quickly noting that the type declaration "double(1)" declares a *vector* argumnet, whereas the arguments to your distribution:

dEGPDpower(sigma[i], xi[i], kappa[i])

appear to be scalars.  Also the usage of these arguments in the line beginning

logProb <- ...

appears to treat them as scalars.

To declare these as *scalars* in nimble, you would use either "double()", or equivalently "double(0)".

I haven't tried this change, but wanted to quickly point it out.

--
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 on the web visit https://groups.google.com/d/msgid/nimble-users/b2e7b00a-735e-4294-893e-13e43862b981n%40googlegroups.com.

[Liz Lawler]

I had tried that when I was trying to troubleshoot my issue before posting, but it didn't work. Now that you've mentioned that the scalar piece is likely is the issue, I went in and changed the arguments to "double()" again. It initially didn't work even though I cleared my environment. This did work, however, once I restarted my R session. 

Lesson learned for me to restart R after wiping my environment, as it seems the wrong distribution lingered despite emptying my environment. Thank you for your help!

[Chris Paciorek]

Yeah, we stored user-defined distributions in such a way that to overwrite them, you need to use `deregisterDistributions`. 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/167ce783-ff2f-4131-bccd-26d6cf3650a8n%40googlegroups.com.