Heavy Tails In Python

Below is an exploration of heavy tails using Python, and some of the problems they present for analysis. Heavy tails are distributions with extremely “fat tails”, they have very high likelihood of extreme values relative to a normal bell curve or even a log normal distribution.

Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.


NYC Traffic Accidents Heatmap

A simple chart in one place.

This figure was a nice-looking variant to a paper that was ultimately accepted in EEJ. This figure itself didn’t make it, but it is a really good-looking one. This is a heat map of all accidents in NYC, from July 2012 through March 2019. Black areas are few accidents, and the brighter areas are more accidents. The roads are highlighted in pink to match. It’s fairly obvious that the bulk of accidents are in Manhattan and only a few are in Staten Island.

Mammen, K., Shim, H. S., & Weber, B. S. (2019). Vision Zero: Speed Limit Reduction and Traffic Injury Prevention in New York City. Eastern Economic Journal, 1-19.

Made in QGIS.


A (Quick) Breakdown of SSCAIT Win Rates by AI Race.

Overall, I find that after controlling for individual bot skill level, bots tend to under-perform against Protoss (P) bots and over-perform against Terran (T) bots, and possibly Zerg (Z) bots as well.  Disclaimers, methodology, robustness, conclusion to follow.

This does not mean that the race P is superior or the race T is inferior. It only means bots perform worse than their past performance would suggest against P opponents, and better against T. It is entirely possible that the exploits against P bots are harder to identify, or the bot meta favors strategies that are played by P, (eg, late game death-ball style plays). This prediction is also honed to the average skill level- at tail ends of the skill distribution, these results will vary.


First, let’s dump the data from here into Excel->CSV->R. I had to manually add the races, though.  I’m not really convinced that data types other than CSV aught to even exist, but that’s just me.

Here’s the meat of it:

  1. Get each bots overall win rate for each match. (Here it is an average of the match results, which could be either 100%, 50%, or 0%)
  2. Assume this is win percentage a good proxy for the bot’s overall skill level. This needs a bit of explanation. It is already known the Locutus family is strong, for example, and they are likely to win most match-ups regardless of race (90%+).  We also know they always play Protoss. Let’s try to decompose the win rates into the portion which is approximately due to bot skill level and the part which is due to race match-ups. We need some proxy for bot overall skill level since we are trying to compare the win rates of each race after excluding any effect of skill.
  3. Estimate the following:

MatchWinPct = \beta_p * ProtossOpponent + \beta_t * TerranOpponent + \beta_z * ZergOpponent + \beta_s*SkillProxy + \epsilon

Where the coefficients of interest are going to be: \beta_p, \beta_t,\ and\  \beta_z. Each of these three coefficient represents how different from its average win percentage a given bot is expected to do against an opponent of the given race.  Here are the results:

lm(formula = match_pct ~ skill_level + factor(opp_race) - 1, 
    data = melted_data)

    Min      1Q  Median      3Q     Max 
-90.911 -27.440   1.096  29.510  93.088 

                  Estimate Std. Error t value Pr(>|t|)    
skill_level        0.99908    0.03827  26.107   <2e-16 ***
factor(opp_race)P -5.57628    2.37514  -2.348   0.0190 *  
factor(opp_race)T  4.62700    2.38244   1.942   0.0523 .  
factor(opp_race)Z  1.32735    2.47510   0.536   0.5918    
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 37.52 on 1976 degrees of freedom
Multiple R-squared:  0.6814,	Adjusted R-squared:  0.6808 
F-statistic:  1057 on 4 and 1976 DF,  p-value: < 2.2e-16

First, we see skill level is strongly significant and nearly exactly 1. If a bots win percentage goes up by 1%, we expect it to win each match by precisely 1% more.  This is to be expected.

Next, the coefficient on P is -5.6%. That suggests if a given bot (A) is playing a Protoss opponent, bot A’s overall match score is about 5.6% lower than its typical skill level would suggest, and this is difference is significantly different than 0.

Next, the coefficient on T is +4.6%. That suggests if a given bot (A) is playing a Terran opponent, bot A’s overall match score is about +4.6% higher than its typical skill level would suggest, and this is difference is also significantly different than 0.

