Programming, Statistics, Teaching Materials

Simulating Optimization Problems and Production Possibility Frontiers


Basic constrained optimization problems are as follows:


y=p_1 x_1 + p_2 x_2




#Checking any options inside the constraint or on the constraint

out<-mesh(x1, x2, U)

plot(x1, U)


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:


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




Using these 3 conditions, we can find the equations:

y- p_1 x_1 -p_2 x_2 =0


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.
  scale_color_discrete(name = '') +
  scale_color_discrete(name = '') +

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))+
    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))+
    scale_fill_manual(values = poke_palette) +
    ggtitle("Pokemon Portfolio With Very Strong Risk Aversion")
  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")

Multiple Linear Regression in R

This slideshow requires JavaScript.

In the previous exercise: Why do we need N-2?, I show a simple 1 dimensional regression by hand, which is followed by an examination of sample standard errors.  Below I make more extensive use of R (and an additional package) to plot what linear regression looks like in multiple dimensions. This generates the images above, (along with several others).  This illustrates that linear regression remains flat even in N dimensions, the surface of the regression is linear in coefficients.

As a class exercise, I ask that you consider different pairs dependent variables that are  functions of one another. What happens if the function is linear? What happens if the function is nonlinear, for example, cos(x_1)=x_2? Examine what happens to the surface of your regression as compared to the shape of the relationship you are investigating.  Is there a way you can contort the regression estimate into a curved surface to better match?  Why or why not?

install.packages(“plot3D”) # we need 3d plotting
library(“plot3D”, lib.loc=”~/R/win-library/3.1″) #Load it into R’s current library, may vary by computer.

set.seed(2343) #ensures replicatation. Sets seed of random number generators.
n<-25 #number of samples
x_1<-rnorm(n) #Our x’s come from a random sampling of X’s.
b_2<-(-3) # But student transit programs are associated with a decline in crime.
y<-b_0+b_1*x_1+b_2*x_2+u #This is defining our true Y. The true relationship is linear.

#look at data in each dimension
#look at data overall
points3D(x_1,x_2,y,xlab=”x_1″,ylab=”x_2″,zlab=”y”,phi=5) #look at data. phi/theta is tilt.

fit<-lm(y~x_1+x_2)  #fit it with a linear model, regressing y on x_1, x_2

#Make a surface
x_1.pred <- seq(min(x_1), max(x_1), length.out = n)
x_2.pred <- seq(min(x_2), max(x_2), length.out = n)
xy <- expand.grid(x_1=x_1.pred, x_2=x_2.pred)
y.pred <- matrix (nrow = n, ncol = n, data = predict(fit, newdata = data.frame(xy), interval = “prediction”))

summary(fit) #view output of variables.

fitpoints<-predict(fit)  #get predicted points, needed to make a surface.

scatter3D(x_1,x_2,y,xlab=”x_1″,ylab=”x_2″,zlab=”y”,phi=5 , surf=list(x = x_1.pred, y = x_2.pred, z = y.pred, facets = NA, fit = fitpoints)) #look at data. phi/theta is tilt.
scatter3D(x_1,x_2,y,xlab=”x_1″,ylab=”x_2″,zlab=”y”,phi=45, surf=list(x = x_1.pred, y = x_2.pred, z = y.pred, facets = NA, fit = fitpoints)) #From straight on it is a flat plane, residuals are highlighted
scatter3D(x_1,x_2,y,xlab=”x_1″,ylab=”x_2″,zlab=”y”,phi=30, surf=list(x = x_1.pred, y = x_2.pred, z = y.pred, facets = NA, fit = fitpoints)) #From other angles it is clear it is somewhat straight.
scatter3D(x_1,x_2,y,xlab=”x_1″,ylab=”x_2″,zlab=”y”,phi=60, surf=list(x = x_1.pred, y = x_2.pred, z = y.pred, facets = NA, fit = fitpoints)) #look at data. phi/theta is tilt.


Programming, Statistics, Teaching Materials

Why do we need n-2? An example in R



####Coding in R


n<-25 #number of samples





plot(x,hat_u) # Let’s see our residuals
hist(hat_u) # Let’s see our histogram


(sigma_sq-sigma_sq_naive) > (sigma_sq-sigma_sq_hat)








Below is the final forecast, announced 5/16:
















5_16_acf 5_16_pacf


Commentary on Other predictions Types and Styles:




Historical Forecasts




Revisiting the Forecast (5/3) -May 16th



Initial Forecast (4/24) 
 – May 8th.


#Load Dependencies



#Source of data:

#Get opening and closing time to sequence data.


#let us collapse this
#Need objects as xts:
t15 min.

# we will add a lowest observed badge marker.
#let’s get these factors as dummy variables.

#Seasons matter. I prefer Fourier Series:
fourier {
n X for(i in 1:terms)
colnames(X) return(X)



#forecast value of button:
#size of forecast
w n<-7*24*60/t



Programming, Statistics







detrended<-data-xhat #also, norm.
for(n in 2:N){



Programming, Statistics

Music and Math



