# Functions and helper functions to simulate mark recapture data

# example of simulating data from model Mb, then computing
#   summary statistics for Mb, M0, and Mt models


# simulate data under model Mb with N=100, p=0.3, c=0.5

# obs <- simMb3(c(100,0.3,0.5))

# your estimation functions probably want the sufficient statistics:
#    they are calculated by the stats functions

#  obs.Mb <- statsMb3(obs)	# gives t, Mt+1, Mdot and mdot 
#  obs.M0 <- statsM03(obs)	# gives t, Mt+1, ndot
#  obs.Mt <- statsMt3(obs)	# gives t, Mt+1, n1, n2, n3

# you should then use obs.Mb, obs.M0, and obs.Mt to estimate the 
#   mle's for the three models.


rmulti <- function(p,totn) {
# generate one realization of a multinomial with the specified 
#   vector of p's and a total count of n

nclass <- length(p)
cump <- cumsum(p)
cump[nclass] <- 1	# eliminate roundoff errors

x <- runif(totn)
mx <- outer(cump,x,'<')
obs <- apply(mx,2,sum)
table(c(0:(nclass-1),obs))-1
}

simM03 <- function(N, p) {
# simulate one realization of data with constant P[capture],
#    3 capture occasions

# output is list with two components:
# history: vector of counts for 7 observed capture histories
#  in the order: YYY, YYN, YNY, YNN, NYY, NYN, NNY
# stats: sufficient statistics: vector with t, Mt+1, ndot

Phistory <- c(p^3, p^2*(1-p), p^2*(1-p), p*(1-p)^2, (1-p)*p^2, 
		(1-p)^2*p, (1-p)^2*p, (1-p)^3)
		
Nhistory <- rmulti(Phistory, N)		# simulate a random data set
names(Nhistory) <- names(Phistory)

Nhistory
}

statsM03 <- function(Nhistory) {
# compute summary statistics for M0 model
# input: counts for each capture history
# output: vector of summary statistics (t=3, Mt+1, ndot)

Mt1 <- sum(Nhistory[1:7])		# count number seen at least once
ndot <- 3*Nhistory[1]+2*sum(Nhistory[c(2,3,5)])+sum(Nhistory[c(4,6,7)])
	# and total number of times individuals are seen
c(t=3,Mt1=Mt1,ndot=ndot)
}

simMt3 <- function(N, p) {
# simulate one realization of data with time varying P[capture]
#   3 capture occasions

# output is vector of counts for 7 observed capture histories
#  in the order: YYY, YYN, YNY, YNN, NYY, NYN, NNY

Phistory <- c(p[1]*p[2]*p[3], p[1]*p[2]*(1-p[3]), p[1]*(1-p[2])*p[3],
		p[1]*(1-p[2])*(1-p[3]), (1-p[1])*p[2]*p[3],
		(1-p[1])*p[2]*(1-p[3]), (1-p[1])*(1-p[2])*p3,
		(1-p[1])*(1-p[2])*(1-p3) )

Nhistory <- rmulti(Phistory,N)
names(Nhistory) <- names(Phistory)

Nhistory
}

statsMt3 <- function(Nhistory) {
# input: counts for each capture history
# output: vector of summary statistics (t=3, Mt+1, n1, n2, n3

Mt1 <- sum(Nhistory[1:7])
n1 <- sum(Nhistory[1:4])
n2 <- sum(Nhistory[c(1,2,5,6)])
n3 <- sum(Nhistory[c(1,3,5,7)])

c(t=3,Mt1=Mt1,n1=n1,n2=n2,n3=n3)
}
		
				
simMb3 <- function(N,p,c) {
# simulate one realization of data with behavioral 
#   heterogeneity, 3 capture occasions

# output is vector of counts for 7 observed capture histories
#    in the order: YYY, YYN, YNY, YNN, NYY, NYN, NNY


Phistory <- c(YYY=p*c^2, YYN=p*c*(1-c), YNY=p*(1-c)*c, YNN=p*(1-c)^2, 
		NYY=(1-p)*p*c, NYN=(1-p)*p*(1-c), NNY=(1-p)^2*p, 
		NNN=(1-p)^3)
		
Nhistory <- rmulti(Phistory, N)		# simulate a random data set
names(Nhistory) <- names(Phistory)

Nhistory
}

statsMb3 <- function(Nhistory) {
# compute summary statistics for 3 period behavioural het. model
# input: counts for each capture history
# output: vector of summary statistics (t=3, Mt+1, Mdot, mdot)

Mt1 <- sum(Nhistory[1:7])		# count number seen at least once
M2 <- sum(Nhistory[1:4])		# marked at time 2
M3 <- sum(Nhistory[1:6])		# marked at time 3
m2 <- Nhistory[1]+Nhistory[2]		# marked caught at time 2
m3 <- Nhistory[1]+Nhistory[3]+Nhistory[5]	# ditto at time 3

c(t=3, Mt1=Mt1, Mdot=M2+M3, mdot=m2+m3)
}