The coefficient on Z is positive but small and insignificant. Bots tend to do about average, maybe slightly better than average against Z opponents.  (It may be significantly different than P or T though, I am dropping the constant and comparing coefficients is outside the scope of this quick post.)

A note about R2: R squared only measure predictive power. We do not care about predictive power here in this context, only the estimated effect of the coefficients.  If that cofficient is accurate to the real world, a good deal of accuracy can be sacrificed.

I have also redone this with Logit, which is slightly more accurate on the far tails. (In a journal article on this subject, we’d probably focus more on this result) But here we can simply make some examples. The estimated race effect will be larger near the bottom of the skill distribution and smaller near the top.  Example: Given a bot with the skill level of Purplewave (90.09%), the expected win rate moves about +5% when transitioning from a P opponent to a T opponent.  Given a bot with the skill level of SijiaXu (40.91%), its expected win rate moves about +12% when transitioning from a P to a T opponent.


After controlling for skill level of individual bots (which is the vast and overwhelming factor in game wins), we find some modest suggestion that race is a relevant factor in matchup performance. In general bots under-perform against opponents of the race P and over-perform against the race T, but this factor is small at the higher end of the skill distribution.




#Get that data – copy/paste from online table. Note spaces in names cause some rudeness in R-studio with autocomplete. I did a find and replace in the csv before hand.
SSCAITdata <-“C:/Users/Bryan/Desktop/SSCAITdata.csv”))
# View(SSCAITdata)

