A Dynamic Simulation of a Zombie Apocalypse

# Zombies vs Humans Agent Based Simulation (the r script file in case blogger mangled my code)

# I also wrote a Spatial Simulation of a zombie infection in Stata previously.

# This simulation is a simple repeated matching simulation in which each agent is matched with a random different agent.

# I believe this is an appropriate way of modelling a human zombie exchange as portrayed in the movies.  Generally speaking, each encounter can be thought of as a probabilistic draw in which either the human becomes zombified or the zombie is permanently killed.

# Agents all start out as humans (except a random percent which are initially zombified).  Each human is

defined on a scale of 0 to 90 percentile with increments of 10 which reflect the chance of that human vanquishing a zombie that the human encountered or being in turn vanquished.

# If a human is vanquished a zombie is added to the zombie population.

# number of humans (approximated due to rounding issues)
start.humans = 100000

# matrix of human types
htypes = seq(0,90,10)

# Frequency of each time of human from 0 to 90 percentile.
# This number is only relative to the scale of the other frequencies so long as it is positive.
# Thus if only type of human had a 5 then it would be 5 times more likely than a type of human with a frequency of 1.
freq = rep(1,length(htypes))

# frequency is the most important parameter choice in the model as will be seen below.

# Bind the information into a single matrix.
human.types = cbind(htypes,freq)

# Now we calculate what percentage of our start.humans are each type.
perc = round((freq/sum(freq))*start.humans)

# Finally we generate our moving things data.
# Initially the only thing moving is humans.
walking.things = rep(htypes,perc)

# Looking good
walking.things

# Some percentage of the humans are initially infected.
infected.per = .025

# Calculate the  number initially that become infected.
nselected = round(start.humans*infected.per)

# Now we randomly select the humans infected initially.
initial.zombies = sample(1:length(walking.things), nselected)

walking.things=c(walking.things[-initial.zombies], rep(-77,nselected))
           
# -77 Is the number for a zombie.
   
# After the intial infection phase peopole get their guns out and start acting defensively.
walking.things

# Percent zombies
perc.zombies = mean(walking.things==-77)

# Total population (living and dead)
nthings.vector = nthings = length(walking.things)

# Count the number of zombies
nzombies = sum(walking.things==-77)

# Count the number of humans
nhumans = sum(walking.things!=-77)

  nhumans0  = sum(walking.things==00)
  nhumans10 = sum(walking.things==10)
  nhumans20 = sum(walking.things==20)
  nhumans30 = sum(walking.things==30)
  nhumans40 = sum(walking.things==40)
  nhumans50 = sum(walking.things==50)
  nhumans60 = sum(walking.things==60)
  nhumans70 = sum(walking.things==70)
  nhumans80 = sum(walking.things==80)
  nhumans90 = sum(walking.things==90)


# This command pairs up a vector.  It is used to match humans with zombies.
pairup = function(x, unmatched.self=T) {

  # Calculate the length of the x vector.
  xleng = length(x)

  # This checks if the input vector is a scalar.
  if (xleng==1) x = 1:(xleng=x)

  # Half the length of x rounded down.
  hleng = floor(xleng/2)

  # Randomize x
  x = x[order(runif(xleng))]
  pairs = cbind(x[1:hleng],x[(hleng+1):(2*hleng)])

  # If there is a odd number of xs then this will match the remaining unmatched x with itself if unmatched.self is T.
  if ((unmatched.self)&(xleng/2!=hleng)) pairs=rbind(pairs, c(x[2*hleng+1]))

  return (pairs)
}

#
max.rounds = 45

# Let's start the simulation:
n = 1
while  (nzombies[n]>0 & nhumans[n]>0 & n  n = n+1

  # This calls the previously defined function pairup to match two different individuals together.
  # This matches them by position in the walking.things vector.
  encounter=pairup(nthings)

  # This assigns to the matrix the values
  types = cbind(walking.things[encounter[,1]],walking.things[encounter[,2]])

  # Create a vector of terminated or zombified things
  conflict = types*0
     # 0 Unresolved
     # 1 Zombified
     # 2 Permenent Death
     # 3 No conflict
     # 4 win conflict
 
   # This code will check if a zombie is in the right column and human in the left.
   hvz = (types[,2]==-77)&(types[,1]>=0)
     # If so, the human and zombie places will be switched.
     types.temp = types
     types[hvz,1]=types.temp[hvz,2]
     types[hvz,2]=types.temp[hvz,1]
   
     encounter.temp = encounter
     encounter[hvz,1]=encounter.temp[hvz,2]
     encounter[hvz,2]=encounter.temp[hvz,1]
   
   # Zombie encounters human
   zvh = (types[,1]==-77)&(types[,2]>=0)
     # Calculate the win count of the conflict
     win.zvh = (runif(sum(zvh))>types[zvh,2]/100)
   
     # Translate a zombie win onto the conflict map
     conflict[zvh,1][win.zvh]=4
     conflict[zvh,2][win.zvh]=1  
   
     # Translate a human win onto the conflict map
     conflict[zvh,1][!win.zvh]=2
     conflict[zvh,2][!win.zvh]=4

     # Resolve non-conflict. Zombies don't fight zombies and humans don't fight humans.
     conflict[types[,1]==types[,2],] = 3
     conflict[(types[,1]>=0)&(types[,2]>=0),] = 3
   
   # Finally, adjust the walking.things vector to adjust for the changes.
   # Zombify some
   walking.things[encounter[conflict==1]] = -77
   
   # Remove others
   walking.things=walking.things[-encounter[conflict==2]]

  # Store stats
  # Percent zombies
  perc.zombies = c(perc.zombies, mean(walking.things==-77))

  # Total population (living and dead)
  nthings = length(walking.things)
  nthings.vector = c(nthings.vector, nthings)

  # Count the number of zombies
  nzombies = c(nzombies, sum(walking.things==-77))

  # Count the number of humans and save them to vectors
  nhumans = c(nhumans,  sum(walking.things!=-77))

  nhumans0  = c(nhumans0,   sum(walking.things==0))
  nhumans10 = c(nhumans10,  sum(walking.things==10))
  nhumans20 = c(nhumans20,  sum(walking.things==20))
  nhumans30 = c(nhumans30,  sum(walking.things==30))
  nhumans40 = c(nhumans40,  sum(walking.things==40))
  nhumans50 = c(nhumans50,  sum(walking.things==50))
  nhumans60 = c(nhumans60,  sum(walking.things==60))
  nhumans70 = c(nhumans70,  sum(walking.things==70))
  nhumans80 = c(nhumans80,  sum(walking.things==80))
  nhumans90 = c(nhumans90,  sum(walking.things==90))

}

