Problem using mulitonial

Modeling
1

Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

model {
  vector[N] mu = alpha + beta * x;
  y ~ normal(mu, sigma); 
} 

I must not be understanding how to  properly code the multinomial in Stan. I use Rstan. Here's a simple example. After compiling, when I run the model I get an error message to the effect  that is am attempting to assign reals to a vector. Perhaps someone can identify the error in my code. Here's the example:

//simple_cjs_test.stan
// Jume 11 2018
// Basic cormack-jolly-seber mark-recapture model for one initial marking and release ocassion/site (site 1)
// and two recapture/resight ocassions (sites 2 and 3). Multinomial model:
// Cases: case1: total number marked (M), that survive to site2, are detected (resighted) at site2, 
// and survive and are detected at site 3 (=m11).
// Case2: total number marked that survive to site 2, are not detected at site 2 but do survive
// to site 3 and are detected at site 3 (=m01). Note: survival from site 2 to 3 
// and probability of detection  at site 3 are inseparable, hence one parameter for the joint
// probability of survival from 2 to 3 and detection at 3.
// Case3: total number marked that survive to site 2, are resighted at site 2 and not resigthed thereafter (=m10).
// Case4: total number marked that are never resighted (=m00).
// M = m11 + m01 + m10 + m00; m00 = M - (m11 + m01 + m10.
// Likelihood(dropping the combinatorics: 
// p(m11, m01, m10, m00| s,p,v) = spv^m11 * s(1-p)v^m01 * sp(1-v)^m10 *(s(1-p)(1-v)+(1-s))^m00.

data {
  int<lower=0> n_occasions;     // Mark and resight data 
  int<lower=0> marr[n_occasions]; 
}

parameters {
  vector<lower=0, upper=1>[3] phi; 
  // phi[1]:Survival from inital marking at site1 to first resighting location at site2 (Alternative: s) 
  // phi[2]:Probability of resigthing at site 2 (Alternative: p)
  // phi[3]:Probability of survival from site 2 to site 3 and recapture at site 3. (Alternative: v)
}

transformed parameters {
// simplex[n_occasions] pr[n_occasions];
   vector<lower=0, upper=1>[4] pr;//(Alternative? - I am unsure about this).

  // Define the cell probabilities of pr:
      pr[1] = phi[1] * phi[2];
      pr[2] = phi[1] *(1-phi[2])*phi[3];
  	  pr[3] = phi[1]*phi[2]*(1-phi[3]);
      pr[4] = 1-(pr[1]+pr[2]+pr[3]);
	  
//  Alternative coding:
//  pr[1] = s*p;
//  pr[2] = s*(1-p)*v;
//  pr[3] = s*p*(1-v);
//  pr[4] = 1-(pr[1]+pr[2]+pr[3]);  	
    
}

model {
  // Priors
  // Uniform priors for phi are implicitly defined.

  // Define the multinomial likelihood
  for (t in 1:n_occasions)
    marr[t] ~ multinomial(pr[t]);
}

#simple_cjs_test.R.
#NG June 11 2018

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(123)

## Read data
sdata<- read_rdump("simple_cjs_test_data.R")
n_occasions<-length(sdata$marr)

stan_data = list(sdata,n_occasions)

## Parameters monitored
params <- c("phi")

#Alternative: params<-("s", "p", "v")

## MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 1

## Call Stan from R
simple_test_ex1  <- stan("simple_cjs_test.stan",
                     data = stan_data, pars = params,
                     chains = nc, iter = ni, warmup = nb, thin = nt,
                     seed = 1,
                     open_progress = FALSE)

## Summarize posteriors
print(simple_test_ex1, digits = 3)

marr <-
structure(c(350, 15, 250, 385))

Much thanks for any feedback
2

The random variable for a multinomial is an array of integers with counts of each outcome. If you just have one integer between 1 and K, then you want to use the categorical distribution. See the manual.

3

Bob,
Thanks. But I thought that the Rdump data format is used for the variable marr produces an array of integers. hence my puzzlement regarding the run-time error. How should I code these four integers so that the multinomial sampling statement can read it?

4

In marr[t] ~ multinomial(pr[t]);, the expression marr[t] denotes an integer and pr[t] a real. The multinomial wants an integer array and vector. So the previous sampling statement should probably just be marr ~ multinomial(pr).

5

Bob,

My apologies for not being able to respond sooner. Thanks for the suggestion. I dropped the for t … in the likelihood and recompiled and ran the model. I received an error that marr was not an integer. I eliminated this error by dropping the use of r_dump for the data in the .R file and replaced it with a column vector text file with head ‘marr’,. , then read the data in using read.table. Then it ran fine.

6

Thanks for reporting back!