NIMBLE for joint model

[nealg...@gmail.com]

Hi,

For marked spatial point process, I need to develop a joint likelihood for mark and spatial point process. How can I write joint likelihood in nimble.

Best

Neal 

[nealg...@gmail.com]

require(nimble)

require(spatstat)

lamb <- function(x,y,a,b) {2*exp(a * x+b*y)}

pp <- rpoispp(lamb, a=3,b=1)

pp_iten=2*exp(3*pp$x+1*pp$y)

z1=runif(length(pp_iten),0,1)

z2=runif(length(pp_iten),0,1)

mark=pp_iten/100*2+z1*0.5+z2*1+rnorm(length(pp_iten),0,1)

code <- nimbleCode({

  beta1~dnorm(0,sd=10)

  beta2~dnorm(0,sd=10)

  alpha1~dnorm(0,sd=10)

  alpha2~dnorm(0,sd=10)

  alpha3~dnorm(0,sd=10)

  lambda0~dgamma(1.5,1.5)

  tau2~dgamma(0.01,0.01)

  for (i in 1:100){

    for (j in 1:100){

      lambda_D[i,j]<-lambda0*exp(beta1*0.01*i+beta2*0.01*j)

    }

  }

  for (i in 1:N){

    lambda[i]<-lambda0*exp(beta1*x[i]+beta2*y[i])

    mark_mean[i]<-alpha1*lambda[i]/100+alpha2*z1[i]+alpha3*z2[i]

    ll[i]<-exp(-tau2*(mark[i]-mark_mean[i])^2/2)*sqrt(tau2)/sqrt(2*3.14)

  }

  s_ll<-mean(lambda_D[1:100,1:100])

  L<-prod(lambda[1:N])*exp(-s_ll)*prod(ll[1:N])

})

constants <- list(N = length(pp_iten))

data <- list(

  x=pp$x,

  y=pp$y,

  mark=mark,

  z1=z1,

  z2=z2

)

inits <- list( beta1 = 0, beta2=0,lambda0=1,alpha1=0,alpha2=0,alpha3=0,tau2=1)

Rmodel <- nimbleModel(code=code, constants=constants, data=data, inits=inits, check = FALSE)

Rmcmc <- buildMCMC(Rmodel)

Cmodel <- compileNimble(Rmodel)

Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

Cmcmc$run(20000)

samples <- as.matrix(Cmcmc$mvSamples[10001:20000,])

summary(as.mcmc(samples))

在 2018年4月12日星期四 UTC-4下午5:48:33，nealg...@gmail.com写道：

[nealg...@gmail.com]

I modify my code by using Customized log-likelihood evaluations:

Have following error:

compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.

Error: Number of dimensions 1 of the return() argument does not match number 0 given in the returnType() statement.

 This occurred for: return(ll)

 This was part of the call:  {

  ll_m <- model_ll_m

  s_ll <- model_s_ll

  log_ll <- model_log_ll

  ll <- ll_m - s_ll + ll_m

  return(ll)

}

My code is :

# data generation

require(nimble)

require(spatstat)

lamb <- function(x,y,a,b) {2*exp(a * x+b*y)}

pp <- rpoispp(lamb, a=3,b=1)

pp_iten=2*exp(3*pp$x+1*pp$y)

z1=runif(length(pp_iten),0,1)

z2=runif(length(pp_iten),0,1)

mark=pp_iten/100*2+z1*0.5+z2*1+rnorm(length(pp_iten),0,1)

code <- nimbleCode({

  beta1~dnorm(0,sd=10)

  beta2~dnorm(0,sd=10)

  alpha1~dnorm(0,sd=10)

  alpha2~dnorm(0,sd=10)

  alpha3~dnorm(0,sd=10)

  lambda0~dgamma(1.5,1.5)

  tau2~dgamma(0.01,0.01)

  for (i in 1:100){

    for (j in 1:100){

      lambda_D[i,j]<-lambda0*exp(beta1*0.01*i+beta2*0.01*j)

    }

  }

  s_ll<-mean(lambda_D[1:100,1:100])

  for (i in 1:N){

    lambda[i]<-lambda0*exp(beta1*x[i]+beta2*y[i])

    mark_mean[i]<-alpha1*lambda[i]/100+alpha2*z1[i]+alpha3*z2[i]

    logmark[i]<--0.5*log(2*3.14)-0.5*log(1/tau2)-tau2/2*pow((mark[i]-mark_mean[i]), 2)

  }

  ll_m<-sum(logmark[1:N])

  log_ll<-log(prod(lambda[1:N]))

})