#Seperate the Win-losses
melted_data$wins<-as.numeric(sapply(strsplit(as.character(melted_data$value),”\\-“), “[“, 1))
melted_data$losses<-as.numeric(sapply(strsplit(as.character(melted_data$value),”\\-“), “[“, 2))

melted_data$opp_race<-rep(SSCAITdata$Race[present_elements],each = sum(present_elements)) #make sure the races line up
melted_data<-melted_data[!$value),] #drop games with no opponent.

melted_data$match_pct<-melted_data$wins/(melted_data$wins+melted_data$losses) * 100
skill_level<-aggregate(match_pct ~ Bot, data = melted_data, mean)
melted_data<-merge(melted_data, skill_level, all.x=TRUE)

melted_data$matchup_type = paste0(melted_data$own_race,”v”,melted_data$opp_race)
# View(melted_data)

aggregate(outsized_performance ~ matchup_type, data = melted_data, mean) #compare to below:
summary(lm(match_pct~skill_level+factor(own_race):factor(opp_race)-1, data=melted_data))

#with intercept terms like most people expect
summary(lm(match_pct~skill_level+factor(opp_race), data=melted_data))
summary(lm(match_pct~skill_level+I(skill_level^2)+I(skill_level^3)+factor(opp_race), data=melted_data))
summary(glm(I(match_pct/100)~skill_level+factor(opp_race), data=melted_data, family=binomial(link=”logit”)))

#without intercept terms because of clean presentation
summary(lm(match_pct~skill_level+factor(opp_race)-1, data=melted_data))
summary(lm(match_pct~skill_level+I(skill_level^2)+I(skill_level^3)+factor(opp_race)-1, data=melted_data))
summary(glm(I(match_pct/100)~skill_level+factor(opp_race)-1, data=melted_data, family=binomial(link=”logit”)))

melted_data$prediction<-predict(glm(I(match_pct/100)~skill_level+factor(opp_race)-1, data=melted_data, family=binomial(link=”logit”)), melted_data, type = “response”)


Undermind Interview

I talk with a few experts in programming about CUNYBot and the Starcraft-BW programming scene below, as well as a discussion about causal inference in games (forthcoming paper in IEEE-CoG):

Undermind Interview

Undermind Interview #2

Currently CUNYbot has been undergoing a lot of revisions. The code for many sections has been scuttled and redesigned from the ground up. Much “technical debt” has been paid, and the code is approaching clarity. Speed is improved and many already established systems have been incorporated.

Programming, Statistics

Python and Zero-Inflated Models

While this is hardly a tutorial, I’ve been spending a good deal of time working with zero-inflated data for a forthcoming project, and have worked with it extensively in the past.

The point of zero-inflated models is there are ultimately two sources of zeros, zeros can come from the primary model (usually Poisson), or they can be injected into the total model from some other source. For example, casualties from terror attacks may be zero inflated. Consider the idea that either there is no attack on a particular day (0 casualties), or the attack occurs but no one is hurt (also 0 casualties). However, we may want to model the daily casualties as a single model, so we need both parts in one area. (Seminal ZIP paper,  my applied paper on ZIP)

Here’s some code creating ZIP data, and then testing the ZIP model in Python to see if it properly recovers the coefficients I am looking for. Indeed, it does.  Since the ZIP models in python are not designed for the same things economists are typically looking for (ex, no output of standard errors and coefficient estimates in one table), this is a good exploration in to ZIP models, and a way to look around Python, the code du jour.


Load Prerequisites
from matplotlib import pyplot as plt
import numpy as np
import math as math
import pandas as pd
from scipy import stats
import statsmodels.discrete.count_model as reg_models
Generate Mock Data
The key point here in zero inflated (ZI) processes is that there is TWO ways of generating zeros.   The zero can be generated either through the (ZI) or through another process, usually Poisson (P).   Common examples include assembly line failure, the number of crimes in a neighborhood in a given hour.    Critically here was the challenge of indexing Python appropriately.


N = 100000
x_1= np.random.normal(5,1,size=N)
x_2= np.random.normal(5,1,size=N)
poisson_part = np.zeros(N)
zi_part = np.zeros(N)

for i, item in enumerate(x_1):
    poisson_part[i] = np.random.poisson(math.exp(0.2*x_1[i]-0.1*x_2[i])) #needed to initialize the test object. Note the poisson parameter is of the form e^(Bx), ln(lambda) = Bx
for i, item in enumerate(x_1):
    zi_part[i] = np.random.logistic(0.3*x_1[i]-0.2*x_2[i]) > 0 #needed to initialize the test object.

y = zip_model_data = poisson_part * zi_part

Let’s take a look at the data:
          0         1
0  7.212902  5.344758
1  7.128398  4.654915
2  6.841711  4.715646
3  5.082382  5.770086
4  5.858964  3.736356
5  4.173984  3.839683
6  6.157271  5.560635
7  6.375915  7.051355
8  5.943028  4.532561
9  5.847871  5.960908
[0. 1. 4. 3. 0. 2. 0. 2. 3. 0.]
[1. 0. 1. 0. 1. 0. 1. 1. 1. 1.]
[0. 0. 4. 0. 0. 0. 0. 2. 3. 0.]
Generate the model-object itself.
out=reg_models.ZeroInflatedPoisson(y,x,x, inflation='logit')
Fit it:
fit_regularized=out.fit_regularized(maxiter = 100) #essentially forces convergence by penalizing. Biases estimates.
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 1.479267848993134
            Iterations: 10
            Function evaluations: 12
            Gradient evaluations: 10
Let’s take a look at what our (regularized) fit was:
fit_regularized.params ##notice that these are regularized values, not the true values.  The ordinal scale of the variables 
inflate_0   -3.339664
inflate_1   -3.144663
0            0.237893
1           -0.232087
dtype: float64
Let’s try a REAL fit, no regularization, no biasing the estimates.'bfgs', maxiter = 1000) # May need more than the default 35 iterations, very small number!
..../python3.6/site-packages/statsmodels/base/ ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
Optimization terminated successfully.
         Current function value: 1.374978
         Iterations: 13
         Function evaluations: 16
         Gradient evaluations: 16
How did it look?
inflate_0   -0.305607
inflate_1    0.208904
0            0.197614
1           -0.096841
dtype: float64

Nailed it! Regularization had had a noticeable effect on the process, but these are the true coefficients I was attempting to recover.

Note that ZIP in python does not automatically assume the existence of an intercept term. Nor does it automatically calculate the standard errors. I’ll be bootstrapping those soon enough.

Programming, Statistics

AI and Economics

I have just released my first version of an AI to play the classic 1998 Starcraft: Broodwar.  My goal with this project is to combine economic fundamentals and behavioral tendencies with existing machine learning techniques to create a versatile and lifelike opponent.

