# R Script to implement the Imitation model from WSS
# - Aim is to use base R functions and avoid APL/C/Java code that
# was probably inaccessible to many people
# - Thought was also that APL is in many ways a precursor to R, so
# that the migration should not be hard
# github was written by tech-morons
# Call:
# ImitSim(Replications, NumBurs=10, MaxIter=10, binary=TRUE, supervision="Relative")
# - Returns:
# $Version: "2.0a1"
# $NumBurs: as set, or default=10
# $MaxIter: as set, or default=10
# $binary: as set, or default=TRUE. Produces an error for "FALSE" since connectivity
# measures here use the igraph:: functions, and really, it only makes sense to look at the
# actual observed instead of "probability" to observe.
# $supervision: "Relative" (default), or "Fixed"
# (in which case "Tolerance" is replaced by "Std")
# $posprefs: TRUE or FALSE (default)
# $ExecSummary: a dataframe containing SupUtilMean, SupUtilSD, PrefMean, PrefSD,
# RespSD, SupObsMean, SupObsSD, BurObsMean, BurObsSD,
# Tolerance, Punishment, Connectivity, ActualResponseMean, ActualResponseSD,
# BurUtilMean, BurUtilSD
# $PerformanceRecord: a matrix of (repl_ct, iter, response that iteration for each bureaucrat)
# $ObstyRecord: a matrix of (repl_ct, obsty matrix for each trial)
# $SeenRecord: a matrix of (repl_ct, iter, vector of whether each bureaucrat was seen that iter and trial)
# - e.g.,
# wss <- BigLoop1a(1000)
# conducts a simulation with default parameters at 1000 Replications (replications)
# Defined in ImitSim
# Replications - Number of Replications (or replications of an original simulation)
# ResponseParms - Response Parameters (mean, sd ~ unif)
# SupObsParms - Observability of each bureaucrat (mean, sd ~ unif)
# BurObsParms - Observability of each bureaucrat to one another (mean, sd ~ unif) #Needed???
# Obsty - Symmetric (NumBurs x NumBurs) matrix of probabilities of observability ~ unif
# SupObsty - Vector of probabilities that supervisor sees each bureaucrat ~ normal #? why normal??
# PrefParms - Preferences of the policy (~unif(-1,1))
# Tolerance - # of Std Deviations from others that supervisor permits ~ unif(-2,0)
# Std - Minimum acceptable performance (~unif)
# Punishment - How much punishment is levied (unif(-2,0))
# Performance - Record of each bureaucrat's performance (MaxIter x NumBurs)
# supervision - "Relative" (1a) or Fixed" (1b/1c/2b/2c)
# posprefs - FALSE (default) or TRUE (prefs > 0 (2a/2b/2c))
# omniscient - FALSE (default) or TRUE (supervisor can see everyone, 1c, 2c)
# quiet - FALSE (default) or TRUE (no output)
#' ImitSim
#'
#' @param Replications integer
#' @param NumBurs integer
#' @param MaxIter integer
#' @param supervision character
#' @param posprefs logical
#' @param omniscient logical
#' @param Tolerance double
#' @param Std double
#' @param Punishment double
#' @param ... additional parameters
##' @param SupObsParms vector(2)
##' @param ResponseParms vector(2)
##' @param PrefParms vector(2)
##' @param BurObsParms vector(2)
##' @param quiet logical
##' @param debug logical
#'
#' @return simulation object
#' @export
#' @import tidyverse
#'
ImitSim <- function(Replications,
NumBurs=10,
MaxIter=10,
supervision="Relative",
posprefs=FALSE,
omniscient=FALSE,
Tolerance=-99,
Std=-99,
Punishment=-99,
SupObsParms=c(-99,-99),
PrefParms=c(-99,-99),
ResponseParms=c(-99,-99),
BurObsParms=c(-99,-99),
quiet=FALSE,
debug=FALSE) {
# This routine randomly draws the assorted parameters
# and stores the final mean supoutc as a dependent var
# Version 1a calls policy 1a (relative punishment) based on punrate
if (!quiet) {
cat("Running ImitSim over", Replications, "Replications\n")
cat(" with NumBurs", NumBurs, "\n")
cat(" over MaxIter", MaxIter, "\n")
cat("Other parameters:\n")
cat(" supervision:", supervision, "\n")
if (!posprefs) cat(" with random preferences\n")
else cat(" with positive preferences\n")
if (!omniscient) cat(" without omniscient supervisor\n")
else cat(" with omniscient supervisor\n")
}
if (SupObsParms[1] != -99 & !quiet) cat("SupObsParms", SupObsParms, "\n")
Init_SupObsParms <- SupObsParms
if (Tolerance != -99 & !quiet) cat("Tolerance", Tolerance, "\n")
Init_Tolerance <- Tolerance
if (Punishment != -99 & !quiet) cat("Punishment", Punishment, "\n")
Init_Punishment <- Punishment
if (Std != -99 & !quiet) cat("Std", Std, "\n")
Init_Std <- Std
if (ResponseParms[1] != -99 & !quiet) cat("ResponseParms", ResponseParms, "\n")
Init_ResponseParms <- ResponseParms
if (PrefParms[1] != -99 & !quiet) cat("PrefParms", PrefParms, "\n")
Init_PrefParms <- PrefParms
if (BurObsParms[1] != -99 & !quiet) cat("BurObsParms", BurObsParms, "\n")
Init_BurObsParms <- BurObsParms
if (debug) cat("DEBUGGING ACTIVE\n")
ExecSummary <- array(data=rep(Replications*17,0), dim=c(Replications, 17))
if (supervision == "Relative")
colnames(ExecSummary) <- c("SupUtilMean", "SupUtilSD", "PrefMean", "PrefSD", "RespMean", "RespSD",
"SupObsMean", "SupObsSD", "BurObsMean", "BurObsSD",
"Tolerance", "Punishment", "Connectivity", "ActualResponseMean", "ActualResponseSD",
"BurUtilMean", "BurUtilSD")
else if (supervision == "Fixed")
colnames(ExecSummary) <- c("SupUtilMean", "SupUtilSD", "PrefMean", "PrefSD", "RespMean", "RespSD",
"SupObsMean", "SupObsSD", "BurObsMean", "BurObsSD",
"Std", "Punishment", "Connectivity", "ActualResponseMean", "ActualResponseSD",
"BurUtilMean", "BurUtilSD")
# Store the results overall in Performance ## Do I want Performance Record to be global?
Performance <- array(NA, dim=c(Replications*MaxIter,2+NumBurs))
ObstyRecord <- array(NA, dim=c(0, 1+NumBurs))
SeenRecord <- array(NA, dim=c(0, 2+NumBurs))
# Loop starts here ("go")
for (repl_ct in 1:Replications) {
if (ResponseParms[1] == -99) {
ResponseParms <- runif(2)
ResponseParms[1] <- -1+2*ResponseParms[1]
}
# draw SupObsParms randomly for policy 1a/1b/2a/2b
if (SupObsParms[1] == -99) SupObsParms <- runif(2)
# uncomment for policy 1c/2c
if (BurObsParms[1] == -99) BurObsParms <- runif(2)
### Probably a bug from the original code: the Obsty matrix was always hyperconnected,
### which really shouldn't be the case.
###
### observability matrix (for bureaucrats) dist'd unif(0,1)
### symmetric (x sees y sees x equally), where each entry is
### probability of seeing
### (obsty has to be global for the connectivity check)
###Obsty <- array(runif(NumBurs^2),dim=c(NumBurs,NumBurs)) ## seems like numburs must be global??
# observability matrix (for bureaucrats) dist'd normal(mean=BurObsParms[1], sd=BurObsBarms[2])
# but then converted to PROBABILITIES
# symmetric (x sees y sees x equally), where each entry is
# probability of seeing
# (obsty has to be global for the connectivity check)
Obsty <- array(rnorm(NumBurs^2, mean=BurObsParms[1], sd=BurObsParms[2]),dim=c(NumBurs,NumBurs))
Obsty <- pnorm(Obsty)
Obsty[lower.tri(Obsty)] <- t(Obsty)[lower.tri(Obsty)]
Obsty <- Obsty>.5
Obsty[Obsty==0] <- NA
diag(Obsty) <- 1
ObstyRecord <- rbind(ObstyRecord,cbind(rep(repl_ct,NumBurs),Obsty))
# Add to ObstyRecord
#ObstyRecord[(repl_ct-1):(repl_ct-1)+NumBurs, ] <- Obsty
# supervisor observability (see and be seen, here) dist'd unif(0,1) ## I say unif in the code, but seems norm??
SupObsty <- rnorm(NumBurs, mean=SupObsParms[1], sd=SupObsParms[2])
SupObsty[SupObsty < 0] <- 0
SupObsty[SupObsty > 1] <- 1
# Set the range of Preferences (why did I call it popparms??) ## NOT IN ORIGINAL CODE
if (posprefs == FALSE & Init_PrefParms[1] == -99) PrefParms <- c(runif(1, min=-1, max=1), runif(1))
else if (posprefs == TRUE & Init_PrefParms[1] == -99) PrefParms <- c(runif(1, min=0, max=1), runif(1))
### I'm not sure on the is.null call here. I *know* it's not null
# "Relative" is for 1a/2a; "Fixed" is for 1b/1c/2b/2c
if (supervision == "Relative" & Init_Tolerance != -99) Tolerance <- runif(1, min = -2, max=2)
else if (supervision == "Fixed") Std <- runif(1, min=-1, max=1)
if (Init_Punishment == -99) Punishment <- runif(1,min=0, max=2)
Connectivity <- igraph::vertex_connectivity(igraph::graph_from_adjacency_matrix(Obsty))
# this version of the policy routine
# assigns outcomes on the basis of dowhat * response
# unless the saboteur is caught
# Response - What Bureaucrats do this round (~ truncated normal ([-1,1], resp_parms))
# Catch - Do I catch anyone? (dummy)
# Prefs - Preferences of each Bureaucrat (~ trunc normal ([-1,1], prefs_parms))
# BurUtil - Bureaucrats' Utility
# SupUtil - Supervisor's Utility (mean and sd of bureaucrats' performance)
# SupSeen - Who does supervisor see this round?
# BurUtilSeen - Utility of Bureaucrats who are Seen by Bureaucrats
# Envy - Which bureaucrat does each bureaucrat envy (highest utility)?
# iter - iteration counter
#
# get initial conditons
# burs respond over -1 to 1, dist'd normally
Response <- rnorm(NumBurs, mean=ResponseParms[1], sd=ResponseParms[2])
Response[Response < -1] <- -1
Response[Response > 1] <- 1
# change call below to switch policy variation
# set counters to zero, open result matrices
Catch <- 0
BurUtil <- array(0, dim=c(MaxIter, NumBurs)) ## Do I want this to be global?
SupUtil <- array(0, dim=c(MaxIter, 2))
#getpolicy(mnpop) ## divergence!!! scope prevents this from working
Prefs <- rnorm(NumBurs, mean=PrefParms[1], sd=PrefParms[2]) ## whoa! was I allowing Prefs outside of -1,1???
if (posprefs == TRUE) {
Prefs[Prefs < 0] = 0
Prefs[Prefs > 1] = 1
}
if (debug) {
cat("\n============\nReplication:", repl_ct, "\n============\n")
cat("Obsty:\n")
print(Obsty)
}
# play: go from here
for (iter in 1:MaxIter) {
# Do ranges from -1 (complete defection) to 1 (complete compliance)
Do <- Response
# Record Performance
Performance[(repl_ct-1)*MaxIter+iter, 1:2] <- c(repl_ct, iter)
Performance[(repl_ct-1)*MaxIter+iter, 3:(NumBurs+2)] <- Response
# outcome equals desires times do
if (debug) {
cat("Prefs:\n")
print(Prefs)
cat("Do:\n")
print(Do)
}
BurUtil[iter,] <- Prefs*Do
if (debug) {
cat("BurUtil (before supervision)\n")
print(BurUtil[iter,])
}
# supervisor enforces against those she sees defecting
SupUtil[iter,] <- c(mean(Do), sd(Do))
Seen <- trunc(SupObsty+runif(1))
if (debug) {
cat("Seen:\n")
print(Seen)
}
SeenRecord <- rbind(SeenRecord, cbind(repl_ct, iter, t(Seen)))
# Relative tolerance => deviants are those whose behavior is below "Mean - Tolerance*SD"
# Fixed supervision => deviants are less than Std units (on a -1,1 scale)
if (supervision == "Relative") deviants <- Do < (SupUtil[iter,1] - (Tolerance*SupUtil[iter,2]))
else if (supervision == "Fixed") deviants <- Do < Std
if (debug) {
cat("deviants:\n")
print(deviants)
}
# ERROR???still have to be deviant even if omniscient
# FIX? use "logical()" above to test for deviants
if (!omniscient) seendeviants <- (deviants * as.integer(Seen))>0
else seendeviants <- deviants>0
if (debug) {
cat("seendeviants:\n")
print(seendeviants)
}
WhoCaught <- which(as.logical(seendeviants))
if (debug) {
cat("WhoCaught:\n")
print(WhoCaught)
}
if(length(WhoCaught)!=0) {
Catch <- Catch + 1
BurUtil[iter,WhoCaught] <- BurUtil[iter,WhoCaught] - Punishment
}
if (debug) {
cat("---\nBurUtil (after supervision)\n")
print(BurUtil[iter,])
}
# adapt:
BurUtilAll <- array(data=BurUtil[iter,], dim=c(NumBurs, NumBurs))
if (debug) {
cat("BurUtilAll\n:")
print(BurUtilAll)
}
BurUtilSeen <- array(NA, dim=c(NumBurs, NumBurs))
for (i in 1:NumBurs) BurUtilSeen[i,] <- BurUtilAll[i,] * Obsty[i,]
if (debug) {
cat("\nBurUtilSeen:\n")
print(BurUtilSeen)
}
# choose response of the person you saw who did best
Envy <- apply(BurUtilSeen, 2, which.max)
if (debug) {
cat("\nEnvy:\n")
print(Envy)
}
if (debug) cat("\nOld Response:\n", Do)
Response <- Do[Envy]
if (debug) cat("\nNew Response:\n", Response)
# play again, unless this iteration exceeds maxiter
}
ExecSummary[repl_ct, 1:10] <- c(SupUtil[MaxIter,], PrefParms, ResponseParms, SupObsParms, BurObsParms)
# will need to change punrate in line below to std for 1b/1c/2b/2c
if (supervision == "Relative")
ExecSummary[repl_ct, 11:15] <- c(Tolerance, Punishment, Connectivity, mean(Response), sd(Response))
else if (supervision == "Fixed")
ExecSummary[repl_ct, 11:15] <- c(Std, Punishment, Connectivity, mean(Response), sd(Response))
ExecSummary[repl_ct, 16:17] <- c(mean(BurUtil[MaxIter,]), sd(BurUtil[MaxIter,]))
if (Init_SupObsParms[1] == -99) SupObsParms <- c(-99,-99)
if (Init_Tolerance == -99) Tolerance <- Init_Tolerance
if (Init_Std == -99) Std <- -99
if (Init_Punishment == -99) Punishment <- -99
if (Init_PrefParms[1] == -99) PrefParms <- c(-99, -99)
if (Init_ResponseParms[1] == -99) ResponseParms <- c(-99, -99)
if (Init_BurObsParms[1] == -99) BurObsParms <- c(-99,-99)
}
BigList <- list("Version" = "2.02",
"NumBurs" = NumBurs,
"MaxIter" = MaxIter,
"Replications" = Replications,
"supervision" = supervision,
"posprefs" = posprefs,
"omniscient" = omniscient,
"ExecSummary" = as.data.frame(ExecSummary),
"Performance" = Performance,
"ObstyRecord" = ObstyRecord,
"SeenRecord" = SeenRecord)
BigList
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.