Box Cox transformation

50 views
Skip to first unread message

nealg...@gmail.com

unread,
Apr 3, 2019, 9:53:43 PM4/3/19
to nimble-users
Hi 
I'm doing a project for bayesian inference based box-cox transformation data. I define a box-cox transformation function with nimble function, and give priors on box-cox transformation parameter lambda, and use loglikelihood control the sampling. The code cannot be complied into C++. I attached the code in this email. I think it is something wrong with boxcox transfomation code. Could you help me with this problem?
Best
Guanyu





boxcox<- nimbleFunction( 
  run=function(raw=double(1),lambda=double(1),constant=double(1)){
    returnType(double(1))
    if (lambda == 0) {
      return(log(raw - constant))
    } else {
      return(((raw - constant) ^ lambda - 1) / (lambda))
    }
  })


boxcoxfit<- nimbleCode({
  for (i in 1:n){
    y_star[i]<- boxcox(y[i],lambda,constant) 
    mu[i]<-beta0+beta1*x1[i]+beta2*x2[i]+beta3*x3[i]
    #y_star[i]~dnorm(mu[i],var=sigma2)
    logll[i]<-0.5*log(2*3.14)-0.5*log(sigma2)-1/sigma2*pow((y_star[i]-mu[i]),2)
  }
  ll<-sum(logll[1:n])
  sigma2~dinvgamma(1,1)
  beta0~dnorm(0,sd=10)
  beta1~dnorm(0,sd=10)
  beta2~dnorm(0,sd=10)
  beta3~dnorm(0,sd=10)
  lambda~dnorm(0,sd=10)
})
llFun <- nimbleFunction(
  setup <- function(model) {},
  
  run <- function() {
    ll<-model$ll
    returnType(double())
    return(ll[1])
  }
)
n=100
y=rnorm(100,0,1)
x1=rnorm(100,0,1)
x2=rnorm(100,0,1)
x3=rnorm(100,0,1)


data <- list(
  y =y ,x1=x1,x2=x2,x3=x3
)
constants <- list(
  n = n,constant=-10
)
inits <- list(lambda=0,sigma2=1,beta0=1,beta1=1,beta2=1,beta3=1)

Rmodel <- nimbleModel(boxcoxfit, data = data, constants = constants,inits=inits)

RllFun <- llFun(Rmodel)
mcmcConf <- configureMCMC(Rmodel, nodes = NULL)
mcmcConf$addMonitors(target = c("beta0","beta1", "beta2", "beta3", "sigma2","lambda"))
mcmcConf$addSampler(target = "beta0", type = 'RW_llFunction',
                   control = list(llFunction = RllFun,
                                  includesTarget = FALSE))

mcmcConf$addSampler(target = "beta1", type = 'RW_llFunction',
                    control = list(llFunction = RllFun,
                                   includesTarget = FALSE))

mcmcConf$addSampler(target = "beta2", type = 'RW_llFunction',
                    control = list(llFunction = RllFun,
                                   includesTarget = FALSE))

mcmcConf$addSampler(target = "beta3", type = 'RW_llFunction',
                    control = list(llFunction = RllFun,
                                   includesTarget = FALSE))

mcmcConf$addSampler(target = "sigma2", type = 'RW_llFunction',
                    control = list(llFunction = RllFun,
                                   includesTarget = FALSE))

mcmcConf$addSampler(target = "lambda", type = 'RW_llFunction',
                    control = list(llFunction = RllFun,
                                   includesTarget = FALSE))

Rmcmc <- buildMCMC(mcmcConf)
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
## MCMC run
Cmcmc$run(1000)

Chris Paciorek

unread,
Apr 4, 2019, 10:45:17 AM4/4/19
to Guanyu Hu, nimble-users
Hi Guanyu,

It looks like there are some problems with the dimensions in your
boxcox() function and in your model.
You specify that 'raw', 'lambda', and 'constant' are all vectors in
your nimbleFunction arguments but then
you use 'lambda' as a scalar. Also, in the model, all three of them
are scalars even though they are passed
in to the boxcox() function.

Note that if you run printErrors() after getting a compilation error,
you can often get useful information
that would guide you in determining what the problem is.

As a side note, it might be more straightforward to create a
user-defined boxcox distribution rather than
using the "RW_llFunction" sampler and let NIMBLE choose what sampler
to use (or you could
override NIMBLE's default and choose a different sampler, such as a
slice sampler). But it's up to you.

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 post to this group, send email to nimble...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/8f5cbfec-67f9-4cc0-bb86-a3a1c815e496%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Guanyu Hu

unread,
Apr 8, 2019, 2:46:25 PM4/8/19
to Chris Paciorek, nimble-users
Thanks for careful reply. I will try that.
Best

Guanyu Hu

Perry de Valpine

unread,
Apr 8, 2019, 2:47:57 PM4/8/19
to Guanyu Hu, Chris Paciorek, nimble-users
Normally, this list needs no moderation, but Google's spam filters somehow placed this message on hold for moderator review.  This led to a delay in posting it.

-Perry


Reply all
Reply to author
Forward
0 new messages