This is a very active area of research, indeed, a publication was made on this topic yesterday.  Further literature is available in the collections listed here, here, and here.  Major competitors in this field include Facebook, dozens of PhD projects from major universities, as well as some of the best private sector programmers and developers from around the world. I am currently competing in the largest tournament of AI bots, SSCAIT, where you can see my bot play live 24/7 for your viewing pleasure and academic interest.

My AI was also selected for review in the weekly report, and this review gives an excellent demonstration of its macro-economic leanings. In this match, it makes aggressive and early investment into its long-term workforce, and has makes minimal investment in the short-run returns of combat units.   The effect is an an overwhelming crush of units in the mid-to-late game.

In this next match, I play a similar macro style, but aggressively manage complex cloaking and air attacks at the same time.

Regardless, this AI learns from every game, and every opponent is different, so each match is a unique experience.  Expect lots more fun game play from it soon!

Features of the AI:

  • Uses macroeconomic modeling to evaluate the state of its own internal economy.  Maximizes its total output.
  • Learns over time to tune its economy to match each opponent.  For this, I use a genetic algorithm with a long memory and a low mutation rate to slowly move in the right direction as the bot plays dozens and dozens of games.
  • Uses no preexisting map libraries or open source bots, written entirely from scratch in C++.
  • High performance even in full combat situations with many units on both sides. Few lag-outs or crashes (indeed, all on record are my fault for accidentally uploading the bot’s memory files incorrectly).

You can see my bot’s page directly on SSCAIT here.

Other assorted press of my bot are below:


Behavioral Economics and In-class Exercizes

I wanted to spend a moment to report on the incredible success I have been having using in-class exercises. Students have been engaging the content more directly, understanding of the immediate subject matter has improved, and much to my surprise, interest continues onto subjects beyond that immediate exercise.  This continuation of interest for the entire 2-hour class period suggests more than mere amusement is involved, actual interest in the subject has been created. Perhaps the activities have moved their reference point to a more engaged and interested one.

I would like to share some of the excellent exercises that have been done with our student body.  This is not a limiting set of these activities, there are many more available. Many thanks to Virgina Econ Lab for its excellent work in digitizing some of these classical class exercises.

The first exercise, best as an opening exercise, describes the opening idea that individuals are well-described by the rational model, and shows how the rational class naturally arrives at the results that we expect in classical economics.   We have some of the class as suppliers, some of the class as buyers, and random values for the good are assigned. They enter a mock market place and trade goods among themselves.   Output of the program is very good visually (although a little too colorful for my tastes). Learning is very visible.

Possible Output of Supply and Demand Exercise

Total time was maybe 15-20 minutes of execution, 5-10 minutes of setting up and logging into the computer, costs that have reduced after multiple exercises.

I thought it was best to lead the exercise with:

  • Review of PPFs
  • Refresh basic supply and demand model.
  • Review of basic rational model for firms (MC=MR)\ or\ (M\pi =0) .

I thought the gains from the exercise included:

  • More confidence in the idea of rational agents having validity.
    • More interest in rational agents
  • Students put idea of rationality into contrast with behavioral models and procedurally rational models.
  • Strong retention and application of supply and demand, rational model for firms.

We have, since then, conducted several simulations where students interacted directly with rational (and procedurally rational) models. I have programmed these simulations myself and they are available here: Simulating Optimization Problems and Production Possibility Frontiers  and Simulations of Optimization over Time. I would recommend these to any undergraduate body that has easy access to a PC lab, and more detailed interaction with the models is suitable for earlier graduates. I note that the programs do not generate exact correct answers, so one can still ask homework questions on the subject.

Most recently, we have conducted auctions to examine how individuals handle different auction schemes.   Rational results suggest that individuals bid under their valuation in first-price sealed-bid auctions, and bid their valuation in second price sealed-bid auctions. But people do not always perform to those expectations, and it should be illustrated to students.

Some sample output from this program might look like below:

Possible Output of Auction Exercise

We would hope it would clearly show variation between the two auctions, like the plot above.  The second-price auctions lead to significantly higher prices overall, and will often lead students bidding higher than their valuation.