llFun <- nimbleFunction(

  setup = function(model) {},

  

  run=function(){

    ll_m<-model$ll_m

    s_ll<-model$s_ll

    log_ll<-model$log_ll

    ll<-ll_m-s_ll+ll_m

  returnType(double())

  return(ll)

  }

)

constants <- list(N = length(pp_iten))

data <- list(

  x=pp$x,

  y=pp$y,

  mark=mark,

  z1=z1,

  z2=z2

)

inits <- list( beta1 = 0, beta2=0,lambda0=1,alpha1=0,alpha2=0,alpha3=0,tau2=1)

Rmodel <- nimbleModel(code=code, constants=constants, data=data, inits=inits, check = FALSE)

RllFun<-llFun(Rmodel)

mcmcConf <- configureMCMC(Rmodel,nodes=NULL)

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='tau2',type='RW_llFunction',control=list(llFunction=RllFun,includesTarget=FALSE))

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

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

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

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

Rmcmc <- buildMCMC(mcmcConf)

Cmodel<-compileNimble(Rmodel)

Cmcmc <- compileNimble(Rmcmc, project = Rmodel)

在 2018年4月12日星期四 UTC-4下午5:50:34，nealg...@gmail.com写道：

[Chris Paciorek]

Hi, I'm not too sure what your original question was - in general for
a joint likelihood you would write the likelihoods for each component
using BUGS code (if they factorize in some helpful way) or write a
joint likelihood using a single (multivariate) distribution, either
one that NIMBLE provides or a user-defined distribution.

With regard to your error, it appears that the output of functions
like prod() and sum() are treated as 1-length vectors internally, but
your returnType specifies the output as a scalar. The simple fix is to
do something like this:

return(ll[1])

as that forces the output to be treated as a scalar, which is then
consistent with the returnType being specified as a scalar.

A few other comments;

- In general, prod() is not numerically stable, so you would want to do:
sum(log(lambda[1:N]))

- I'm not fully following the details of what your calculations are
doing without looking more carefully, but note that you may want to
move some of the calculations into your llFun run code rather than
defining the computations in BUGS code. For example, you could
probably move all of the deterministic calculations into the llFun
(see below). That will remove a bunch of nodes from the model graph
and allow NIMBLE to set up the model more quickly and do the
calculations during the MCMC more quickly.

## all this could probably be done in llFun

for (i in 1:100){
for (j in 1:100){
lambda_D[i,j]<-lambda0*exp(beta1*0.01*i+beta2*0.01*j)
}
}
s_ll<-mean(lambda_D[1:100,1:100])
for (i in 1:N){
lambda[i]<-lambda0*exp(beta1*x[i]+beta2*y[i])
mark_mean[i]<-alpha1*lambda[i]/100+alpha2*z1[i]+alpha3*z2[i]
logmark[i]<--0.5*log(2*3.14)-0.5*log(1/tau2)-tau2/2*pow((mark[i]-mark_mean[i]),
2)
}
ll_m<-sum(logmark[1:N])
log_ll<-log(prod(lambda[1:N]))

-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/2f1ad3c1-1af7-4520-80c9-a9623ed7405a%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.

[nimble project]

I think the underlying issue is that model variables are always vectors.  Model scalars are represented as length-1 vectors, but then in a nimbleFunction those are still vectors.

> email to nimble-users+unsubscribe@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/2f1ad3c1-1af7-4520-80c9-a9623ed7405a%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.

--
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+unsubscribe@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/CAG%3DM9Lp6s7HbwOBKvuzRsr_%3DCu30tWjeCcVTyMaUQ0Gd1_Wptg%40mail.gmail.com.