### SIMULATE a DATASET under the SITE CONFIRMATION DESIGN
### MULTIPLE-DETECTION-STATES
rm(list = ls())	# Clean the workspace
set.seed(741)

### POPULATION PARAMETERS
psi <- 0.6	# occupancy probability
p11 <- 0.7	# site-level detection probability (from ambiguous method)
p10 <- 0.1	# site-level false positive probability (from ambiguous method)
b <- 0.2	# site-level detection probability (from unambiguous method)

### SAMPLE CHARACTERSITICS
 K <- 5	# number of survey occasions
R <- 500	# number of sites

### OCCUPANCY STATES
Z <- rbinom(R, 1, psi)

### SIMULATE MULTI-STATE-DETECTION DATA (data Y)
 Y <- matrix(NA,R,3)				### Create a matrix to store the mutlistate observation data 
  for (i in 1:R){					### Loop over all sites 'i'
pr.X0 <- Z[i]*(1-p11) + (1-Z[i])*(1-p10)	### Define the probability of observation-state y = 0 (non-detection)
pr.X1 <- Z[i]*(1-b)*p11 + (1-Z[i])*p10	### Define the probability of observation-state y = 1 (ambiguous-detection)
pr.X2 <- Z[i]*b*p11				### Define the probability of observation-state y = 2 (unambiguous-detection)
  pi <- c( pr.X0 , pr.X1 , pr.X2 )		### Define the vector of probabilities of the Categorical distribution
Y[i,] <- as.vector(rmultinom(1,K,pi))	### Site-level data Y [mutli-state]
   }							### End Loop over sites

### Save the data used in the model
data.SiteConf.MS <- Y						### Observation DATA
occ.SiteConf.MS <- K						### number of survey occasions
TrueValues.SiteConf.MS <- cbind(psi, p11, p10, b)	### Save the true parameter values for comparison wih estimates

# Save these obsjects as a 'RData' file
save(file = "DataSiteConf-MS.RData", list = c("data.SiteConf.MS", "occ.SiteConf.MS", "TrueValues.SiteConf.MS"))