I thought it was best to lead the exercise with:

  • Complete teview of behavioral ideas, particularly:
    • Transactional Utility
    • Risk Aversion
    • Reference Points
    • Mental Accounting

I thought the gains from the exercise included:

  • Better identification among each auction type:
    • Vickry
    • English
    • Dutch
  • Students correctly identified explanations for behavior.
  • Students internalized the optimal bidding value for first-price sealed-bid auctions over uniform distributions, \frac{n-1}{n}

Simulations of Optimization over Time

Behavioral Economics class incorporates a discussion about consumption smoothing, the idea that people prefer gradual changes in consumption over time, and so agents have a smooth relationship between consumption in one period and the next. In fact, it can be shown (with simple log utility functions) that the ratio between consumption today (x_t) and consumption tomorrow (x_{t+1}) is:


Where delta is the “discount rate”, the relative value of today vs tomorrow for the agent.

Below is a simulation I designed to help demonstrate this concept to my class.   It shows an agent struggling to optimize a three-period consumption model.   We always pause to note how the marginal value of consumption smoothly declines every period, and that the discounted marginal utility is nearly the same each period.   (The model simulation residuals are surprisingly large, but nevertheless illustrative.)   The simulation output indicates the consumption in each period with the red line, and a good estimate of the other possible consumption points in black.


#Define Terms
#Checking any options inside the constraint or on the constraint
x1<-runif(100000, min=0, max=y/p1)
x2<-runif(100000, min=0, max=(y-p1*x1)/p2)
#Checking only on the constraint.  Assumes no leftover resources.

#Typical Utility Function

U<-log(x1)+ delta * log(x2) + delta^2 * log(x3)
U2<-delta * log(x2) # undiscounted utility of period 2.
U3<-delta^2 * log(x3) # undiscounted utility of period 3.

plot(x1, U1, ylim=c(0,2.5))
abline(v=x1[which(U==max(U))], col=”red”)
plot(x2, delta*U2, ylim=c(0,2.5))
abline(v=x2[which(U==max(U))], col=”red”)
plot(x3, delta^3*U3, ylim=c(0,2.5))
abline(v=x3[which(U==max(U))], col=”red”)



#Marginal Utility
1/log(x1_star); 1/log(x2_star); 1/log(x3_star);

#Discounted Marginal Utility
1/log(x1_star); delta*1/log(x2_star); delta^2*1/log(x3_star); #Discounted marginal utilities are nearly identical.

Programming, Statistics, Teaching Materials

Simulating Optimization Problems and Production Possibility Frontiers

In teaching Behavioral Economics, optimization problems require some intuition. This intuition can be opaque without calculus literacy.  Below is a simulation to demonstrate that the process for constrained optimization works. It has the added benefit of showing isoquants (by colored stripes in the image below), and the strict boundary condition of the efficiency frontier.

Basic constrained optimization problems are as follows:


y=p_1 x_1 + p_2 x_2

I have made code in R to simulate the results for a two-part optimization process. The example uses U=sqrt(x_1) + sqrt(x_2) as the functional form.


x1<-runif(25000, min=0, max=y/p1)

#Checking any options inside the constraint or on the constraint
x2<-runif(25000,min=0, max=(y-p1*x1)/p2)

out<-mesh(x1, x2, U)
points3D(x1, x2, U, xlab=”x1″, ylab=”x2″, zlab=”Utility”, phi=-90, theta=90)

plot(x1, U)
abline(v=x1[which(U==max(U))], col=”red”)

And it outputs the following plots (with minor variation). Note that the colored bands represent utility curves, isoquants. The end of the colored points represents the efficiency frontier.

This slideshow requires JavaScript.


The actual solution is found by:


subject to:

y=p_1 x_1+p_2 x_2

The Lagrangian is then:

L=sqrt(x_1)+sqrt(x_2) + \lambda (y - p_1 x_1 - p_2 x_2)

Leading to the first order conditions (derivatives of L):

L_1 : 0.5 x_1^{-0.5} - \lambda p_1=0

L_2: 0.5 x_2^{-0.5} - \lambda p_2=0

L_{\lambda} : y- p_1 x_1 -p_2 x_2 =0

