R/mcmc.conCatSMR.natural.R

Defines functions mcmc.conCatSMR.natural

Documented in mcmc.conCatSMR.natural

#' Fit the categorical spatial mark resight model for natural marks or other situation where
#' the number of marked individuals is unknown
#' @param data a data list as formatted by sim.conCatSMR(). See description for more details.
#' @param niter the number of MCMC iterations to perform
#' @param nburn the number of MCMC iterations to discard as burnin
#' @param nthin the MCMC thinning interval. Keep every nthin iterations.
#' @param M1 the level of data augmentation for marked individuals
#' @param M2 the level of data augmentation for unmarked individuals
#' @param inits a list of initial values for lam0, sigma, gamma, psi1, and psi2. The list element for 
#' gamma is itself a list with ncat elements. See the example below.
#' @param obstype a character string indicating the observation model, "bernoulli" or "poisson".
#' @param nswap an integer indicating how many samples for which the latent identities
#' are updated on each iteration.
#' @param propars a list of proposal distribution tuning parameters for lam0, sigma, s, and st, for the
#' the activity centers of untelemetered and telemetered individuals, respectively. The tuning parameter
#' should be smaller for individuals with telemetry and increasingly so as the number of locations per
#' individual increases.
#' @param storeLatent a logical indicator for whether or not the posteriors of the latent individual identities, z, and s are
#' stored and returned
#' @param storeGamma a logical indicator for whether or not the posteriors for gamma are stored and returned
#' @param IDup a character string indicating whether the latent identity update is done by Gibbs or Metropolis-
#' Hastings, "Gibbs", or "MH". For obstype="bernoulli", only "MH" is available because the full conditional is not known.
#' @description This function fits the conventional categorical spatial mark resight model when the number of
#' marked individuals is unknown.The distribution of marked individuals across the landscape is assumed to be
#' spatially uniform. This can be relaxed by modeling the marking process in the generalized spatial mark resight
#' samplers. An extension of this sampler that allows detection function parameters to vary by the
#' levels of one categorical covariate is located in mcmc.conCatSMR.df(). A version of this sampler
#' for a known number of marked individuals is in mcmc.conCatSMR().
#' 
#' the data list should be formatted to match the list outputted by sim.conCatSMR(), but not all elements
#' of that object are necessary. y.sight.marked, y.sight.unmarked, G.marked, and G.unmarked are necessary
#' list elements. y.sight.x and G.x for x=unk and marke.noID are necessary if there are samples
#' of unknown marked status or samples from marked samples without individual identities.
#' 
#' An element "X", a matrix of trap coordinates, an element "K", the integer number of occasions, and
#' an element n.marked, the integer number of marked individuals are necessary.
#'
#' IDlist is a list containing elements ncat and IDcovs. ncat is an integer for the number
#' of categorical identity covariates and IDcovs is a list of length ncat with elements containing the
#' values each categorical identity covariate may take.
#' 
#' An element "locs", an n.marked x nloc x  2 array of telemetry locations is optional. This array can
#' have missing values if not all individuals have the same number of locations and the entry for individuals
#' with no telemetry should all be missing values (coded NA).
#' 
#' I will write a function to build the data object with "secr-like" input in the near future.
#' 
#' If you need category level detection functions with an unknown number of marks,
#' email Ben. He hasn't written this one, yet, but it will be easy to do.
#' @author Ben Augustine
#' @examples
#' \dontrun{library(coda)
#' N=50
#' n.marked=12
#' lam0=0.35
#' sigma=0.50
#' K=10 #number of occasions
#' buff=3 #state space buffer
#' X<- expand.grid(3:11,3:11) #make a trapping array
#' pMarkID=c(.8,.8)#probability of observing marked status of marked and unmarked individuals
#' pID=.8 #Probability marked individuals are identified
#' ncat=3  #number of ID categories
#' gamma=IDcovs=vector("list",ncat) #population frequencies of each category level. Assume equal here.
#' nlevels=rep(2,ncat) #number of levels per IDcat
#' for(i in 1:ncat){
#'   gamma[[i]]=rep(1/nlevels[i],nlevels[i])
#'   IDcovs[[i]]=1:nlevels[i]
#' }
#' pIDcat=rep(1,ncat)#category observation probabilities
#' tlocs=0 #telemetry locs/marked individual
#' obstype="poisson" #observation model, count or presence/absence?
#' marktype="natural" #premarked or natural ID (marked individuals must be captured)?
#' data=sim.conCatSMR(N=N,n.marked=n.marked,lam0=lam0,sigma=sigma,K=K,X=X,buff=buff,obstype=obstype,ncat=ncat,
#'                       pIDcat=pIDcat,gamma=gamma,IDcovs=IDcovs,pMarkID=pMarkID,tlocs=tlocs,pID=pID,marktype=marktype)
#' #MCMC
#' inits=list(lam0=lam0,sigma=sigma,gamma=gamma,psi1=0.5,psi2=0.5) #start at simulated values
#' proppars=list(lam0=0.1,sigma=0.08,s=0.5,st=0.08) #st only for telemetered inds. Should be smaller than s.
#' M1=50 #marked data augmentation
#' M2=100 #unmarked data augmentation
#' storeLatent=TRUE
#' storeGamma=FALSE
#' niter=1000
#' nburn=0
#' nthin=1
#' IDup="Gibbs"
#' out=mcmc.conCatSMR.natural(data,niter=niter,nburn=nburn, nthin=nthin, M1 = M1,M2=M2, inits=inits,obstype=obstype,
#'                               proppars=proppars,storeLatent=TRUE,storeGamma=TRUE,IDup=IDup)
#' 
#' 
#' plot(mcmc(out$out))
#' 
#' 1-rejectionRate(mcmc(out$out)) #shoot for 0.2 - 0.4 for lam0 and sigma. If too low, raise proppar. If too high, lower proppar.
#' 1-rejectionRate(mcmc(out$s1xout)) #marked activity center acceptance in x dimension. Shoot for min of 0.2
#' 1-rejectionRate(mcmc(out$s2xout)) #unmarked activity center acceptance in x dimension. Shoot for min of 0.2
#' 
#' #Example for regular conventional SMR (no identity covariates)
#' N=50
#' n.marked=12
#' lam0=0.35
#' sigma=0.50
#' K=10 #number of occasions
#' buff=3 #state space buffer
#' X<- expand.grid(3:11,3:11) #make a trapping array
#' pMarkID=c(.8,.8)#probability of observing marked status of marked and unmarked individuals
#' pID=.8 #Probability marked individuals are identified
#' ncat=1  #number of ID categories
#' gamma=IDcovs=vector("list",ncat) #population frequencies of each category level. Assume equal here.
#' nlevels=rep(1,ncat) #number of levels per IDcat
#' for(i in 1:ncat){
#'   gamma[[i]]=rep(1/nlevels[i],nlevels[i])
#'   IDcovs[[i]]=1:nlevels[i]
#' }
#' pIDcat=rep(1,ncat)#category observation probabilities
#' tlocs=0 #telemetry locs/marked individual
#' obstype="poisson" #observation model, count or presence/absence?
#' marktype="natural" #premarked or natural ID (marked individuals must be captured)?
#' data=sim.conCatSMR(N=N,n.marked=n.marked,lam0=lam0,sigma=sigma,K=K,X=X,buff=buff,obstype=obstype,ncat=ncat,
#'                       pIDcat=pIDcat,gamma=gamma,IDcovs=IDcovs,pMarkID=pMarkID,tlocs=tlocs,pID=pID,marktype=marktype)
#' #MCMC
#' inits=list(lam0=lam0,sigma=sigma,gamma=gamma,psi1=0.5,psi2=0.5) #start at simulated values
#' proppars=list(lam0=0.1,sigma=0.08,s=0.5,st=0.08) #st only for telemetered inds. Should be smaller than s.
#' M1=50 #marked data augmentation
#' M2=100 #unmarked data augmentation
#' storeLatent=TRUE
#' storeGamma=FALSE
#' niter=1000
#' nburn=0
#' nthin=1
#' IDup="Gibbs"
#' out=mcmc.conCatSMR.natural(data,niter=niter,nburn=nburn, nthin=nthin, M1 = M1,M2=M2, inits=inits,obstype=obstype,
#'                               proppars=proppars,storeLatent=TRUE,storeGamma=TRUE,IDup=IDup)
#' 
#' 
#' plot(mcmc(out$out))
#' }
#' @export
mcmc.conCatSMR.natural <-
  function(data,niter=2400,nburn=1200, nthin=5, M1 = 30,M2=200, inits=NA,obstype="poisson",nswap=NA,
           proppars=list(lam0=0.05,sigma=0.1,sx=0.2,sy=0.2),
           storeLatent=TRUE,storeGamma=TRUE,IDup="Gibbs",tf=NA,priors=NA){
    ###
    library(abind)
    y.sight.marked=data$y.sight.marked
    y.sight.unmarked=data$y.sight.unmarked
    X<-as.matrix(data$X)
    J<-nrow(X)
    K<- dim(y.sight.marked)[3]
    ncat=data$IDlist$ncat
    IDcovs=data$IDlist$IDcovs
    buff<- data$buff
    marked.guys=1:dim(y.sight.marked)[1]
    G.marked=data$G.marked
    G.unmarked=data$G.unmarked
  
    if(!is.matrix(G.marked)){
      G.marked=matrix(G.marked)
    }
    if(!is.matrix(G.unmarked)){
      G.unmarked=matrix(G.unmarked)
    }
    if(!is.list(IDcovs)){
      stop("IDcovs must be a list")
    }
    nIDcovs=unlist(lapply(IDcovs,length))
    if(ncol(G.marked)!=ncat){
      stop("G.marked needs ncat number of columns")
    }
    if(ncol(G.unmarked)!=ncat){
      stop("G.unmarked needs ncat number of columns")
    }
    if(!all(is.na(priors))){
      if(!is.list(priors))stop("priors must be a list" )
      if(!all(names(priors)==c("sigma")))stop("priors list element name must be sigma")
      warning("Using Gamma prior for sigma")
      usePriors=TRUE
    }else{
      warning("No sigma prior entered, using uniform(0,infty).")
      usePriors=FALSE
    }
    #Are there unknown marked status guys?
    useUnk=FALSE
    if("G.unk"%in%names(data)){
      if(!is.na(data$G.unk[1])){
        G.unk=data$G.unk
        if(!is.matrix(G.unk)){
          G.marked=matrix(G.unk)
        }
        if(ncol(G.unk)!=ncat){
          stop("G.unk needs ncat number of columns")
        }
        y.sight.unk=data$y.sight.unk
        useUnk=TRUE
      }
    }
    #Are there marked no ID guys?
    useMarkednoID=FALSE
    if("G.marked.noID"%in%names(data)){
      if(!is.na(data$G.marked.noID[1])){
        G.marked.noID=data$G.marked.noID
        if(!is.matrix(G.marked.noID)){
          G.marked.noID=matrix(G.marked.noID)
        }
        if(ncol(G.marked.noID)!=ncat){
          stop("G.marked.noID needs ncat number of columns")
        }
        y.sight.marked.noID=data$y.sight.marked.noID
        useMarkednoID=TRUE
      }
    }

    #data checks
    if(length(dim(y.sight.marked))!=3){
      stop("dim(y.sight.marked) must be 3. Reduced to 2 during initialization")
    }
    if(length(dim(y.sight.unmarked))!=3){
      stop("dim(y.sight.unmarked) must be 3. Reduced to 2 during initialization")
    }
    if(useUnk){
      if(length(dim(y.sight.unk))!=3){
        stop("dim(y.sight.unk) must be 3. Reduced to 2 during initialization")
      }
    }
    if(useMarkednoID){
      if(length(dim(y.sight.marked.noID))!=3){
        stop("dim(y.sight.marked.noID) must be 3. Reduced to 2 during initialization")
      }
    }

    if(!IDup%in%c("MH","Gibbs")){
      stop("IDup must be MH or Gibbs")
    }
    if(IDup=="MH"){
      # stop("MH not implemented, yet")
    }
    if(obstype=="bernoulli"&IDup=="Gibbs"){
      stop("Must use MH IDup for bernoulli data")
    }
    
    #If using polygon state space
    if("vertices"%in%names(data)){
      vertices=data$vertices
      useverts=TRUE
      xlim=c(min(vertices[,1]),max(vertices[,1]))
      ylim=c(min(vertices[,2]),max(vertices[,2]))
    }else if("buff"%in%names(data)){
      buff<- data$buff
      xlim<- c(min(X[,1]),max(X[,1]))+c(-buff, buff)
      ylim<- c(min(X[,2]),max(X[,2]))+c(-buff, buff)
      vertices=cbind(xlim,ylim)
      useverts=FALSE
    }else{
      stop("user must supply either 'buff' or 'vertices' in data object")
    }
    if(!any(is.na(tf))){
      if(any(tf>K)){
        stop("Some entries in tf are greater than K.")
      }
      if(is.null(dim(tf))){
        if(length(tf)!=J){
          stop("tf vector must be of length J.")
        }
        K2D.M=matrix(rep(tf,M1),nrow=M1,ncol=J,byrow=TRUE)
        K2D.UM=matrix(rep(tf,M2),nrow=M2,ncol=J,byrow=TRUE)
      }else{
        stop("2D tf not accepted by this function since marks are natural so unlikely deaths are known")
        # if(!all(dim(tf)==c(M,J))){
        #   stop("tf must be dim M by J if tf varies by individual")
        # }
        # K2D=tf
        # warning("Since 2D tf entered, assuming individual exposure to traps differ")
      }
    }else{
      tf=rep(K,J)
      K2D.M=matrix(rep(tf,M1),nrow=M1,ncol=J,byrow=TRUE)
      K2D.UM=matrix(rep(tf,M2),nrow=M2,ncol=J,byrow=TRUE)
    }
    ##pull out initial values
    if(!(("psi1"%in%names(inits))|("psi2"%in%names(inits)))){
      stop("Must supply psi1 and psi2 inits.")
    }
    psi1<- inits$psi1
    psi2<- inits$psi2
    lam0=inits$lam0
    sigma<- inits$sigma
    gamma=inits$gamma
    if(!is.list(gamma)){
      stop("inits$gamma must be a list")
    }
    
    if(useUnk&!useMarkednoID){
      G.use=rbind(G.unmarked,G.unk)
      status=c(rep(2,nrow(G.unmarked)),rep(0,nrow(G.unk)))
      G.use=cbind(G.use,status)
      G.marked=cbind(G.marked,rep(1,nrow(G.marked)))
      ncat=ncat+1
      y.sight.latent=abind(y.sight.unmarked,y.sight.unk,along=1)
    }else if(!useUnk&useMarkednoID){
      G.use=rbind(G.unmarked,G.marked.noID)
      status=c(rep(2,nrow(G.unmarked)),rep(1,nrow(G.marked.noID)))
      G.use=cbind(G.use,status)
      G.marked=cbind(G.marked,rep(1,nrow(G.marked)))
      ncat=ncat+1
      y.sight.latent=abind(y.sight.unmarked,y.sight.marked.noID,along=1)
    }else if(useUnk&useMarkednoID){
      G.use=rbind(G.unmarked,G.unk,G.marked.noID)
      status=c(rep(2,nrow(G.unmarked)),rep(0,nrow(G.unk)),rep(1,nrow(G.marked.noID)))
      G.use=cbind(G.use,status)
      G.marked=cbind(G.marked,rep(1,nrow(G.marked)))
      ncat=ncat+1
      nIDcovs=c(nIDcovs,2)
      y.sight.latent=abind(y.sight.unmarked,y.sight.unk,y.sight.marked.noID,along=1)
    }else{
      G.use=G.unmarked
      y.sight.latent=y.sight.unmarked
    }
    n.samp.latent=nrow(y.sight.latent)
    if(is.na(nswap)){
      nswap=n.samp.latent/2
      warning("nswap not specified, using n.samp.latent/2")
    }
    
    #make constraints for data initialization
    constraints=matrix(1,nrow=n.samp.latent,ncol=n.samp.latent)
    for(i in 1:n.samp.latent){
      for(j in 1:n.samp.latent){
        guys1=which(G.use[i,]!=0)
        guys2=which(G.use[j,]!=0)
        comp=guys1[which(guys1%in%guys2)]
        if(any(G.use[i,comp]!=G.use[j,comp])){
          constraints[i,j]=0
        }
      }
    }
    #If bernoulli data, add constraints that prevent y.true[i,j,k]>1
    binconstraints=FALSE
    if(obstype=="bernoulli"){
      idx=t(apply(y.sight.latent,1,function(x){which(x>0,arr.ind=TRUE)}))
      for(i in 1:n.samp.latent){
        for(j in 1:n.samp.latent){
          if(i!=j){
            if(all(idx[i,1:2]==idx[j,1:2])){
              constraints[i,j]=0 #can't combine samples from same trap and occasion in binomial model
              constraints[j,i]=0
              binconstraints=TRUE
            }
          }
        }
      }
    }
    
    #Build y.sight.true
    y.sight.marked.true=array(0,dim=c(M1,J,K))
    y.sight.unmarked.true=array(0,dim=c(M2,J,K))
    y.sight.marked.true[marked.guys,,]=y.sight.marked
    ID=matrix(NA,nrow=n.samp.latent,ncol=2) #second column: 1 indicates assigned to marked, 2 to unmarked
    idx=1
    for(i in 1:n.samp.latent){
      if(useMarkednoID){
        if(status[i]==1)next
      }
      if(idx>M2){
        stop("Need to raise M2 to initialize y.true")
      }
      traps=which(rowSums(y.sight.latent[i,,])>0)
      y.sight.unmarked.true2D=apply(y.sight.unmarked.true,c(1,2),sum)
      if(length(traps)==1){
        cand=which(y.sight.unmarked.true2D[,traps]>0)#guys caught at same traps
      }else{
        cand=which(rowSums(y.sight.unmarked.true2D[,traps])>0)#guys caught at same traps
      }
      if(length(cand)>0){
        if(length(cand)>1){#if more than 1 ID to match to, choose first one
          cand=cand[1]
        }
        #Check constraint matrix
        cands=which(ID[,1]%in%cand)#everyone assigned this ID
        if(all(constraints[i,cands]==1)){#focal consistent with all partials already assigned
          y.sight.unmarked.true[cand,,]=y.sight.unmarked.true[cand,,]+y.sight.latent[i,,]
          ID[i,]=c(cand,2)
        }else{#focal not consistent
          y.sight.unmarked.true[idx,,]=y.sight.latent[i,,]
          ID[i,]=c(idx,2)
          idx=idx+1
        }
      }else{#no assigned samples at this trap
        y.sight.unmarked.true[idx,,]=y.sight.latent[i,,]
        ID[i,]=c(idx,2)
        idx=idx+1
      }
    }
    if(useMarkednoID){#Need to initialize these guys to marked guys
      fix=which(status==1)
      meanloc=matrix(NA,nrow=n.marked,ncol=2)
      for(i in marked.guys){
        trap=which(rowSums(y.sight.marked[i,,])>0)
        locs2=matrix(0,nrow=0,ncol=2)
        if(length(trap)>0){
          locs2=rbind(locs2,X[trap,])
        }
        if(nrow(locs2)>1){
          meanloc[i,]=colMeans(locs2)
        }else if(nrow(locs2)>0){
          meanloc[i,]=locs2
        }
      }
      for(i in fix){
        trap=which(rowSums(y.sight.latent[i,,])>0)
        dists=sqrt((X[trap,1]-meanloc[,1])^2+(X[trap,2]-meanloc[,2])^2)
        ID[i,]=c(which(dists==min(dists,na.rm=TRUE))[1],1)
        y.sight.marked.true[ID[i],,]=y.sight.marked.true[ID[i],,]+y.sight.latent[i,,]
      }
    }
    
    #Check assignment consistency with constraints
    checkID=unique(ID)
    checkID=checkID[!checkID%in%marked.guys]
    for(i in 1:length(checkID)){
      idx=which(ID==checkID[i])
      if(!all(constraints[idx,idx]==1)){
        stop("ID initialized improperly")
      }
    }
   
    y.sight.marked.true=apply(y.sight.marked.true,c(1,2),sum)
    y.sight.unmarked.true=apply(y.sight.unmarked.true,c(1,2),sum)
    y.sight.latent=apply(y.sight.latent,c(1,2),sum)
    
    known.vector.marked=1*(rowSums(y.sight.marked.true)>0)
    known.vector.unmarked=1*(rowSums(y.sight.unmarked.true)>0)
    
    #Initialize zs
    z1=1*(known.vector.marked>0)
    z2=1*(known.vector.unmarked>0)
    add1=M1*(0.5-sum(z1)/M1)
    if(add1>0){
      z1[sample(which(z1==0),add1)]=1 #switch some uncaptured z's to 1.
    }
    add2=M2*(0.5-sum(z2)/M2)
    if(add2>0){
      z2[sample(which(z2==0),add2)]=1 #switch some uncaptured z's to 1.
    }
    
    #Optimize starting locations given where they are trapped.
    s1<- cbind(runif(M1,xlim[1],xlim[2]), runif(M1,ylim[1],ylim[2])) #assign random locations
    idx=which(rowSums(y.sight.marked.true)>0) #switch for those actually caught
    for(i in idx){
      trps<- matrix(X[y.sight.marked.true[i,]>0,1:2],ncol=2,byrow=FALSE)
      if(nrow(trps)>1){
        s1[i,]<- c(mean(trps[,1]),mean(trps[,2]))
      }else{
        s1[i,]<- trps
      }
    }
    s2<- cbind(runif(M2,xlim[1],xlim[2]), runif(M2,ylim[1],ylim[2])) #assign random locations
    idx=which(rowSums(y.sight.unmarked.true)>0) #switch for those actually caught
    for(i in idx){
      trps<- matrix(X[y.sight.unmarked.true[i,]>0,1:2],ncol=2,byrow=FALSE)
      if(nrow(trps)>1){
        s2[i,]<- c(mean(trps[,1]),mean(trps[,2]))
      }else{
        s2[i,]<- trps
      }
    }
    if(useverts==TRUE){
      inside=rep(NA,nrow(s1))
      for(i in 1:nrow(s1)){
        inside[i]=inout(s1[i,],vertices)
      }
      idx=which(inside==FALSE)
      if(length(idx)>0){
        for(i in 1:length(idx)){
          while(inside[idx[i]]==FALSE){
            s1[idx[i],]=c(runif(1,xlim[1],xlim[2]), runif(1,ylim[1],ylim[2]))
            inside[idx[i]]=inout(s1[idx[i],],vertices)
          }
        }
      }
      inside=rep(NA,nrow(s2))
      for(i in 1:nrow(s2)){
        inside[i]=inout(s2[i,],vertices)
      }
      idx=which(inside==FALSE)
      if(length(idx)>0){
        for(i in 1:length(idx)){
          while(inside[idx[i]]==FALSE){
            s2[idx[i],]=c(runif(1,xlim[1],xlim[2]), runif(1,ylim[1],ylim[2]))
            inside[idx[i]]=inout(s2[idx[i],],vertices)
          }
        }
      }
    }
    #Initialize G.true
    G.true.marked=matrix(0,nrow=M1,ncol=ncat)
    G.true.marked[marked.guys,]=G.marked
    G.true.unmarked=matrix(0,nrow=M2,ncol=ncat)
    for(i in unique(ID[,1])){
      idx=which(ID[,1]==i&ID[,2]==2)
      if(length(idx)==1){
        G.true.unmarked[i,]=G.use[idx,]
      }else{
        if(ncol(G.use)>1){
          G.true.unmarked[i,]=apply(G.use[idx,],2, max) #consensus
        }else{
          G.true.unmarked[i,]=max(G.use[idx,])
        }
      }
    }
    if(useUnk|useMarkednoID){#augmented guys are unmarked.
      if(max(ID)<M2){
        G.true.unmarked[,ncol(G.true.unmarked)]=2
      }
      unkguys=which(G.use[,ncol(G.use)]==0)
    }
    
    G.latent.marked=G.true.marked==0#Which genos can be updated?
    G.latent.unmarked=G.true.unmarked==0#Which genos can be updated?
    if(!(useUnk|useMarkednoID)){
      for(j in 1:(ncat)){
        fix=G.true.marked[,j]==0
        G.true.marked[fix,j]=sample(IDcovs[[j]],sum(fix),replace=TRUE,prob=gamma[[j]])
        fix=G.true.unmarked[,j]==0
        G.true.unmarked[fix,j]=sample(IDcovs[[j]],sum(fix),replace=TRUE,prob=gamma[[j]])
      }
    }else{
      for(j in 1:(ncat-1)){
        fix=G.true.marked[,j]==0
        G.true.marked[fix,j]=sample(IDcovs[[j]],sum(fix),replace=TRUE,prob=gamma[[j]])
        fix=G.true.unmarked[,j]==0
        G.true.unmarked[fix,j]=sample(IDcovs[[j]],sum(fix),replace=TRUE,prob=gamma[[j]])
      }
      #Split marked status back off
      Mark.obs=G.use[,ncat]
      # Mark.status=G.true[,ncat]
      ncat=ncat-1
      G.use=G.use[,1:ncat]
      G.true.marked=G.true.marked[,1:ncat]
      G.true.unmarked=G.true.unmarked[,1:ncat]
    }
    if(!is.matrix(G.use)){
      G.use=matrix(G.use,ncol=1)
    }
    if(!is.matrix(G.true.marked)){
      G.true.marked=matrix(G.true.marked,ncol=1)
    }
    if(!is.matrix(G.true.unmarked)){
      G.true.unmarked=matrix(G.true.unmarked,ncol=1)
    }
    # some objects to hold the MCMC output
    nstore=(niter-nburn)/nthin
    if(nburn%%nthin!=0){
      nstore=nstore+1
    }
    out<-matrix(NA,nrow=nstore,ncol=10)
    dimnames(out)<-list(NULL,c("lam0","sigma","n.m","n.um","n.tot","N.m","N.um","N.all","psi1","psi2"))
    if(storeLatent){
      s1xout<- s1yout<- z1out<-matrix(NA,nrow=nstore,ncol=M1)
      s2xout<- s2yout<- z2out<-matrix(NA,nrow=nstore,ncol=M2)
      IDout=array(NA,dim=c(nstore,nrow(ID),2))
    }
    iteridx=1 #for storing output not recorded every iteration
    if(storeGamma){
      gammaOut=vector("list",ncat)
      for(i in 1:ncat){
        gammaOut[[i]]=matrix(NA,nrow=nstore,ncol=nIDcovs[i])
        colnames(gammaOut[[i]])=paste("Lo",i,"G",1:nIDcovs[i],sep="")
      }
    }
    if(!is.na(data$locs[1])){
      uselocs=TRUE
      locs=data$locs
      telguys=which(rowSums(!is.na(locs[,,1]))>0)
      #update starting locations using telemetry data
      for(i in telguys){
        s1[i,]<- c(mean(locs[i,,1]),mean(locs[i,,2]))
      }
      ll.tel=matrix(0,nrow=length(telguys),ncol=dim(locs)[2])
      for(i in telguys){
        ll.tel[i,]=dnorm(locs[i,,1],s1[i,1],sigma,log=TRUE)+dnorm(locs[i,,2],s1[i,2],sigma,log=TRUE)
      }
      ll.tel.cand=ll.tel
    }else{
      uselocs=FALSE
      telguys=c()
    }
    
    D1=e2dist(s1, X)
    D2=e2dist(s2, X)
    lamd.sight.marked<- lam0*exp(-D1*D1/(2*sigma*sigma))
    lamd.sight.unmarked<- lam0*exp(-D2*D2/(2*sigma*sigma))
    ll.y.sight.marked=array(0,dim=c(M1,J))
    ll.y.sight.unmarked=array(0,dim=c(M2,J))
    if(obstype=="bernoulli"){
      pd.sight.marked=1-exp(-lamd.sight.marked)
      pd.sight.unmarked=1-exp(-lamd.sight.unmarked)
      pd.sight.marked.cand=pd.sight.marked
      pd.sight.unmarked.cand=pd.sight.unmarked
      ll.y.sight.marked=dbinom(y.sight.marked.true,K2D.M,pd.sight.marked*z1,log=TRUE)
      ll.y.sight.unmarked=dbinom(y.sight.unmarked.true,K2D.UM,pd.sight.unmarked*z2,log=TRUE)
    }else if(obstype=="poisson"){
      ll.y.sight.marked=dpois(y.sight.marked.true,K2D.M*lamd.sight.marked*z1,log=TRUE)
      ll.y.sight.unmarked=dpois(y.sight.unmarked.true,K2D.UM*lamd.sight.unmarked*z2,log=TRUE)
    }
    lamd.sight.marked.cand=lamd.sight.marked
    lamd.sight.unmarked.cand=lamd.sight.unmarked
    ll.y.sight.marked.cand=ll.y.sight.marked
    ll.y.sight.unmarked.cand=ll.y.sight.unmarked
    y.sight.marked.cand=y.sight.marked.true
    y.sight.unmarked.cand=y.sight.unmarked.true
    if(!is.finite(sum(ll.y.sight.marked)))stop("Starting likelihood not finite. 
                                  Try raising lam0 and/or sigma inits.")
    if(!is.finite(sum(ll.y.sight.unmarked)))stop("Starting likelihood not finite. 
                                  Try raising lam0 and/or sigma inits.")
    for(iter in 1:niter){
      #Update both observation models
      if(obstype=="bernoulli"){
        #Update lam0
        llysightsum=sum(ll.y.sight.marked)+sum(ll.y.sight.unmarked)
        lam0.cand<- rnorm(1,lam0,proppars$lam0)
        if(lam0.cand > 0){
          lamd.sight.marked.cand<- lam0.cand*exp(-D1*D1/(2*sigma*sigma))
          lamd.sight.unmarked.cand<- lam0.cand*exp(-D2*D2/(2*sigma*sigma))
          pd.sight.marked.cand=1-exp(-lamd.sight.marked.cand)
          pd.sight.unmarked.cand=1-exp(-lamd.sight.unmarked.cand)
          ll.y.sight.marked.cand= dbinom(y.sight.marked.true,K2D.M,pd.sight.marked.cand*z1,log=TRUE)
          ll.y.sight.unmarked.cand= dbinom(y.sight.unmarked.true,K2D.UM,pd.sight.unmarked.cand*z2,log=TRUE)
          llysightcandsum=sum(ll.y.sight.marked.cand)+sum(ll.y.sight.unmarked.cand)
          if(runif(1) < exp(llysightcandsum-llysightsum)){
            lam0<- lam0.cand
            lamd.sight.marked=lamd.sight.marked.cand
            lamd.sight.unmarked=lamd.sight.unmarked.cand
            pd.sight.marked=pd.sight.marked.cand
            pd.sight.unmarked=pd.sight.unmarked.cand
            ll.y.sight.marked=ll.y.sight.marked.cand
            ll.y.sight.unmarked=ll.y.sight.unmarked.cand
            llysightsum=llysightcandsum
          }
        }
        #Update sigma
        sigma.cand<- rnorm(1,sigma,proppars$sigma)
        if(sigma.cand > 0){
          lamd.sight.marked.cand<- lam0*exp(-D1*D1/(2*sigma.cand*sigma.cand))
          lamd.sight.unmarked.cand<- lam0*exp(-D2*D2/(2*sigma.cand*sigma.cand))
          pd.sight.marked.cand=1-exp(-lamd.sight.marked.cand)
          pd.sight.unmarked.cand=1-exp(-lamd.sight.unmarked.cand)
          ll.y.sight.marked.cand= dbinom(y.sight.marked.true,K2D.M,pd.sight.marked.cand*z1,log=TRUE)
          ll.y.sight.unmarked.cand= dbinom(y.sight.unmarked.true,K2D.UM,pd.sight.unmarked.cand*z2,log=TRUE)
          llysightcandsum=sum(ll.y.sight.marked.cand)+sum(ll.y.sight.unmarked.cand)
          if(uselocs){
            for(i in telguys){
              ll.tel.cand[i,]=dnorm(locs[i,,1],s1[i,1],sigma.cand,log=TRUE)+dnorm(locs[i,,2],s1[i,2],sigma.cand,log=TRUE)
            }
          }else{
            ll.tel.cand=ll.tel=0
          }
          if(usePriors){
            prior.curr=dgamma(sigma,priors$sigma[1],priors$sigma[2],log=TRUE)
            prior.cand=dgamma(sigma.cand,priors$sigma[1],priors$sigma[2],log=TRUE)
          }else{
            prior.curr=prior.cand=0
          }
          if(runif(1) < exp((llysightcandsum+sum(ll.tel.cand)+prior.cand)-(llysightsum+sum(ll.tel)+prior.curr))){
            sigma<- sigma.cand
            lamd.sight.marked=lamd.sight.marked.cand
            lamd.sight.unmarked=lamd.sight.unmarked.cand
            pd.sight.marked=pd.sight.marked.cand
            pd.sight.unmarked=pd.sight.unmarked.cand
            ll.y.sight.marked=ll.y.sight.marked.cand
            ll.y.sight.unmarked=ll.y.sight.unmarked.cand
            ll.tel=ll.tel.cand
          }
        }
      }else{#poisson
        #Update lam0
        llysightsum=sum(ll.y.sight.marked)+sum(ll.y.sight.unmarked)
        lam0.cand<- rnorm(1,lam0,proppars$lam0)
        if(lam0.cand > 0){
          lamd.sight.marked.cand<- lam0.cand*exp(-D1*D1/(2*sigma*sigma))
          lamd.sight.unmarked.cand<- lam0.cand*exp(-D2*D2/(2*sigma*sigma))
          ll.y.sight.marked.cand= dpois(y.sight.marked.true,K2D.M*lamd.sight.marked.cand*z1,log=TRUE)
          ll.y.sight.unmarked.cand= dpois(y.sight.unmarked.true,K2D.UM*lamd.sight.unmarked.cand*z2,log=TRUE)
          llysightcandsum=sum(ll.y.sight.marked.cand)+sum(ll.y.sight.unmarked.cand)
          if(runif(1) < exp(llysightcandsum-llysightsum)){
            lam0<- lam0.cand
            lamd.sight.marked=lamd.sight.marked.cand
            lamd.sight.unmarked=lamd.sight.unmarked.cand
            ll.y.sight.marked=ll.y.sight.marked.cand
            ll.y.sight.unmarked=ll.y.sight.unmarked.cand
            llysightsum=llysightcandsum
          }
        }
        # Update sigma
        sigma.cand<- rnorm(1,sigma,proppars$sigma)
        if(sigma.cand > 0){
          lamd.sight.marked.cand<- lam0*exp(-D1*D1/(2*sigma.cand*sigma.cand))
          lamd.sight.unmarked.cand<- lam0*exp(-D2*D2/(2*sigma.cand*sigma.cand))
          ll.y.sight.marked.cand= dpois(y.sight.marked.true,K2D.M*lamd.sight.marked.cand*z1,log=TRUE)
          ll.y.sight.unmarked.cand= dpois(y.sight.unmarked.true,K2D.UM*lamd.sight.unmarked.cand*z2,log=TRUE)
          llysightcandsum=sum(ll.y.sight.marked.cand)+sum(ll.y.sight.unmarked.cand)
          if(uselocs){
            for(i in telguys){
              ll.tel.cand[i,]=dnorm(locs[i,,1],s1[i,1],sigma.cand,log=TRUE)+dnorm(locs[i,,2],s1[i,2],sigma.cand,log=TRUE)
            }
          }else{
            ll.tel.cand=ll.tel=0
          }
          if(usePriors){
            prior.curr=dgamma(sigma,priors$sigma[1],priors$sigma[2],log=TRUE)
            prior.cand=dgamma(sigma.cand,priors$sigma[1],priors$sigma[2],log=TRUE)
          }else{
            prior.curr=prior.cand=0
          }
          if(runif(1) < exp((llysightcandsum+sum(ll.tel.cand)+prior.cand)-(llysightsum+sum(ll.tel)+prior.curr))){
            sigma<- sigma.cand
            lamd.sight.marked=lamd.sight.marked.cand
            lamd.sight.unmarked=lamd.sight.unmarked.cand
            ll.y.sight.marked=ll.y.sight.marked.cand
            ll.y.sight.unmarked=ll.y.sight.unmarked.cand
            ll.tel=ll.tel.cand
          }
        }
      }
      
      # ID update
      if(IDup=="Gibbs"){
        #Update y.sight.true from full conditional canceling out inconsistent combos with constraints.
        up=sample(1:n.samp.latent,nswap,replace=FALSE)
        for(l in up){
          nj=which(y.sight.latent[l,]>0)
          #Can only swap if IDcovs match
          idx=which(G.use[l,]!=0)
          if(length(idx)>1){#multiple loci observed
            possible.M=which(z1==1&apply(G.true.marked[,idx],1,function(x){all(x==G.use[l,idx])}))
            possible.UM=which(z2==1&apply(G.true.unmarked[,idx],1,function(x){all(x==G.use[l,idx])}))
          }else if(length(idx)==1){#single loci observed
            possible.M=which(z1==1&G.true.marked[,idx]==G.use[l,idx])
            possible.UM=which(z2==1&G.true.unmarked[,idx]==G.use[l,idx])
          }else{#fully latent G.obs
            possible.M=which(z1==1)#Can match anyone
            possible.UM=which(z2==1)#Can match anyone
          }
          if(!(useUnk|useMarkednoID)){#mark status exclusions handled through G.true
            possible.M=c() #Can't swap to a marked guy
          }else{
            if(Mark.obs[l]==2){#This is an unmarked sample
              possible.M=c()#Can't swap to a marked guy
            }
            if(Mark.obs[l]==1){#This is a marked sample
              possible.UM=c()#Can't swap to an unmarked guy
            }
          }
          if(length(c(possible.M,possible.UM))==0)next
          njprobs.M=lamd.sight.marked[,nj]
          njprobs.M[setdiff(1:M1,possible.M)]=0
          njprobs.UM=lamd.sight.unmarked[,nj]
          njprobs.UM[setdiff(1:M2,possible.UM)]=0
          njprobs=c(njprobs.M,njprobs.UM)
          njprobs=njprobs/sum(njprobs)
          
          newID=sample(1:(M1+M2),1,prob=njprobs)
          if(newID<=M1){
            newID=c(newID,1)
          }else{
            newID=c(newID-M1,2)
          }
          
          if(!all(ID[l,]==newID)){ #did you change ID number or marked/unmarked status?
            swapped=c(ID[l,1],newID[1])
            if(ID[l,2]==1&newID[2]==1){
              samptype="M2M"
            }else if(ID[l,2]==1&newID[2]==2){
              samptype="M2UM"
            }else if(ID[l,2]==2&newID[2]==2){
              samptype="UM2UM"
            }else if(ID[l,2]==2&newID[2]==1){
              samptype="UM2M"
            }
            if(samptype=="M2M"){#marked guy
              y.sight.marked.true[ID[l,1],]=y.sight.marked.true[ID[l,1],]-y.sight.latent[l,]
              y.sight.marked.true[newID[1],]=y.sight.marked.true[newID[1],]+y.sight.latent[l,]
              if(obstype=="bernoulli"){
                ll.y.sight.marked[swapped,]= dbinom(y.sight.marked.true[swapped,],K2D.M[swapped,],pd.sight.marked[swapped,],log=TRUE)
              }else{
                ll.y.sight.marked[swapped,]= dpois(y.sight.marked.true[swapped,],K2D.M[swapped,]*lamd.sight.marked[swapped,],log=TRUE)
              }
            }else if(samptype=="UM2UM"){#unmarked guy
              y.sight.unmarked.true[ID[l,1],]=y.sight.unmarked.true[ID[l,1],]-y.sight.latent[l,]
              y.sight.unmarked.true[newID[1],]=y.sight.unmarked.true[newID[1],]+y.sight.latent[l,]
              if(obstype=="bernoulli"){
                ll.y.sight.unmarked[swapped,]= dbinom(y.sight.unmarked.true[swapped,],K2D.UM[swapped,],pd.sight.unmarked[swapped,],log=TRUE)
              }else{
                ll.y.sight.unmarked[swapped,]= dpois(y.sight.unmarked.true[swapped,],K2D.UM[swapped,]*lamd.sight.unmarked[swapped,],log=TRUE)
              }
            }else if(samptype=="M2UM"){
              y.sight.marked.true[ID[l,1],]=y.sight.marked.true[ID[l,1],]-y.sight.latent[l,]
              y.sight.unmarked.true[newID[1],]=y.sight.unmarked.true[newID[1],]+y.sight.latent[l,]
              if(obstype=="bernoulli"){
                ll.y.sight.marked[swapped[1],]= dbinom(y.sight.marked.true[swapped[1],],K2D.M[swapped[1],],pd.sight.marked[swapped[1],],log=TRUE)
                ll.y.sight.unmarked[swapped[2],]= dbinom(y.sight.unmarked.true[swapped[2],],K2D.UM[swapped[2],],pd.sight.unmarked[swapped[2],],log=TRUE)
              }else{
                ll.y.sight.marked[swapped[1],]= dpois(y.sight.marked.true[swapped[1],],K2D.M[swapped[1],]*lamd.sight.marked[swapped[1],],log=TRUE)
                ll.y.sight.unmarked[swapped[2],]= dpois(y.sight.unmarked.true[swapped[2],],K2D.UM[swapped[2],]*lamd.sight.unmarked[swapped[2],],log=TRUE)
              }
            }else if(samptype=="UM2M"){
              y.sight.unmarked.true[ID[l,1],]=y.sight.unmarked.true[ID[l,1],]-y.sight.latent[l,]
              y.sight.marked.true[newID[1],]=y.sight.marked.true[newID[1],]+y.sight.latent[l,]
              if(obstype=="bernoulli"){
                ll.y.sight.marked[swapped[2],]= dbinom(y.sight.marked.true[swapped[2],],K2D.M[swapped[2],],pd.sight.marked[swapped[2],],log=TRUE)
                ll.y.sight.marked[swapped[1],]= dbinom(y.sight.unmarked.true[swapped[1],],K2D.UM[swapped[1],],pd.sight.unmarked[swapped[1],],log=TRUE)
              }else{
                ll.y.sight.marked[swapped[2],]= dpois(y.sight.marked.true[swapped[2],],K2D.M[swapped[2],]*lamd.sight.marked[swapped[2],],log=TRUE)
                ll.y.sight.unmarked[swapped[1],]= dpois(y.sight.unmarked.true[swapped[1],],K2D.UM[swapped[1],]*lamd.sight.unmarked[swapped[1],],log=TRUE)
              }
            }
            ID[l,]=newID
          }
        }
      }else{#MH update
        up=sample(1:n.samp.latent,nswap,replace=FALSE)
        for(l in up){
          nj=which(y.sight.latent[l,]>0)
          #Can only swap if IDcovs match
          idx=which(G.use[l,]!=0)
          if(length(idx)>1){#multiple loci observed
            possible.M=which(z1==1&apply(G.true.marked[,idx],1,function(x){all(x==G.use[l,idx])}))
            possible.UM=which(z2==1&apply(G.true.unmarked[,idx],1,function(x){all(x==G.use[l,idx])}))
          }else if(length(idx)==1){#single loci observed
            possible.M=which(z1==1&G.true.marked[,idx]==G.use[l,idx])
            possible.UM=which(z2==1&G.true.unmarked[,idx]==G.use[l,idx])
          }else{#fully latent G.obs
            possible.M=which(z1==1)#Can match anyone
            possible.UM=which(z2==1)#Can match anyone
          }
          if(!(useUnk|useMarkednoID)){#mark status exclusions handled through G.true
            possible.M=c() #Can't swap to a marked guy
          }else{
            if(Mark.obs[l]==2){#This is an unmarked sample
              possible.M=c()#Can't swap to a marked guy
            }
            if(Mark.obs[l]==1){#This is a marked sample
              possible.UM=c()#Can't swap to an unmarked guy
            }
          }
          if(length(c(possible.M,possible.UM))==0)next
          njprobs.M=lamd.sight.marked[,nj]
          njprobs.M[setdiff(1:M1,possible.M)]=0
          njprobs.UM=lamd.sight.unmarked[,nj]
          njprobs.UM[setdiff(1:M2,possible.UM)]=0
          njprobs=c(njprobs.M,njprobs.UM)
          njprobs=njprobs/sum(njprobs)
          newID=ID
          newID[l,1]=sample(1:(M1+M2),1,prob=njprobs)
          if(newID[l,1]<=M1){
            newID[l,]=c(newID[l,1],1)
          }else{
            newID[l,]=c(newID[l,1]-M1,2)
          }
          if(all(ID[l,]==newID[l,]))next #did you change ID number or marked/unmarked status?
          swapped=c(ID[l,1],newID[l,1])
          #proposal probabilities for proposing to this individual
          swapped2=swapped
          if(ID[l,2]==2){
            swapped2[1]=swapped2[1]+M1
          }
          if(newID[l,2]==2){
            swapped2[2]=swapped2[2]+M1
          }
          
          propprob=njprobs[swapped2[2]]
          backprob=njprobs[swapped2[1]]
          if(ID[l,2]==1&newID[l,2]==1){
            samptype="M2M"
          }else if(ID[l,2]==1&newID[l,2]==2){
            samptype="M2UM"
          }else if(ID[l,2]==2&newID[l,2]==2){
            samptype="UM2UM"
          }else if(ID[l,2]==2&newID[l,2]==1){
            samptype="UM2M"
          }
          if(samptype=="M2M"){#marked guy
            y.sight.marked.cand[ID[l,1],]=y.sight.marked.true[ID[l,1],]-y.sight.latent[l,]
            y.sight.marked.cand[newID[l,1],]=y.sight.marked.true[newID[l,1],]+y.sight.latent[l,]
            if(obstype=="bernoulli"){
              ll.y.sight.marked.cand[swapped,]= dbinom(y.sight.marked.cand[swapped,],K2D.M[swapped,],pd.sight.marked[swapped,],log=TRUE)
            }else{
              ll.y.sight.marked.cand[swapped,]= dpois(y.sight.marked.cand[swapped,],K2D.M[swapped,]*lamd.sight.marked[swapped,],log=TRUE)
            }
            llysum=sum(ll.y.sight.marked[swapped,])
            llycandsum=sum(ll.y.sight.marked.cand[swapped,])
            #proposal probabilities for selecting this sample to update
            focalprob=(sum(ID[l,1]==ID[,1]&ID[l,2]==ID[,2])/n.samp.latent)*(y.sight.marked.true[ID[l,1],nj]/sum(y.sight.marked.true[ID[l,1],]))
            focalbackprob=(sum(newID[l,1]==newID[,1]&newID[l,2]==newID[,2])/n.samp.latent)*(y.sight.marked.cand[newID[l,1],nj]/sum(y.sight.marked.cand[newID[l,1],]))
            
          }else if(samptype=="UM2UM"){#unmarked guy
            y.sight.unmarked.cand[ID[l,1],]=y.sight.unmarked.true[ID[l,1],]-y.sight.latent[l,]
            y.sight.unmarked.cand[newID[l,1],]=y.sight.unmarked.true[newID[l,1],]+y.sight.latent[l,]
            if(obstype=="bernoulli"){
              ll.y.sight.unmarked.cand[swapped,]= dbinom(y.sight.unmarked.cand[swapped,],K2D.UM[swapped,],pd.sight.unmarked[swapped,],log=TRUE)
            }else{
              ll.y.sight.unmarked.cand[swapped,]= dpois(y.sight.unmarked.cand[swapped,],K2D.UM[swapped,]*lamd.sight.unmarked[swapped,],log=TRUE)
            }
            llysum=sum(ll.y.sight.unmarked[swapped,])
            llycandsum=sum(ll.y.sight.unmarked.cand[swapped,])
            #proposal probabilities for selecting this sample to update
            focalprob=(sum(ID[l,1]==ID[,1]&ID[l,2]==ID[,2])/n.samp.latent)*(y.sight.unmarked.true[ID[l,1],nj]/sum(y.sight.unmarked.true[ID[l,1],]))
            focalbackprob=(sum(newID[l,1]==newID[,1]&newID[l,2]==newID[,2])/n.samp.latent)*(y.sight.unmarked.cand[newID[l,1],nj]/sum(y.sight.unmarked.cand[newID[l,1],]))
          }else if(samptype=="M2UM"){
            
            y.sight.marked.cand[ID[l,1],]=y.sight.marked.true[ID[l,1],]-y.sight.latent[l,]
            y.sight.unmarked.cand[newID[l,1],]=y.sight.unmarked.true[newID[l,1],]+y.sight.latent[l,]
            if(obstype=="bernoulli"){
              ll.y.sight.marked.cand[swapped[1],]= dbinom(y.sight.marked.cand[swapped[1],],K2D.M[swapped[1],],pd.sight.marked[swapped[1],],log=TRUE)
              ll.y.sight.unmarked.cand[swapped[2],]= dbinom(y.sight.unmarked.cand[swapped[2],],K2D.UM[swapped[2],],pd.sight.unmarked[swapped[2],],log=TRUE)
            }else{
              ll.y.sight.marked.cand[swapped[1],]= dpois(y.sight.marked.cand[swapped[1],],K2D.M[swapped[1],]*lamd.sight.marked[swapped[1],],log=TRUE)
              ll.y.sight.unmarked.cand[swapped[2],]= dpois(y.sight.unmarked.cand[swapped[2],],K2D.UM[swapped[2],]*lamd.sight.unmarked[swapped[2],],log=TRUE)
            }
            llysum=sum(ll.y.sight.marked[swapped[1],])+sum(ll.y.sight.unmarked[swapped[2],])
            llycandsum=sum(ll.y.sight.marked.cand[swapped[1],])+sum(ll.y.sight.unmarked.cand[swapped[2],])
            #proposal probabilities for selecting this sample to update
            focalprob=(sum(ID[l,1]==ID[,1]&ID[l,2]==ID[,2])/n.samp.latent)*(y.sight.marked.true[ID[l,1],nj]/sum(y.sight.marked.true[ID[l,1],]))
            focalbackprob=(sum(newID[l,1]==newID[,1]&newID[l,2]==newID[,2])/n.samp.latent)*(y.sight.unmarked.cand[newID[l,1],nj]/sum(y.sight.unmarked.cand[newID[l,1],]))
          }else if(samptype=="UM2M"){
            
            y.sight.unmarked.cand[ID[l,1],]=y.sight.unmarked.true[ID[l,1],]-y.sight.latent[l,]
            y.sight.marked.cand[newID[l,1],]=y.sight.marked.true[newID[l,1],]+y.sight.latent[l,]
            if(obstype=="bernoulli"){
              ll.y.sight.marked.cand[swapped[2],]= dbinom(y.sight.marked.cand[swapped[2],],K2D.M[swapped[2],],pd.sight.marked[swapped[2],],log=TRUE)
              ll.y.sight.unmarked.cand[swapped[1],]= dbinom(y.sight.unmarked.cand[swapped[1],],K2D.UM[swapped[1],],pd.sight.unmarked[swapped[1],],log=TRUE)
            }else{
              ll.y.sight.marked.cand[swapped[2],]= dpois(y.sight.marked.cand[swapped[2],],K2D.M[swapped[2],]*lamd.sight.marked[swapped[2],],log=TRUE)
              ll.y.sight.unmarked.cand[swapped[1],]= dpois(y.sight.unmarked.cand[swapped[1],],K2D.UM[swapped[1],]*lamd.sight.unmarked[swapped[1],],log=TRUE)
            }
            llysum=sum(ll.y.sight.marked[swapped[2],])+sum(ll.y.sight.unmarked[swapped[1],])
            llycandsum=sum(ll.y.sight.marked.cand[swapped[2],])+sum(ll.y.sight.unmarked.cand[swapped[1],])
            #proposal probabilities for selecting this sample to update
            focalprob=(sum(ID[l,1]==ID[,1]&ID[l,2]==ID[,2])/n.samp.latent)*(y.sight.unmarked.true[ID[l,1],nj]/sum(y.sight.unmarked.true[ID[l,1],]))
            focalbackprob=(sum(newID[l,1]==newID[,1]&newID[l,2]==newID[,2])/n.samp.latent)*(y.sight.marked.cand[newID[l,1],nj]/sum(y.sight.marked.cand[newID[l,1],]))
          }
          
          if(runif(1)<exp(llycandsum-llysum)*(backprob/propprob)*(focalbackprob/focalprob)){
            if(samptype=="M2M"){#marked guy
              y.sight.marked.true[swapped,]=y.sight.marked.cand[swapped,]
              ll.y.sight.marked[swapped,]=ll.y.sight.marked.cand[swapped,]
            }else if(samptype=="UM2UM"){#unmarked guy
              y.sight.unmarked.true[swapped,]=y.sight.unmarked.cand[swapped,]
              ll.y.sight.unmarked[swapped,]=ll.y.sight.unmarked.cand[swapped,]
            }else if(samptype=="M2UM"){
              y.sight.marked.true[ID[l,1],]=y.sight.marked.cand[ID[l,1],]
              y.sight.unmarked.true[newID[l,1],]=y.sight.unmarked.cand[newID[l,1],]
              ll.y.sight.marked[swapped[1],]=ll.y.sight.marked.cand[swapped[1],]
              ll.y.sight.unmarked[swapped[2],]=ll.y.sight.unmarked.cand[swapped[2],]
            }else if(samptype=="UM2M"){
              y.sight.unmarked.true[ID[l,1],]=y.sight.unmarked.cand[ID[l,1],]
              y.sight.marked.true[newID[l,1],]=y.sight.marked.cand[newID[l,1],]
              ll.y.sight.marked[swapped[2],]= ll.y.sight.marked.cand[swapped[2],]
              ll.y.sight.unmarked[swapped[1],]=ll.y.sight.unmarked.cand[swapped[1],]
            }
            ID[l,]=newID[l,]
          }
        }
      }
      #update known.vector
      known.vector.marked=1*(rowSums(y.sight.marked.true)>0)
      known.vector.unmarked=1*(rowSums(y.sight.unmarked.true)>0)
      
      #Update G.true.marked
      G.true.marked.tmp=matrix(0, nrow=M1,ncol=ncat)
      G.true.marked.tmp[marked.guys,]=1
      for(i in unique(ID[ID[,2]==1,1])){
        idX=which(ID[,1]==i&ID[,2]==1)
        if(length(idX)==1){
          G.true.marked.tmp[i,]=G.use[idX,]
        }else{
          if(ncol(G.use)>1){
            G.true.marked.tmp[i,]=apply(G.use[idX,],2, max) #consensus
          }else{
            G.true.marked.tmp[i,]=max(G.use[idX,]) #consensus
          }
        }
      }
      G.latent.marked=G.true.marked.tmp==0
      for(j in 1:ncat){
        swap=G.latent.marked[,j]
        G.true.marked[swap,j]=sample(IDcovs[[j]],sum(swap),replace=TRUE,prob=gamma[[j]])
      }
      
      #Update G.true.unmarked
      G.true.unmarked.tmp=matrix(0, nrow=M2,ncol=ncat)
      for(i in unique(ID[ID[,2]==2,1])){
        idX=which(ID[,1]==i&ID[,2]==2)
        if(length(idX)==1){
          G.true.unmarked.tmp[i,]=G.use[idX,]
        }else{
          if(ncol(G.use)>1){
            G.true.unmarked.tmp[i,]=apply(G.use[idX,],2, max) #consensus
          }else{
            G.true.unmarked.tmp[i,]=max(G.use[idX,]) #consensus
          }
        }
      }
      G.latent.unmarked=G.true.unmarked.tmp==0
      for(j in 1:ncat){
        swap=G.latent.unmarked[,j]
        G.true.unmarked[swap,j]=sample(IDcovs[[j]],sum(swap),replace=TRUE,prob=gamma[[j]])
      }
      
      #update genotype frequencies
      for(j in 1:ncat){
        x=rep(NA,nIDcovs[[j]])
        for(k in 1:nIDcovs[[j]]){
          x[k]=sum(G.true.marked[z1==1,j]==k)+sum(G.true.unmarked[z2==1,j]==k)#genotype freqs in pop
        }
        gam=rgamma(rep(1,nIDcovs[[j]]),1+x)
        gamma[[j]]=gam/sum(gam)
      }
      
      ## probability of not being captured in a trap AT ALL by either method
      if(obstype=="poisson"){
        pd.sight.marked=1-exp(-lamd.sight.marked)
        pd.sight.unmarked=1-exp(-lamd.sight.unmarked)
      }
      pbar.sight.marked=(1-pd.sight.marked)^K2D.M
      pbar.sight.unmarked=(1-pd.sight.unmarked)^K2D.UM
      prob0.sight.marked= exp(rowSums(log(pbar.sight.marked)))
      prob0.sight.unmarked= exp(rowSums(log(pbar.sight.unmarked)))
      fc.marked<- prob0.sight.marked*psi1/(prob0.sight.marked*psi1 + 1-psi1)
      fc.unmarked<- prob0.sight.unmarked*psi2/(prob0.sight.unmarked*psi2 + 1-psi2)
      z1[known.vector.marked==0]<- rbinom(sum(known.vector.marked ==0), 1, fc.marked[known.vector.marked==0])
      z2[known.vector.unmarked==0]<- rbinom(sum(known.vector.unmarked ==0), 1, fc.unmarked[known.vector.unmarked==0])
      if(obstype=="bernoulli"){
        ll.y.sight.marked= dbinom(y.sight.marked.true,K2D.M,pd.sight.marked*z1,log=TRUE)
        ll.y.sight.unmarked= dbinom(y.sight.unmarked.true,K2D.UM,pd.sight.unmarked*z2,log=TRUE)
      }else{
        ll.y.sight.marked= dpois(y.sight.marked.true,K2D.M*lamd.sight.marked*z1,log=TRUE)
        ll.y.sight.unmarked= dpois(y.sight.unmarked.true,K2D.UM*lamd.sight.unmarked*z2,log=TRUE)
      }
      psi1=rbeta(1,1+sum(z1),1+M1-sum(z1))
      psi2=rbeta(1,1+sum(z2),1+M2-sum(z2))
      
      ## Now we have to update the activity centers
      for (i in 1:M1) {#marked guys first
        if(i%in%telguys){
          Scand <- c(rnorm(1, s1[i, 1], proppars$st), rnorm(1, s1[i, 2], proppars$st))
        }else{
          Scand <- c(rnorm(1, s1[i, 1], proppars$s), rnorm(1, s1[i, 2], proppars$s))
        }
        if(useverts==FALSE){
          inbox <- Scand[1] < xlim[2] & Scand[1] > xlim[1] & Scand[2] < ylim[2] & Scand[2] > ylim[1]
        }else{
          inbox=inout(Scand,vertices)
        }
        if (inbox) {
          dtmp <- sqrt((Scand[1] - X[, 1])^2 + (Scand[2] - X[, 2])^2)
          lamd.sight.marked.cand[i,]<- lam0*exp(-dtmp*dtmp/(2*sigma*sigma))
          if(obstype=="bernoulli"){
            pd.sight.marked.cand[i,]=1-exp(-lamd.sight.marked.cand[i,])
            ll.y.sight.marked.cand[i,]= dbinom(y.sight.marked.true[i,],K2D.M[i,],pd.sight.marked.cand[i,]*z1[i],log=TRUE)
            if(uselocs&(i%in%telguys)){
              ll.tel.cand[i,]=dnorm(locs[i,,1],Scand[1],sigma,log=TRUE)+dnorm(locs[i,,2],Scand[2],sigma,log=TRUE)
              if (runif(1) < exp((sum(ll.y.sight.marked.cand[i,])+sum(ll.tel.cand[i,])) -
                                 (sum(ll.y.sight.marked[i,])+sum(ll.tel[i,])))) {
                s1[i,]=Scand
                D1[i,]=dtmp
                lamd.sight.marked[i,]=lamd.sight.marked.cand[i,]
                pd.sight.marked[i,]=pd.sight.marked.cand[i,]
                ll.y.sight.marked[i,]=ll.y.sight.marked.cand[i,]            
                ll.tel[i,]=ll.tel.cand[i,]
              }
            }else{
              if (runif(1) < exp(sum(ll.y.sight.marked.cand[i,]) -sum(ll.y.sight.marked[i,]))) {
                s1[i,]=Scand
                D1[i,]=dtmp
                lamd.sight.marked[i,]=lamd.sight.marked.cand[i,]
                pd.sight.marked[i,]=pd.sight.marked.cand[i,]
                ll.y.sight.marked[i,]=ll.y.sight.marked.cand[i,]            
              }
            }
          }else{#poisson
            ll.y.sight.marked.cand[i,]= dpois(y.sight.marked.true[i,],K2D.M[i,]*lamd.sight.marked.cand[i,]*z1[i],log=TRUE)
            if(uselocs&(i%in%telguys)){
              ll.tel.cand[i,]=dnorm(locs[i,,1],Scand[1],sigma,log=TRUE)+dnorm(locs[i,,2],Scand[2],sigma,log=TRUE)
              if (runif(1) < exp((sum(ll.y.sight.marked.cand[i,])+sum(ll.tel.cand[i,])) -
                                 (sum(ll.y.sight.marked[i,])+sum(ll.tel[i,])))) {
                s1[i,]=Scand
                D1[i,]=dtmp
                lamd.sight.marked[i,]=lamd.sight.marked.cand[i,]
                ll.y.sight.marked[i,]=ll.y.sight.marked.cand[i,]            
                ll.tel[i,]=ll.tel.cand[i,]
              }
            }else{
              if (runif(1) < exp(sum(ll.y.sight.marked.cand[i,]) -sum(ll.y.sight.marked[i,]))) {
                s1[i,]=Scand
                D1[i,]=dtmp
                lamd.sight.marked[i,]=lamd.sight.marked.cand[i,]
                ll.y.sight.marked[i,]=ll.y.sight.marked.cand[i,]            
              }
            }
          }
        }
      }
      ## Now unmarked guys
      for (i in 1:M2) {
        Scand <- c(rnorm(1, s2[i, 1], proppars$s), rnorm(1, s2[i, 2], proppars$s))
        if(useverts==FALSE){
          inbox <- Scand[1] < xlim[2] & Scand[1] > xlim[1] & Scand[2] < ylim[2] & Scand[2] > ylim[1]
        }else{
          inbox=inout(Scand,vertices)
        }
        if (inbox) {
          dtmp <- sqrt((Scand[1] - X[, 1])^2 + (Scand[2] - X[, 2])^2)
          lamd.sight.unmarked.cand[i,]<- lam0*exp(-dtmp*dtmp/(2*sigma*sigma))
          if(obstype=="bernoulli"){
            pd.sight.unmarked.cand[i,]=1-exp(-lamd.sight.unmarked.cand[i,])
            ll.y.sight.unmarked.cand[i,]= dbinom(y.sight.unmarked.true[i,],K2D.UM[i,],pd.sight.unmarked.cand[i,]*z2[i],log=TRUE)
            if (runif(1) < exp(sum(ll.y.sight.unmarked.cand[i,]) -sum(ll.y.sight.unmarked[i,]))) {
              s2[i,]=Scand
              D2[i,]=dtmp
              lamd.sight.unmarked[i,]=lamd.sight.unmarked.cand[i,]
              pd.sight.unmarked[i,]=pd.sight.unmarked.cand[i,]
              ll.y.sight.unmarked[i,]=ll.y.sight.unmarked.cand[i,]            
            }
          }else{#poisson
            ll.y.sight.unmarked.cand[i,]= dpois(y.sight.unmarked.true[i,],K2D.UM[i,]*lamd.sight.unmarked.cand[i,]*z2[i],log=TRUE)
            
            if (runif(1) < exp(sum(ll.y.sight.unmarked.cand[i,]) -sum(ll.y.sight.unmarked[i,]))) {
              s2[i,]=Scand
              D2[i,]=dtmp
              lamd.sight.unmarked[i,]=lamd.sight.unmarked.cand[i,]
              ll.y.sight.unmarked[i,]=ll.y.sight.unmarked.cand[i,]            
            }
          }
          
        }
      }
      #Do we record output on this iteration?
      if(iter>nburn&iter%%nthin==0){
        if(storeLatent){
          s1xout[iteridx,]<- s1[,1]
          s1yout[iteridx,]<- s1[,2]
          z1out[iteridx,]<- z1
          s2xout[iteridx,]<- s2[,1]
          s2yout[iteridx,]<- s2[,2]
          z2out[iteridx,]<- z2
          IDout[iteridx,,]=ID
        }
        if(storeGamma){
          for(k in 1:ncat){
            gammaOut[[k]][iteridx,]=gamma[[k]]
          }
        }
        NM=sum(z1)
        NUM=sum(z2)
        nM=length(unique(ID[ID[,2]==1,1])  )
        nUM=length(unique(ID[ID[,2]==2,1]))
        ntot=sum(rowSums(y.sight.marked.true)>0)+sum(rowSums(y.sight.unmarked.true)>0)
        out[iteridx,]<- c(lam0,sigma,nM,nUM,ntot,NM,NUM,NM+NUM,psi1,psi2)
        iteridx=iteridx+1
      }
    }  # end of MCMC algorithm
    
    if(storeLatent&storeGamma){
      list(out=out, s1xout=s1xout, s1yout=s1yout, z1out=z1out,s2xout=s2xout, s2yout=s2yout, z2out=z2out,IDout=IDout,gammaOut=gammaOut)
    }else if(storeLatent&!storeGamma){
      list(out=out, s1xout=s1xout, s1yout=s1yout, z1out=z1out,s2xout=s2xout, s2yout=s2yout, z2out=z2out,IDout=IDout)
    }else if(!storeLatent&storeGamma){
      list(out=out,gammaOut=gammaOut)
    }else{ 
      list(out=out)
    }
  }
benaug/SPIM documentation built on April 28, 2024, 7:27 a.m.