jackson_mixing_model_doesnotwork.r

library(siar)
################################
#  DEFINE MODEL PARAMETERS
################################
p = c(0.7,0.2,0.1)  # vector of proportions following Jackson

N = 10 # this is the number of data points to generate
num.prey = length(p)
num.iso = 2

sigmaW = 0.5 # following Jackson
sigmaResid = 0 # this can be anything, Jackson used (0,1)

u = matrix(0,num.prey,num.iso)  # This is the matrix of means
u[,1] = c(-10,0,10) # means for isotope 1
u[,2] = c(0, 17.32, 0) # means for isotope 2
sources.sd = matrix(sigmaW,num.prey,num.iso)
sources.dat = cbind(u[,1],sources.sd[,1],u[,2],sources.sd[,2])

sources = (sourcesdemo[1:3,])
sources[,2:5] = sources.dat
nameVec =c("Min","Q1","Median","Mean","Q3","Max","Mode")
outputNames=c(paste(nameVec,".1",sep=""),paste(nameVec,".2",sep=""),paste(nameVec,".3",sep=""),paste(nameVec,".4",sep=""))
numSims = 1
summaryResults = matrix(0, numSims, 28)
colnames(summaryResults)=outputNames

for(ii in 1:numSims) {
  # X is the random data matrix for simulations
  X = matrix(0, N, num.iso)
  Z = matrix(rnorm(N*num.iso,0,1),N, num.iso)
  for(i in 1:N) {
  
    #for(j in 1:num.iso) {
      # loop over isotopes, draw a random source value
    #  X[i,j] = rnorm(num.prey,u[,j],rep(sigmaW,num.prey))%*%p + rnorm(1,0,sigmaResid)
    #}
    for(j in 1:num.prey) {
       X[i,1] = X[i,1] + p[j] * (u[j,1] + sources.sd[j,1]*Z[i,1])
       X[i,2] = X[i,2] + p[j] * (u[j,2] + sources.sd[j,2]*Z[i,2])       
    }
    #rand.sources = matrix(rnorm(num.prey*num.iso,c(u),sigmaW),num.prey,num.iso)
    #sources.dat[,1] = rand.sources[,1]
    #sources.dat[,3] = rand.sources[,2]
    #X[i,] = t(t(rand.sources)%*%p) + rnorm(num.iso,0,sigmaResid)
  }
  # add resid variance down here
  X = X + matrix(rnorm(N*num.iso,0,sigmaResid),N,num.iso)
  model = siarmcmcdirichlet(X, sources.dat)
  # for each parameter calculate the median, posterior mode, and quantiles
  modes = 0
  for(i in 1:4) {
    d = density(model$output[,i])
    modes[i] = d$x[which(d$y==max(d$y))]
  }
  summaryResults[ii,] = c(as.numeric(c(summary(model$output[,1]),modes[1])),as.numeric(c(summary(model$output[,2]),modes[2])),as.numeric(c(summary(model$output[,3]),modes[3])),as.numeric(c(summary(model$output[,4]),modes[4])))
}