# Count the number of rounds completed.
nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: The number at the end of each line is the probability that an individual\nin this population group will kill a zombie when encountering one. Z is the zombie population"
     , main="Population During a Zombie Attack\nWith a well armed human population a zombie apocalypse is easily prevented")
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# However, this equal proportion of highely effective zombie killers to very ineffective zombie killers is nonrepresentational of a typical zombie movies or games.

# I will instead rerun the simulation with a larger percentage of low ability humans.
htypes = seq(0,90,10)
freq = 10:1

# .... using same code as above but with new human population proportions

nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: Even the best trained individual will be overwhelmed eventually\nif there are too many easy zombie victims. Z is the zombie population"
     , main="Population During a Zombie Attack\nA poorly armed population is ill-equiped to survive a zombie attack")
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# Let's try one more variant with still more weak humans but a little better ratios.

htypes = seq(0,90,10)
freq = seq(5,1, length.out=10)

# .... using same code as above but with new human population proportions

nrounds.completed = length(nhumans0)

plot(c(1,nrounds.completed), c(0,max(nhumans0,nzombies)), type="n",
     ylab="Population", xlab="Round\n*Note: In this scenario only the top 10 to 20% most effective zombie killers survive."
     , main="Population During a Zombie Attack-When the population of dangerous humans\n is sufficiently large, it is possible humanity survives, just barely")
   
for (i in seq(0,90,10)) {
  lines(get(paste("nhumans",i,sep="")))
  text(nrounds.completed+1, get(paste("nhumans",i,sep=""))[nrounds.completed], i)
}
lines(nzombies, lwd=3)
  text(nrounds.completed+1, nzombies[nrounds.completed], "Z")

# Overall conclusion? Mandating zombie defense classes is the only way to be certain humanity will survive.

Concerto v4.b web-based Content Management System for R Users

Concerto is a open source R Content Management System for developing online applications.  The primary purpose of Concerto is for the designing of online tests.  However, the hard working programmers at the University of Cambridge Psychometrics Centre have created a program that is flexible enough to address a wide range of applications such as the development of statistical learning tools, client management systems, as well as large scale survey deployment.  Future plans for Concerto would include the integration of it and the Psychometric's Center's facebook App which has collected detailed information for millions of respondents spanning a number of fascinating surveys.

The real beauty of Concerto is that it provides an interface that builds on the vast resources of R in order to seamlessly integrate a powerful environment capable of managing the enormous contributions of web developers.

Well, enough advertising.  Let's see how v4b works.

In this post I will map out the logic of Concerto environment as I understand it now.  I will flesh out more in future posts as I get a better grasp of the Concerto environment.

For this post and all future posts I will look only at version 4b of Concerto or higher.

For the test developer the base operating system is R.  For now I will use pseudo R-like code.  Pseudo code is useful at explaining coding concepts.

First we start with an initialization of a “test”.  Tests are made up of R code.

# START – pseudo code random item survey
# The first thing a test might do is have a user input demographic information.
1.> intro = Concerto.html.tempate(Introduction)

# The responses from the Introduction page are returned as a list assigned to the Intro object.
# This responses might include any kind of information like: age, name, gender, ect.  We can target the values the same way as any list object.
2> print(intro$age)

# Interestingly this command does not actually do anything as far as the user is concerned because the screen output from R is suppressed.  However in debugging mode screen returns from R are displayed in the testing window.
# We may write the information we have gathered so far to our user information table which is a sql table. Concerto R is able to send commands to sql.  Remember this is only pseudo-code.
3> Concerto.sql.send(user_table, user_name=intro$user.name, user_age=intro$user.age)

# Now let’s administer our test items.  Let’s say we have all of our items in a table.
4> items = concerto.sql.get(items)

# We can figure out how many we are by counting the rows.
5> nitems = nrows(items)

# We want to select 20 items randomly so that almost nobody shares item sets.
6> item.choice = sample(1:nitems, 20)