Using these 3 conditions, we can find the equations:

\frac{0.5 x_1^{-0.5} }{0.5 x_2^{-0.5}} = \frac{p_1}{p_2}
y- p_1 x_1 -p_2 x_2 =0

Where, if y=10, p_1=1, p_2=2 then we can solve for the theoretical solutions: x_1=6.66, x_2=1.66

These indeed match very closely with the real solutions.

Programming, Statistics

The Optimum Pokemon Portfolio and Principal Component Decomposition (PCD) using R

I have very recently completed the Stanford Lagunita online course on Statistical Learning, and Tibrishani & Hastie have taught me a great deal about Principal Components.  No learning is complete without exercises, however, so I have found a wonderful data set that seems popular, the attacks and weaknesses of Pokemon.  (I am, admittedly, not a pokemon player, so I have had to ask others to help me understand some of the intricacies of the game.)

Principal Component Decomposition:

First and foremost, principal component decomposition finds the direction that maximizes variation in the data.  At the same time, this can be said to be the eigenvalue of the data, the direction which best describes the direction of the data.
For example, if there is a spill of dirt on a white tile floor, the direction of the spill (eigenvalue) would always be the direction the dirt is most widely spread (principal component).

After looking at the beautiful charts used in the link above, I realized this would be very interesting to do a PCD on. What Pokemon are most similar and which are most different in terms of strengths and weaknesses? To find out we will break it into its principal components, and find out in which directions the data is spread out.

Pokemon can vary along 18 dimensions of strengths and weaknesses, since there are 18 types of Pokemon. This means there can be up to 18 principal components. We are not sure which principal components are useful without investigation. We show below how much variation is explained by each type of Pokemon. There doesn’t appear to be any clear point where there the principal components drop off in their usefulness, perhaps the first 3 or the first 5 seem to capture the most variation.  The amount of variation captured by each principal component is outlined below.
Let us now look at the principal components of the Pokemon attack/weakness chart directly.  We can visualize them in a biplot, where the arrows show the general attacking direction of the pokemon and the black labels show the defending labels.  The distance from the center of biplot shows the deviation of that pokemon type from the central eigenvalue/principal component.  Labels that are close together are more similar than those further apart.


So for example, Ghost attacks (arrows) are closely aligned with Ghost defence (black label) and Dark defence (black label).  In general, the Pokemon that are most different in defence is Fighting and Ghost, and still again distinct from Flying and Ground defence.  This suggests that if you wanted a Pokemon portfolio that would be very resilient to attack, you would want Fighting/Ghost types.  If you want a variety of attacks, you might want to look into Ghost/Normal types or Grass/Electric.

Keep in mind together these only explain about 35.5% of the variation of Pokemon types, there are other dimensions in which Pokemon vary.  I expected fire and water to be more clearly different (and they are very distinct, they go opposite directions for a long distance from the center!), but they are less distinct than ghost/normal.

The Optimum Pokemon Portfolio:

This lead me to wonder what type of pokemon portfolio would be best against the world, something outside the scope of the Statistical Learning course but well within my reach as an economist.  Since I don’t know what the pokemon-world looks like, I assumed the pokemon that show up are of a randomly and evenly selected type. (This is a relatively strong assumption, it is likely the pokemon encounters are not evenly distributed among the types).  The question is then, what type of pokemon should we collect to be the best against a random encounter, assuming we simply reach into our bag and grab the first pokemon we see to fight with?

First, I converted the matrix of strengths and weaknesses above into one that describes the spread of the strength-weakness gap, that is to say, if Water attacks Fire at 200% effectiveness, and defends at 50% effectiveness, a fight between the Water and Fire is +150% more effective than a regular pokemon attack (say Normal to Normal or Ice to Ice). Any bonuses a pokemon may have against its own type was discarded, because it would be pointless.  The chart for this, much like the wonderful link that got me the data in the first place, is here, where red is bad and blue is good: Rplot12

Then I added the strength-weakness gap together for each type of pokemon, which assumes that the pokemon are facing an a opponent of a random type.  According to this then, the most effective type of pokemon are on average:

