### SIMULATE a DATASET under the CALIBRATION DESIGN
rm(list = ls())	# Clean the workspace
set.seed(782)

### PARAMETERS
psi <- 0.6	# occupancy probability 
p11 <- 0.7	# site-level detection probability (conditional on occupancy status)
p10 <- 0.1	# site-level false positive probability (conditional on occupancy status)

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

n1 <- 50 	# number of reference sites OCCUPIED
n0 <- 70	# number of reference sites UNOCCUPIED

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

### SIMULATE AMBIGUOUS-DETECTION/NON-DETECTION DATA Y
 Y <- matrix(NA,R,J)			### Create a matrix to store detection/non-detection data 
for (i in 1:R){				### Loop over all sites 'i'
 p <- Z[i]*p11 + (1-Z[i])*p10		### Define the 'binomial' parameter, conditional on site occupancy status 
  Y[i,] <- rbinom(J, 1, p)		### Site-level data Y [0/1]
}	  					### End Loop over sites
Y.sum <- apply(Y, 1, sum)		### Summarize the data as 'number of detection for each site'

### Evaluation of the Detection Method
X <- rbinom(1, n1, p11)			### Number of (true) detections in reference sites OCCUPIED	
W <- rbinom(1, n0, p10)			### Number of (true) detections in reference sites UNOCCUPIED

### Save the data used in the model
dataField.Cal <- Y.sum				### Ambiguous Field data 
dataLab.Cal <- rbind( c(X,n1),c(W,n0) )	### Known-truth data from reference sites (calibration method)
occ.Cal <- J					### Number of survey occasions
 TrueValues.Cal <- cbind(psi, p11, p10)	### Save the true parameter values for comparison wih estimates

# Save these obsjects as a 'RData' file
save(file = "DataCal.RData", list = c("dataField.Cal", "dataLab.Cal", "occ.Cal", "TrueValues.Cal"))