# Now to administer the items we simply send them to the user by looping through all of them.
7> response.vector = NULL
8> for (i in item.choice) response=c(response, Concerto.html.tempate(item.admin, params=list(text=items$text[i]))

# Now we have two vectors.  A vector of item responses and a vector of items administered.
# All that is left is to save this response pattern
9> Concerto.sql.send(response_table, items=item.choice, responses=response.vector)

# And provide the user with some feedback
10> feedback = feedback.analysis(item.choice, response.vector)

# Now we simply display the results of this feeback analysis (this of course would be a function that we have already coded.
11> concerto.html.template(Feedback)

# END – pseudo code random item survey

In future posts I will strive to post functioning code and functioning tests.  To learn more about Concerto check out the development website:
https://code.google.com/p/concerto-platform/

Block separated commenting or not?

I recently received the following anonymous comment on my most recent post:

"Your blog is hard to read due to the difficulty in separating text from r code. Why not separate the two instead of comments throughout the code?"

I have been writing a response and decided to post it as an independent post since many people may have similar concerns who read the blog.

Let's first define what I think of as the alternative style that most blogs use.  Please don't take offense all of you fantastic bloggers out there.  I think this format as with all blog formats is extremely useful.  I just think that there are reasons why other formats might also be desirable.

------------------------------------------------------------------------------------
Typical posts look this way:

Explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation.

Explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation.

Explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation.

code block:
> some useful code
> some useful code
> some useful code
> some useful code
> some useful code

Explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation explanation.

----------------------------------------------------------------------------------

My Reasons for doing things the way I do them: such as in Tit for Tat - Axelrod tournament style competitive simulation

1. I believe this is a matter of style.  Personally I find blogs in which there are paragraphs of comments followed by small amounts of code difficult to read.  I would like to figure out how to get my syntax color coded.  However, I have not figured out a way to do this which does not result in code boxes and the like.

2. The format of the blog makes it very easy for someone to take an entire post and copy it right into a code editor.  From that point, they can run the code piece by piece.  This helps to address the lack of syntax color coding.  In whatever favorite editor people may have they can see my code color coded.

3. I was uncertain about this style initially as well. However, many more people have visited my blog than I had originally ever hoped.  Upwards of 100,000 page views at this point.  Why mess with what seems to work?

4. There is always a trade-off between quantity and quality.  Some might be bothered by this but I have aimed to put out more quantity than quality.  I don't mean to degrade my work but often there are spelling errors in posts and other minor errors that might bother people.  However, this is not my job and so I must squeeze in posts when I have the time and inclination.  I experimented initially with separating code from comments but found the structure much more time consuming.

5. I post things to my blog not to give people free code that they can use for their own purposes but to show people how to code various things.  From what I have seen, when code is fully commented ex-ante then it is either very small amounts of code such as a 1 to 5 line command or function.  Or it is using a very large amount of code that users are not expected to read but only copy and paste.  I am not interested in coding solutions to problems for others but rather in showing people how to become competent coders.

6. I think separating the code body from the comments inherently limits the ability of a blogger to post complex or lengthy code.  Some have suggested that they don't want complex or lengthy code in posts. That is fine but there are many other blogs out there which already fulfill these requirements and I don't see a reason to change what I am doing.

Finally, I would like to say that I have no problem with the way other people structure their blogs.  It works for them and this works for me.  Having a diversity of options available to satisfy the tastes of many different people is a great thing.

Meta-Analysis Simulation - Monte Carlo

# Meta-Analysis is a method that allows researchers to merge the scientific results from a number of studies into a single analysis.  This is an especially useful method because it allows for a collective estimate of an effect using a pool of data much large than could be achieved with any single study.
n.studies = 50

# This stores the simulation results from all runs varying the true effect.
total.sim.results=matrix(NA, nrow=0, ncol=5)

# This is true we want to detect.
for (true.eff in c(0,.05,.10,.15,.20)) {

ind.sd = 1

s.scalar = 100
# Study scale is drawn from a chi-squared distribution times the scalar + 5

# unit scale. Normal distribution.
mean.scale = 1

# sd scale. Uniform distribution.
sd.scale = 20

# In this analysis I will structure this as a Monte Carlo simulation to test how this method works.  Please note, I am just working from class notes on this so I might be formulating the meta-analysis proceedure incorrectly.
n.simulations = 100
# We will save the results of the Monte Carlo in a matrix.
sim.results = matrix(NA, nrow=n.simulations, ncol=5)

print(paste("True Effect:",true.eff))

for (s in 1:n.simulations) {


# First let's generate our study parameters:

s.size = ceiling(rnorm(n.studies)^2 * s.scalar)+15
s.sd = runif(n.studies)*sd.scale
s.mean = rnorm(n.studies)*mean.scale

# For didactic purposes I will simulate the individual studies though this is not neccessary.

ind.treat = ind.data = matrix(NA,nrow=max(s.size), ncol=n.studies)

est.corr = est.y.sd = est.x.sd = est.effect = rep(NA, n.studies)

for (i in 1:n.studies) {
  ind.treat[1:s.size[i],i] = rnorm(s.size[i])>0
  ind.data[1:s.size[i],i] = (true.eff*ind.treat[1:s.size[i],i] + rnorm(s.size[i])*ind.sd)*s.sd[i] + s.mean[i]

  # Estimate effect for each
  est.effect[i] = lm(ind.data[1:s.size[i],i]~ind.treat[1:s.size[i],i])$coef[2]

  # Estimate the standard deviation for each study
  est.y.sd[i] = sd(ind.data[1:s.size[i],i], na.rm=T)
  est.x.sd[i] = sd(ind.treat[1:s.size[i],i], na.rm=T)

  # Estimate correlation between treatment and outcome
  est.corr[i] = cor(ind.data[1:s.size[i],i], ind.treat[1:s.size[i],i])
}
# Because everybody is using a different scale, the estimated affect is unique to that study's scale.
est.effect

# We therefore divide the estimated effect by the
z.est.effect = est.effect/est.y.sd

z.est.effect

# We can convert our standardized estimated effects to correlations:
z.est.effect*est.x.sd

est.corr

# In order to aid analysis, tranform the data to fisher z
z.corr = .5*log((1+est.corr)/(1-est.corr))

mean(z.corr)

# The variance of the inidividual Fisher zs can be approximated with:
z.var = 1/(s.size-3)

# Calculate the between study variance
between.var = var(z.corr)

# In order to construct the test we need to specify weights which are equal to the inverse of the variance the study plus the variance between studies.
weight = 1/(between.var+z.var)

# Now we can estimate to see if our effect is statistically signficant.
  (sim.results[s,]=c(true.eff,summary(lm(z.corr~1, weights=weight))$coef))
# Name the columns of the results matrix.
  if (s==1) colnames(sim.results)=c("true.eff" ,colnames(summary(lm(z.corr~1, weights=weight))$coef))
# Provide feedback to the user.
  print(paste("Number of repetitions completed:",s))
}
total.sim.results = rbind(total.sim.results,sim.results)
}

# We can use the Monte Carlos results to test several features:

# 1. Test if our estimator is overpowered.
# To do this we set the true effect equal to zero and see how frequently we reject at the 10% and 5% level (true.eff=0)
mean(total.sim.results[total.sim.results[,1]==0,5]<.1)
# If this number is substantially higher than 10% or 5% then we should worry about our estimator giving us results when there actually are no results.

# 2. We can test the power of our analysis to reject the null when there is an effect of a particular size (for example: true.eff=.1)
mean(total.sim.results[total.sim.results[,1]==.1,5]<.1)
# true.eff=.2
mean(total.sim.results[total.sim.results[,1]==.2,5]<.1)

boxplot(total.sim.results[,5]~total.sim.results[,1], main="P values", ylab="P value", xlab="True Effect")
abline(h=.1, col=grey(.75))

The effect of non-convergence on MLE estimates

* Maximum likelihood proceedures have become widely used to solve a variety of econometric problems.

* Unfortunately there is no guarantee that these proceedures will yeild a single solution which satisfies the convergence criteria of the maximizing function.

* This might occur for reasons difficult to detect such as localy flat spots or discontinuous areas.

* Maximization proceedures are usually evaluated based on their 1. efficiency (speed of convergence) or 2. on their robustness at detecting optimal values.

* The problem is that sometimes in simulation we need to limit the time a MLE proceedure takes in attempting to find a solution.

* What effect does that limitation result in and what do we do with estimates that result from non-convergence?
* 1. Keep them or 2. throw them out

* This simulation will explore both options

* In this simulation we will fall back on the widely used estimator which is equivalent when the standard errors are not structurally estimated to the OLS estimator.

* This is the "normal" regression estimator.  IE the MLE maximization that allows for linearly modeled heteroskedasticity.
cap program drop myNormalReg
program define myNormalReg
  args lnlk xb sigma2
  qui replace `lnlk' = -ln(sqrt(`sigma2'*2*_pi)) - ($ML_y-`xb')^2/(2*`sigma2')
end

* First let's generate a sample data set
clear
set obs 300

* I am going to try to make the problem hard to solve by including both addative and multiplicative error.
gen u = (runiform()-.5)
  * I made this error small because actually when the error is small it is harder to estimate the variance of the error.
  * It takes a little work with simulations to generate data which does not converge.

gen v1 = rnormal()
gen x1 = runiform()-.5

gen v2 = rnormal()
gen x2 = runiform()-.5

gen y = 3 + (1+v1)*x1 + (1+v2)*x2 + u

reg y x1 x2

ml model lf myNormalReg (reg: y=x1 x2) (sigma2:)
ml maximize

* This is the more efficient model because it is explicitly modeling the error.
gen x1_2 = x1^2
gen x2_2 = x2^2

ml model lf myNormalReg (reg: y=x1 x2) (sigma2: x1_2 x2_2)
ml maximize

* It seems that typically this model frequently converges.

* Let's see if we can't dilute the maximization:
gen x1x2 = x1*x2

gen x1abs = abs(x1)
gen x2abs = abs(x2)

ml model lf myNormalReg (reg: y=x1 x2 x1_2 x2_2 x1x2 x1abs x2abs) (sigma2: x1 x2 x1_2 x2_2 x1x2 x1abs x2abs)
ml maximize, iterate(100)

* It looks to me about half the time I run this code it does not converge within a 100 iterations.

* Now lets specify our functions we will use to test differences in results depending upon our method of dealing with convergence.

* This program is just a condensation of the above code.
cap program drop sim_converge
program define sim_converge

  clear
  set obs 300

  gen u = (runiform()-.5)

  gen v1 = rnormal()
  gen x1 = runiform()-.5

  gen v2 = rnormal()
  gen x2 = runiform()-.5
  gen y = 3 + (1+v1)*x1 + (1+v2)*x2 + u
  gen x1_2 = x1^2
  gen x2_2 = x2^2

  ml model lf myNormalReg (reg: y=x1 x2) (sigma2: x1_2 x2_2)
  ml maximize

  gen x1x2 = x1*x2

  gen x1abs = abs(x1)
  gen x2abs = abs(x2)

  ml model lf myNormalReg (reg: y=x1 x2 x1_2 x2_2 x1x2 x1abs x2abs) (sigma2: x1 x2 x1_2 x2_2 x1x2 x1abs x2abs)
  ml maximize, iterate(`1')
  * The only difference is that iterate is specified by the user.
end

* Leaving the first argument blank will not specify a maximum convergence iteration.
sim_converge
sim_converge 50

* Let's first define what we would like to save from the MLE.
* Yes, I am going to use a forbidden global :)
gl savelist ic=e(ic)
  * e(ic) is the macro in which the number of iterations used is saved.

 foreach i in reg sigma2 {
   foreach v in x1 x2 x1_2 x2_2 x1x2 x1abs x2abs {
   gl savelist $savelist `i'`v'=[`i']_b[`v']
   }
 }

 * Let's see what our savelist looks like:
 di "${savelist}"

 * looking pretty good.

simulate ${savelist} , rep(100) seed(32): sim_converge 50
tab ic
* In my simulation 51 times the MLE did not converge by the 50th iteration.

* Let's see if there are systematic differences between estimates.
sum if ic==50
sum if ic<50 p="">  * We can see that if the estimator did converge then it is much more precise (smaller sd) than in the cases when it did not converge.
 
  * The mean estimates of regx1 and regx2 and sigmax1_2 and sigmax2_2 are much closer to 1 which is the true parameter values.
 
* Let's try it again setting convergence at a higher bar:
simulate ${savelist} , rep(100) seed(32): sim_converge 250
tab ic
sum if ic==250
sum if ic<250 p="">
* Raising the max iteration does not lead to any of the observations converging.

* This is problematic because we want to know if there is a systematic difference in the draws for the estimates which converged and those that did not.

* By the results so far we might be tempted just to include the results of the iterations that did converge.

* First off let's see if the estimates that converged quickly are better or worse than those that converged more slowly.

recode ic (1/14=0) (15/49=1), gen(grp)

bysort grp: sum regx1 regx2
anova regx1 grp if grp<250 p="">* It seems there is no detectable difference in the means for those observations that converged more quickly than 15 iterations than those that converged more slowly.

* This implies, assuming the results are generalizable that truncating the simulation to only the results that converge might produce unbiased estimates.

* We should run the simulation again with more repetitions in order to confirm this.

simulate ${savelist} , rep(500) seed(32): sim_converge 50
  tab ic
 
/* tab ic

      e(ic) |      Freq.     Percent        Cum.
------------+-----------------------------------
         10 |          1        0.20        0.20
         11 |         18        3.63        3.83
         12 |         24        4.84        8.67
         13 |         48        9.68       18.35
         14 |         57       11.49       29.84
         15 |         50       10.08       39.92
         16 |         24        4.84       44.76
         17 |         10        2.02       46.77
         18 |         18        3.63       50.40
         19 |         10        2.02       52.42
         20 |          4        0.81       53.23
         21 |          6        1.21       54.44
         22 |          1        0.20       54.64
         23 |          3        0.60       55.24
         24 |          4        0.81       56.05
         25 |          2        0.40       56.45
         26 |          3        0.60       57.06
         29 |          1        0.20       57.26
         32 |          1        0.20       57.46
         33 |          1        0.20       57.66
         50 |        210       42.34      100.00
------------+-----------------------------------
      Total |        496      100.00
*/

  sum if ic==50
  sum if ic<50 p="">
  recode ic (1/14=0) (15/49=1), gen(grp)

  bysort grp: sum regx1 regx2

/*
-> grp = 0

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       regx1 |       148    .9918629    .1076053    .732878   1.296374
       regx2 |       148    .9807507     .112081    .684716   1.447506

-----------------------------------------------------------------------------------------
-> grp = 1

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       regx1 |       138    .9967074    .1183514   .7492483   1.318624
       regx2 |       138    .9913069    .1161812   .6711526   1.407439

-----------------------------------------------------------------------------------------
-> grp = 50

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
       regx1 |       210     .814391    .3268973   .0518116   1.590542
       regx2 |       210    .8106436    .3278075   .0707162   1.775239

*/
  anova regx1 grp if grp<50 p="">
/*
                         Number of obs =     286     R-squared     =  0.0005
                           Root MSE      = .112917     Adj R-squared = -0.0031

                  Source |  Partial SS    df       MS           F     Prob > F
              -----------+----------------------------------------------------
                   Model |  .001675983     1  .001675983       0.13     0.7172
                         |
                     grp |  .001675983     1  .001675983       0.13     0.7172
                         |
                Residual |  3.62106459   284  .012750227  
              -----------+----------------------------------------------------
                   Total |  3.62274057   285   .01271137  
*/

* We can see that even when the sample size is larger (about 140 per iteration group) there is no discernable difference between those draws that converge before the first 15 iterations and those that converge after.

* If we assume that for those draws that do not converge within 250 draws are sampling from the same probability of convergence distribution then this evidence suggests that rate of convergence is independent of the actual estimates and therefore it might be safe to exclude observations in which convergence did not occur.


* Now finally what we might interested in is seeing how our estimates change (for the values in which convergence is achieved) as we include more and more estimates by way of increasing our threshold of max iterations.

* How do we do this?

* Well, this might take a little while but we basically loop through the simulation saving the results it terms of means and standard deviations from each run.

forv i = 15(5)35 {
  simulate ${savelist} , rep(100) seed(32): sim_converge `i'
  sum regx1 if ic<`i'
    global mean_x1_`i'=r(mean)
    global var_x1_`i'=r(sd)^2
  sum regx2 if ic<`i'
    global mean_x2_`i'=r(mean)
    global var_x2_`i'=r(sd)^2
}
* By the way this is a highly redundant and inefficient method.

clear
set obs 5
gen mean_x1 = .
  label var mean_x1 "Mean of x1 estimates"
gen mean_x2 = .
  label var mean_x2 "Mean of x2 estimates"
gen var_x1 = .
  label var var_x1 "Variance of x1 estimates"
gen var_x2 = .
  label var var_x2 "Variance of x2 estimates"

gen i = .
  label var i "Max # iterations"

* Save the results as variables
forv i = 1(1)5 {
  local ii = 10+`i'*5
  replace mean_x1 = ${mean_x1_`ii'} if _n==`i'
  replace mean_x2 = ${mean_x2_`ii'} if _n==`i'
  replace var_x1  = ${var_x1_`ii'}  if _n==`i'
  replace var_x2  = ${var_x2_`ii'}  if _n==`i'
  replace i = `ii' if _n==`i'
}

two (connected mean_x1 i, msize(large) lwidth(thick)) ///
    (connected mean_x2 i, msize(large) lwidth(thick)), name(means, replace)
two (connected var_x1  i, msize(large) lwidth(thick)) ///
    (connected var_x2  i, msize(large) lwidth(thick)), name(vars,  replace)

graph combine means vars, col(1) title("Estimates are insensitive to speed of convergence")

* The take away seems to be that it is safe (at least in this simulation) to exclude from your analysis simulations that did not converge in the specified iteration count.

* This simulation also suggests that it is not ideal to include in your results MLE estimates from iterations in which no convergence was achieved.

Tit for Tat - Axelrod tournament style competitive simulation

# I have been reading the Moral Animal by Robert Wright and it has gotten me interested in putting together a simulation which is modeled after the classic Robert Axelrod's tournament style influential simulations.  So here it is!  Hope you enjoy it.

# You can also find a link to the code in full here.

# I would like to simulate multi-agent a micro-economy in which agents perform the prisoner's dimemma repeatedly.

# Games like this have been important as evidence for an evolutionary psychological explanation of social cooperation.

See how differing strategies
perform against each other

# Returns for:
  # cooperate
  coop.coop = 1

  # you cheat but opponent cooperates
  cheat.coop = 2

  # You cooperate but opponent cheats
  coop.cheat = -2

  # Both cheat
  cheat.cheat = -1

# This function turns the parameters into returns.
returns = function(you,opponent) {
  ret = 0*you
  ret[(you==1)&(opponent==1)] = coop.coop
  ret[(you==0)&(opponent==1)] = cheat.coop
  ret[(you==1)&(opponent==0)] = coop.cheat
  ret[(you==0)&(opponent==0)] = cheat.cheat
  ret
}

# For example, if you cheat but your opponent cooperates for the first two rounds
returns(c(0,0,0),c(1,1,0))

# Thus the game is defined as a zero sum game when players interact randomly with each other.  But a greater than zero sum when players are systematically cooperating and less than zero sum when players are systematically cheating.

# Number of rounds:
  n.rounds = 200

# Number of agents. This number must be even to ensure every agent is matched.
  n.agents = 20

# Basic framework.
# Each round agents meet other agents.
# Each agent has access to knowledge of the other agents previous behavior (either cheat or cooperate)

# Let's define the types of agents.
# Agents get fed input h which is the list of prior history with other agent.
# The argument they return is their action.

# The stupid agents are.
always.coop  = function(h) 1
always.cheat = function(h) 0

# Random Tod
Random.Tod = function(h) 1==rbinom(1,1,.5)

# tit4tat says that if the most recent action the opponent took was cheat then cheat otherwise cooperate.
tit4tat = function(h) {
  # When the length of h = 0 then that means there is no history.
  if (length(h)==0) T
  else rev(h)[1]==1
}

# initially agressive tit4tat
ag.tit4tat = function(h) {
  if (length(h)==0) F
  else rev(h)[1]==1
}

# opportunist tit4tat
op.tit4tat = function(h) {
  if (length(h)==0) T
  else ((n  # This is a somewhat tricky action behavior.
  # It states that play tit for tat up to the point near when the game ends.
  # Then start cheating because the game is coming to an end.
}


# fairbot is defined such that fairbot will only cheat if on average the previous actions the opponent cheated more
fairbot = function(h) {
  if (length(h)==0) T
  else mean(h,na.rm=T)>=.5
}

# Initially agressive fair bot
ag.fairbot = function(h) {
  if (length(h)==0) F
  else mean(h,na.rm=T)>=.5
}

# Coward dictator. If opponent on average cooperates then cheat, if opponent on average cheats then cooperate.
dictator = function(h) {
  if (length(h)==0) T
  else mean(h,na.rm=T)<=.5
}

# Aggressive dictator
ag.dictator = function(h) {
  if (length(h)==0) F
  else mean(h,na.rm=T)<=.5
}

# Frenemy
frenemy = function(h) {
  if (length(h)==0) T
  else mean(rev(h)[1:2],na.rm=T)==1
  # Frienemy is friendly until you have two interactions in which you cooperate with it then it cheats before turning friendly again if you cooperate twice then it turns ugly again.
}

agent.list = c("always.coop", "always.cheat", "tit4tat" , "ag.tit4tat",
               "fairbot"    , "ag.fairbot"  , "dictator", "ag.dictator",
               "op.tit4tat" , "Random.Tod"  , "frenemy")
             
# First I will draw a random draw of agents.
agents = sort(sample(agent.list, n.agents, replace=T))

# Now let's generate the agent history matrices for each player.
agent.number = matrix(NA, ncol=n.agents, nrow=n.rounds)

# The agent actions matrices are
agent.actions.recieved = matrix(NA, ncol=n.agents, nrow=n.rounds)
agent.actions.given = matrix(NA, ncol=n.agents, nrow=n.rounds)

# Agent point history (the total is the sum of the columns)
agent.points = matrix(NA, ncol=n.agents, nrow=n.rounds)

# Let's start the simulation:
n = 1
for (n in 1:n.rounds) {

  # This while loop ensures that every agent is matched with another agent
  while (sum(is.na(agent.number[n,]))>0) {
  # In each round each agents gets randomly paired up with another agent.
  unpaired = 1:n.agents

    # Loop through all of the agents
    while (length(unpaired)>1) {
      # Select the first element of the unpaired list as a starting place.
      i = unpaired[1]
      # Remove agent i from the unpaired list
      unpaired = unpaired[unpaired!=i]
      # This randomly selects a single agent to match i to.
        # I threw the repeate function in there because when unpaired gets down to
        # only a single number in length then sample starts functioning differently.
        # Thus the rep keeps it from ever being length 1.
      i.match = sample(rep(unpaired,2), 1)
      # Remove i.match from unpaired list
      unpaired = unpaired[unpaired!=i.match]
      # Assign to agent i and i.match their match values
      agent.number[n,i] = i.match
      agent.number[n,i.match] = i
    } # This loop matches agents
  }
  # There should be no NAs in this list if the previous matching loops worked properly.
  sort(agent.number[n,])
 
  for (i in 1:n.agents) {
    # Define each agent's history from previous encounters with the current counterpart
    a.faced = agent.number[n,i]
    a.numbers = agent.number[,i]
    a.actions = agent.actions.recieved[,i]
 
    # Select h as the history of actions with the current rival
    # excluding this round.
    h = a.actions[(a.numbers==a.faced)&(n>1:n.rounds)]
 
 
 
    # Feed history into response function
    response = get(agents[i])(h)
    print(h)
    print(response)

    # Save the response to the actions
    agent.actions.given[n,i] = response
    agent.actions.recieved[n,a.faced] = response
 
    # print(i)
  }
 
}

# Calculate scores by round:
scores.round = returns(agent.actions.given, agent.actions.recieved)

# Final results are:
cbind(agents,apply(scores.round, 2, sum))

# Let's map out the social networks inherent in this exchange:
require(igraph)

# Let's first create out edge list
# Specify a number less than the total number of rounds to map since the total number of rounds is probably a little too large.
map.rounds = 7

x = cbind(1:n.agents, agent.number[1,])
for (n in 2:map.rounds) {
  x = rbind(x, cbind(1:n.agents, agent.number[n,]))
}

# We will use the statnet package to plot the social network
require(statnet)
g = network(x)

plot(g, label=agents, main="A social network graph does not help")

# A better graph might be one that shows how different agents performed over time.
score.cum = apply(scores.round, 2, cumsum)

agent.color = rainbow(length(agent.list))

plot(x=c(1,n.rounds+25),y=c(min(score.cum),max(score.cum)),type="n",
     main="AI performance over time", xlab="round", ylab="score")
for (i in 1:n.agents) lines(1:n.rounds, score.cum[,i], col=agent.color[agents[i]==agent.list], lwd=2)

# Reorder agents by lowest to highest cumscore at the end
agents.ranked = agents[order(score.cum[n.rounds,])]

yrange = max(score.cum)-min(score.cum)
yinc = (yrange/n.agents)
for (i in 1:n.agents) text(n.rounds+17.5, yinc*i+min(score.cum)-2.5, paste(n.agents-i+1,agents.ranked[i]), col=agent.color[agents.ranked[i]==agent.list])

# Each time you run the simulation you get different results depending upon the random draws and the agents in the pool.  Experiment with what happens as the size of the pool increases.  Cooperative strategies loose out to aggressive strategies.

# We can see that up to between round 50-80 all of the bots perform equally since they do not have any history on each other.

# Notice that the opportunistic tit4tat seems to perform the best frequently.  This is kind of unfair since this bot is using meta-game data.  That is, it knows the game is nearly ended so it starts cheating like crazy because there is probably not time for future reprisals.

flattr as an open-source financing solution

Dear readers, many of you are probably wondering what the flattr button is that has appeared on my blog.  I think I have some reasonable arguments for why flattr is a great idea.

First off let’s look at what flattr does.  From superficial description flattr could be easily confused with many traditional forms of non-centralized donation based fundraising operations.  However, flattr rather than asking for money asks for “flattrs” which are used to distribute funds in a monthly payment system.  Thus the flattrer is paying a fixed amount monthly which is redistributed evenly throughout the recipients of the flattrs.

In this case you choose a and n, a being the amount you give monthly and n being the number of projects that you give to.  gamma is the reasonable fixed value of .1 which flattr takes as commission for providing the service.  Thus if you choose to give $10 a month and choose four flattr “things” which is a general term meant to include anything from software developments, music productions, to blogs.  Thus you would be effectively giving $9/4 which is $2.25 to each of your "things", from now on I will call them projects.  

For the individual:

From an individual perspective I think this choice pattern could be optimal for a number of reasons.  For one it reduces transaction costs.  It reduces the time involved in making many small or large donations which might involve inputting credit card information and unpleasant follow ups by organizations.  This transaction cost I believe can be quite burdensome.  Imagine if you would like to sponsor 100 organizations with a total of $200 over a year.  Using traditional methods assuming it takes 10 minutes to sponsor each organization we are talking about 33 hours spent just trying to pay the organizations.  Additionally, organizations might not find it worth the time to handle such small transactions.  Flattr make these kinds of transactions very possible and efficient.

Flattr circumvents that Kantian Categorical Imperative logic which argues, “if you are giving to this person or organization why aren’t you giving to all organizations?”  This question leads you naturally to a scenario where you would be giving too much money to be fiscally sound or the size of the donations are so small as to make transaction costs larger than the payments.  However, using flattr, there is no penalty for increasing the number of projects you are giving to which can be independently modified without changing the amount you are giving.

As an alternative to advertisements:

For media production organizations Flattr provides an alternative to advertisements.  I know some economists argue that advertisements provide signaling information to consumers about the quality of their product and that information has some value to the consumer.  I think there is some validity in this argument however I believe the cost of providing advertisement information to the consumer in almost all scenarios vastly outweighs the marginal benefit of that information to consumers.  Also, that quality is not always signaled but rather other primal impulses are triggered such greed, lust, and gluttony which do not enrich the human experience in general.  In addition, the direct amount of time and energy that advertisements absorb is huge and costly.  For an obvious example, when watching tv, about 1/3 of time spent is watching commercials.  Yet the benefit to the television station for you using up 20 minutes for each hour of your time watching junk advertisements is probably only a few cents.  Yet they do it because it is the only way to extract payment from the user.

Though I have been getting a fair bit of traffic on this blog, I have refused to post any advertisements because of the high cost to users that these advertisements would impose.  Yet, I would like to have some reasonable amount of revenue that I can use to turn around and invest in other projects that I am interested in starting such as a project to develop a free interactive learning environment for students.  I think Flattr provides a reasonable system to raise funds from interested micro-donors.

As a means of covering overhead:

Often times it is very difficult for organizations to raise money to cover overhead costs yet these costs are critical for organizations to function efficiently.  These overhead resources such as the cost of equipment, personnel, advertising, etc. are essential for the daily operations of many organizations yet this is always one of the hardest things for organizations to raise money for.  I very good talk addressing some issues can be found at:

The decentralization of information

The Economist Friedrich Hayek made quite a splash with his seminal paper “The use of knowledge in society”, 1945.  In the paper he argued that a decentralized system for organizing information largely through use of price mechanisms can be far more effective at efficiently allocating resources than a centralized system.  This paper was a direct affront to communist arguments pushing the efficiencies of centralized control and large scale production.  As it turned out, with the fall of the Soviet Union he was shown extremely right.

However, price systems have not shown to be uniformly effective on all scales in all markets.  Some markets such as media production, book publishing, and reporting have resisted decentralization.  Yes, in most industrialized nations there is some separation between these industries and the government, yet these industries all have natural economies of scale which makes them prone to forming natural monopolies.  For instance, you can be a reporter looking into US military misconduct yet if you want to publish anything you need to go through one of the big publishing organizations which might very well reject your article out of self-preservation.  Looking at the destruction rained down on wikileaks, it is clear that censorship is very alive in the world.  Flattry might present a mechanism for small specialized production organizations to develop projects which would be too controversial for established organizations to participate in.

As funding for open-source development

Finally, I think that flattr has the most potential as a source for the funding of open-source project such as R which year after year must compete with highly funded and well-advertised commercial ventures such as Stata, SAS, SPSS, etc.  R for example is one of the most effective open source collaborations ever established.  It has literally tens of thousands of packages that do any number of operations.  It is also extremely powerful in giving the user near complete control of the objects in the environment.

Yet little of that popularity translates into traditional funding resources which could be used to maintain or expand availability of the software (such as maintaining software repositories or promote the software to new users).  This results in far too many users using extremely expensive proprietary software instead of R which I believe is generally also expensive and inferior outcomes.  Fortunately, R has many university level applications such that universities around the world have established R repositories which host R storage.  Even so, if R users gave a dollar or two a month to advertising and development then I think the user base could be rapidly expanded.

Estimating Reliability - Psychometrics

# Each item has a item characteristic curve (ICC) of:
PL3 = function(theta,a, b, c) c+(1-c)*exp(a*(theta-b))/(1+exp(a*(theta-b)))

# Let's imagine we have a test with parameters a,b,c
nitems = 50

test.a = rnorm(nitems)*.2+1
test.b = rnorm(nitems)
test.c = runif(nitems)*.3

# I classical test theory we have the assumption that total variance of the test results is equal to the variance of the true scores (defined as the expected outcome of the test for a given subject) plus measurement error [Var(X) = Var(T(theta)) + Var(E)].

# We "know" that classical test theory is wrong because we know that the error must be a function of ability (theta) since we know that no matter how good (or bad) a student is, the student will never get above the maximum (or below the minimum) score of the test.

# Let's generate the student pool.
nstud = 5000
stud.theta = rnorm(nstud)

# First let's generate our expected outcomes by student.
prob.correct = matrix(NA, ncol=nitems, nrow=nstud)
for (i in 1:nitems) prob.correct[,i]=PL3(stud.theta, test.a[i], test.b[i], test.c[i])

# Generate true scores
true.score = apply(prob.correct, 1, sum)

# Generate a single sample item responses for test 1.
x1.responses = matrix(runif(nitems*nstud), ncol=nitems, nrow=nstud)
# Generate total scores
x1.score = apply(x1.responses, 1, sum)

# The test sampling seems to be working correctly
cor(true.score, x1.score)

# Let's sample a parrellel test
x2.responses = matrix(runif(nitems*nstud), ncol=nitems, nrow=nstud)
# Generate total scores
x2.score = apply(x2.responses, 1, sum)

# Classical reliability is defined as the correlation between two parrellel test scores (I think):
rel.XX = cor(x1.score, x2.score)
  # Looks like our estimated reliability is about .84 or .85

plot(x1.score, x2.score, main="Highly Reliable Test fit along a 45 degree line")

# In theory we can find the reliability by performing the following calculation:
# rho = var(T)/var(x1)
rel.TX = var(true.score)/var(x1.score)
# This seems to give an estiamte that varies between .79 and .88

# Alternatively according to theory rho = 1 - var(E)/var(X)
E = true.score-x1.score
rel.EX =1-var(E)/var(x1.score)
  # Interestingly, this seems to produce a different estimate which is between .82 and .85 typically.

# An alternative reliability estimate is the Kuder-Richardson Formula 20
# K-R20= k/k-1 * (1-sum(p*q)/var(x1))
# With k defined as number of items
# and p probability correct with q = 1-p

# First let's calculate the estimated p by item.
P = apply(x1.responses, 2, mean)
Q = 1-P

# Now for the K-R20
rel.KR20 = nitems/(nitems-1)*(1-sum(P*Q)/var(x1.score))
  # This estimate seems to be between .84 and .86

# In order to calculate the IRT approximations of the classical test theory reliability we first need to calculate information.
# I = p*q when d=1 and p and q are the true probabilities of getting the correct answer.
# Which is equivalent to the prob correct on each item
I = prob.correct*(1-prob.correct)

# Test information can be found as the sum of the information
x1.I = apply(I, 1, sum)

# The unconditional standard error estimate.
SE = 1/x1.I^.5

# Average (unconditional) standard error
uSE = mean(SE)

# An IRT reliability estimate Thissen (2000) could be formed as:
rel.IRT = 1-uSE^2
  # This estimate is around .89.  However, I have not estimated the student abilities which might lead to a different level of reliability if I needed to.

rel.XX
rel.EX
rel.TX
rel.KR20
rel.IRT