Type                              Effectiveness
Steel                               0.22222222
Fire                                0.11111111
Ground                              0.11111111
Fairy                               0.11111111
Water                               0.08333333
Ghost                               0.08333333
Flying                              0.05555556
Electric                            0.00000000
Fighting                            0.00000000
Poison                             -0.02777778
Rock                               -0.02777778
Dark                               -0.02777778
Ice                                -0.08333333
Dragon                             -0.08333333
Normal                             -0.11111111
Psychic                            -0.11111111
Bug                                -0.11111111
Grass                              -0.19444444

That is to say, Steel pokemon, against a random opponent, will on average be 22% more effective.  (This is the mean, not the median.) And against a random opponent a Grass pokemon will be expected to be 19% less effective than a Fighting pokemon, shockingly low. Amusingly, Normal pokemon are worse than normal (0) against the average pokemon.

This does not mean you ONLY want Steel pokemon because you could come up with an opponent that is strong against Steel. Nor do you want to entirely avoid Grass pokemon, since they are very strong against many things that Steel is weak against. Merely that if you’re willing to roll the dice, a Steel pokemon will probably be your best bet.  Trainers do not want to take strong risks, trainers are risk averse.  You want to maximize your poke-payoff while minimizing how frequently you face negatively stacked fights. The equation for this is:

Maximize: \ \mu * vars - \delta * t(vars) * cov * vars + \lambda*(1- t(ones) * vars) \ wrt. \ vars

Where \mu is your vector of payoffs in the table above, \delta is your risk aversion, cov is the covariance matrix of the differenced pokemon data set, and vars is your portfolio selection which must add up to one hundred percent.

How risk averse are you?  You could be very risk averse and want to never come across a bad pokemon to fight, or you could love rolling the dice and only want one type of pokemon. So I have plotted the optimal portfolio for many levels of risk-tolerance.  It is a little cluttered, so I have labelled them directly as well as in the legend.


The visualization is indeed a little messy, but as you become more risk averse, you add more Electric, Normal, Fire, and Ice pokemon (and more!) to help reduce the chance of a bad engagement.  In order to do this, one reduces the weight we put on Steel, Ground, and Fairy pokemon, but doesn’t eliminate them entirely.  Almost nothing adds Dragon, Ghost, Rock. or Bug pokemon, they are nearly completely dominated by other combinations of pokemon types.

I’ve plotted two interesting portfolios along the spectrum of risk aversion below. They include one with nearly no risk aversion (0.001), and one with high risk aversion (10).

This slideshow requires JavaScript.

Of course, most importantly of all, regardless of your Pokemon and your interest in being “the very best”, you should still pick the coolest Pokemon and play for fun.

Code is included below:

#Data from:
#write.csv(chart, file="/home/bsweber/Documents/poke_chart.csv")
# library(devtools) 
# install_github("vqv/ggbiplot", force=TRUE)

differences <- (poke_chart-1) - (t(poke_chart)-1)
core <- poke_chart

poke_pcd<-prcomp(core, center=TRUE, scale=TRUE)
plot(poke_pcd, type="l", main="Pokemon PCD")

poke_palette<-c("#A8A878", "#EE8130", "#6390F0", "#F7D02C", "#7AC74C", "#96D9D6", "#C22E28", "#A33EA1", "#E2BF65", "#A98FF3", "#F95587", "#A6B91A", "#B6A136", "#735797", "#6F35FC", "#705746", "#B7B7CE", "#D685AD")

ggbiplot(poke_pcd, labels= rownames(core), ellipse = TRUE, circle = TRUE, obs.scale = 1, var.scale = 1) +
  scale_color_discrete(name = '') +
  theme(legend.direction = 'horizontal', legend.position = 'top')
  #Score plot is for rows, attack data. loading lot is for columns, defense data.  So bug and fairy have similar attacks (shown by rays), similar defences (shown by points). Ghost and normal have almost identical defences, but different attacks.
ggbiplot(poke_pcd, labels= colnames(core), ellipse = TRUE, circle = TRUE, obs.scale = 1, var.scale = 1, choice=c(2,3)) +
  scale_color_discrete(name = '') +
  theme(legend.direction = 'horizontal', legend.position = 'top')  #Score plot is for rows, attack data. loading lot is for columns, defense data.
ggbiplot(poke_pcd, labels= colnames(core), ellipse = TRUE, circle = TRUE, obs.scale = 1, var.scale = 1, choice=c(5,6)) +
  scale_color_discrete(name = '') +
  theme(legend.direction = 'horizontal', legend.position = 'top')  #Score plot is for rows, attack data. loading lot is for columns, defense data.
ggbiplot(poke_pcd, labels= colnames(core), ellipse = TRUE, circle = TRUE, obs.scale = 1, var.scale = 1, choice=c(7,8)) +
  scale_color_discrete(name = '') +
  theme(legend.direction = 'horizontal', legend.position = 'top')  #Score plot is for rows, attack data. loading lot is for columns, defense data.

cov_core<- t(differences-mean(differences)) %*% (differences-mean(differences)) #Make the Cov. Matrix of differences. 
cov_core[order(diag(cov_core), decreasing=TRUE),order(diag(cov_core), decreasing=TRUE)]
vars<-as.matrix(rep(1/18, times=18))
mu<-t(as.matrix(apply(differences/18, 1, sum))) #Average rate of return over 18 pokemon types.

data.frame(mu[,order(t(mu), decreasing=TRUE)]) #Table of Pokemon Types

delta<- 1  #risk aversion parameter

  out<- matrix(0, nrow=0, ncol=18)
  for(j in 1:1000){
  Dmat <- cov_core * 2 * delta
  dvec <- mu
  Amat <- cbind(1, diag(18))
  bvec <- c(1, rep(0, 18) )
  qp <- solve.QP(Dmat, dvec, Amat, bvec, meq=1)
  out<-rbind(out, round(pos_answers, digits=3))
  df <- data.frame(x=1:nrow(out))
    df.melted <- melt(out)
      colnames(df.melted)<-c("Risk_Aversion", "Pokemon_Type", "Amount_Used")
 qplot(Risk_Aversion, Amount_Used, data=df.melted, color=Pokemon_Type, geom="path", main="Pokemon % By Risk Aversion") + 
   # ylim(0, 0.175) +
   scale_color_manual(values = poke_palette) +
#   geom_smooth(se=FALSE) +
   geom_text_repel(data=df.melted[df.melted$Risk_Aversion==8.5,], aes(label=Pokemon_Type, size=9, fontface = 'bold'), nudge_y = 0.005, show.legend = FALSE) 

  # Another plot that is less appealing
  # matplot(out, type = "l", lty = 1, lwd = 2, col=poke_palatte)
  # legend( 'center' , legend = colnames(core), cex=0.8,  pch=19, col=poke_palatte)
  pie(head(out, 1), labels= colnames(out), col=poke_palette)
  pie(tail(out, 1), labels= colnames(out), col=poke_palette)
  df_1<-data.frame(matrix(out[1,], ncol=1))
  ggplot(data=df_1, aes(x=Pokemon_Type, y=Percentage, fill=Pokemon_Type))+
    geom_bar(stat="identity", position=position_dodge()) +
    scale_fill_manual(values = poke_palette)+
    ggtitle("Pokemon Portfolio With Almost No Risk Aversion")
  ggplot(data=df_2, aes(x=Pokemon_Type, y=Percentage, fill=Pokemon_Type))+
    geom_bar(stat="identity", position=position_dodge()) +
    scale_fill_manual(values = poke_palette) +
    ggtitle("Pokemon Portfolio With Very Strong Risk Aversion")
  cov_core[order(diag(cov_core), decreasing=TRUE),order(diag(cov_core), decreasing=TRUE)]
  melt_diff$value<- factor(melt_diff$value)
  simplepalette<-colorRampPalette(c("red", "grey", "darkgreen"))
  ggplot(data = melt_diff, aes(x=Var1, y=Var2, fill=value) ) + 
    scale_fill_manual(values=simplepalette(9), breaks=levels(melt_diff$value)[seq(1, N, by=1)], name="Net Advantage" )+
    ggtitle("Net Pokemon Combat Advantage")+
    xlab("Opponent") +
    ylab("Pokemon of Choice")