Tuesday: I only crashed a computer twice today.

Not like I was running stuff that I really wanted finished or anything.

I didn't think I'd have a picture for today, so I took this when I got home.  I totally wasn't trying to avoid the awkward thing of opening my door while my neighbor opens his door.

But, then I saw that this administration continues to be worse than useless, and is planning to investigate to see if white people are being discriminated against in college applications.  Which is dumb, so I did the standard thing of "asking the government about the data, since I'm pretty sure the government has the data" and:

Thanks, competent part of the government!

Lines show fall college enrollment by broad racial category, taken from this table.  Dots show total US demographics, taken from the last table of this wikipedia page, because I'd already parsed and formatted one government table, and this one already had the same broad categories.  The main caveat is that looking at the subsequent section, you can see that younger age categories are less white.  However, I was not able to find a full year/age/race census split, so it's hard to plot a "college age" population.  Still, broadly speaking, white people have historically been over-represented in fall college enrollment, as has the Asian/Pacific Islander category, although that over-representation drops when you note the recently available Two-or-more category.  Conclusion: this administration can't do an hour's worth of work to discover their thesis is likely incorrect.  Alternate conclusion:

I am disappointed that this image makes a return.

R bit (removing all the mistakes, and me learning about the wonder of reshape2.melt, which takes "wide" tables where each column is a category sampled at each index/year, and turns it into a "long" table with each value indexed by the index and category):

p306 = data.frame(read.csv("./tab306_T.P",header=TRUE,sep='\t'))

library(reshape2)

demo = data.frame(read.csv("./demo.dat",header=TRUE,sep='\t'))
demomelt = melt(demo,id.vars="year")

p306X = p306[, c("year", "White", "Black", "Hispanic", "AsianPacificIslander", "AmericanIndianAlaskaNative", "TwoOrMore")]
p306meltX = melt(p306X,id.vars="year")

ggplot()  + geom_line(data=p306meltX,aes(x=year,y=value,color=variable)) + geom_point(data=demomelt,aes(x=year,y=value,color=variable))
ggsave("college_demographics.png")

• Rabbit Rowlet.
• Dogs.
• Weather clustering.  I vaguely remember noting the "large gradient"/"low high temperature" correlation back when weatherspark was a thing.

And tumblr broke xkit again, so my usual evening flow of

1. Read RSS.
2. Save interesting items.
3. Read through saved items, unsaving them as I open in new tab.
4. Like reblog tumblr items as I close them.

doesn't work.

Monday: Ok, but I needed a figure.

I spent a large chunk of the day proofreading a paper (that is much less complete than I thought it would be), and figures 5-12 were all various correlation plots.  Are "correlogram" plots not common (in quotes because google doesn't converge on if that's the right name)?

The NSF data again, which is a terrible data set to use for this, but I needed something.

The idea being that you have N variables, and want to see how they map against each other, so you put them into a big grid so you can see how all the things match against all the other things.  I wrote a quick version years ago, and went with the R version, just because that was reasonably easy to do:

install.packages("corrgram")
library(corrgram)

corrgram(nsf,upper.panel=panel.pts,lower.panel=panel.ellipse,order=TRUE)

png("why.png")

I'm sure there's a snek version as well.  That reminds me that I need to add a "snek bits" tag.

• My initial guess would have been the kid had been a jerk to the squirrel.  I also would have said that the guy might have cut the video to just show the bite, and maybe he did something before that.  However, putting them both together leads me to the same conclusion as the health department; this isn't a normal squirrel, and it's probably sick.  The only time I've ever been bitten by a squirrel I gained the proportional strength and abilities of a squirrel it was a tiny nibble by a squirrel who was really excited about food, and missed when it tried to get some.
• Ok, if this had been any other site telling me to "not order a fancy burger," I'd ignore it and put it in the "don't tell me what to do" category.  The problem is, I watch their youtube channel, and just last month, they had a video that was "this is a super fancy burger, and you should totally order it if you can."
• Also, I'm pretty sure this qualifies as the end state of "don't tell me what to do" in relation to food.
• I like the raincoat one the best.
• I couldn't tell what was wrong when I watched the bootleg copy of the trailer, but it's obvious when someone says it.  Thanos looks weird without his hat.
• Nova and Quasar can both fly, and they're spending their time making horses angry?
• On the other hand, Patsy Walker, Hellcat, playing frisbee with the Hulk is exactly what I would expect.
• This has been all over my Twitter timeline.  I don't understand how hard this is for people.  "Well, we'll have to see what the writers come up with!"  Boom.  You've diffused the ship talk you don't want to deal with, and you don't make something that makes the shippers feel bad.
• In case it isn't clear, I'm talking about the fandom relationship shipping, not the logistics of transporting objects shipping.

Wednesday: No, it's fun to discover a massive problem with a database just when you want to go home.

But I think I was able to fix it with a redo, and changing all the "INSERT"s to "INSERT IGNORE"s.

And I wanted to make a proper version of the garbage scribble at the top of this post.

I really like the way R plots things.

 It's reasonably simple:

b = data.frame(N = seq(0,177))
b$P = choose(177,b$N) * 0.24**b$N * (1 - 0.24)**(177 - b$N)

L = rbind(c(0,0), subset(b, b$N <= 36),  c(36,0))
U = rbind(c(44,0), subset(b, b$N > 43), c(177,0))

ggplot(b,aes(x = N, y = P, color="P_expect = 0.24")) + \
  geom_polygon(data=L,aes(x = N, y = P, fill="Pobs")) + \
  geom_polygon(data=U,aes(x = N, y = P, fill="Pover")) + \
  geom_line() + geom_point()  + \
  scale_colour_manual(values=c("black"),name="",
                      guide=guide_legend(override.aes=aes(fill=NA))) + \
  scale_fill_manual(values=c("red","blue"),name="")

ggsave("/tmp/2017b.png")

The only major issues are that doing the shading requires constructing a polygon object, and that needs to have endpoints set correctly (or it shades the wrong way).  Getting all the colors and legend set was also not super obvious, and that "override.aes" thing is just nonsense.

"But it's not snek."  So then it was a challenge to see how to do this in snek, too.

Snek is simpler, but their documentation is far worse.  Examples should start slow, not alphabetically with "animation."  Why would you do that?

#!/usr/bin/env python3                                                                                                       
import matplotlib.pyplot as plt
import scipy.special
import numpy as np

P_expect = 0.24
N_talks  = 177
N_obs    = 36
N_expect = 44

N = np.arange(0,N_talks,1)
P = scipy.special.binom(N_talks,N) * P_expect**N * (1 - P_expect)**(N_talks - N)

# Lower portion
Nl = np.arange(0,N_obs,1)
Pl = scipy.special.binom(N_talks,Nl) * P_expect**Nl * (1 - P_expect)**(N_talks - Nl)
Zl = Nl * 0.0

# Upper portion
Nu = np.arange(N_expect,N_talks,1)
Pu = scipy.special.binom(N_talks,Nu) * P_expect**Nu * (1 - P_expect)**(N_talks - Nu)
Zu = Nu * 0.0

plt.grid()
label_text = "P_expect = %.2f" % P_expect
plt.fill_between(Nl,y1=Zl,y2=Pl,color="red",label="P_obs")
plt.fill_between(Nu,y1=Zu,y2=Pu,color="blue",label="P_over")
plt.plot(N,P,color="black",label=label_text)
plt.scatter(N,P,s=5,color="black")
plt.legend()
plt.savefig("mpl.png")

Having an explicit fill_between() function saved a lot of time.  The legend() function was also helpful for making it just work.

I'm still behind on my RSS stuff.

• "Oh no!  This new Spider-Man movie confuses the timeline!"  I complain about stupid stuff, but this is too far.  Comic book timelines have been insane forever.  I mean, look at Squirrel Girl.  Doreen was 14 when she defeated Doom with Tony, then she did stuff for a few years with the GLA, then she was kind in the Avengers, then she babysat Luke Cage/Jessica Jones' daughter, and now she's in college.  How old is she?  Why isn't she 39 if the comics follow regular time?  Doesn't matter, she's doing college now, enjoy your wonderful stories.  Comic book time is meaningless.
• Also they cast Doreen, and this seems like it'll work.
• Wonder Woman.
• Whoops.
• Best Spider-Man.
• Buffalo.

Thursday: No seriously, this week is all messed up for me now.

I spent part of the day thinking it was Monday, because we had a meeting, and those are usually on Monday.  But it's not, so tomorrow I have to try and sort out things that I want to finish this week.

Also, I learned that the answer to "I have all this crap in a data frame, and I want answers for each subset, but I don't want to have to do any work" is ddply.  Originally I did this with work data, but using the NSF PhD statistics works as well.

library(plyr)

nsf = data.frame(read.table("./matched.dat",sep='\t',header=TRUE))
nsf$R = nsf$female / (nsf$male + nsf$female)

q = ddply(nsf, .(broad_field), function(x) c(m = mean(x$R),s = sd(x$R), quantile(x$R, c(0.0, 0.5, 0.75, 0.90, 0.98, 1.0))))

                           broad_field         m           s

1                            Education 0.6780287 0.012619293

2                          Engineering 0.2169928 0.014537864

3                  Humanities and arts 0.5035801 0.007278568

4                        Life sciences 0.5384433 0.017871779

5    Mathematics and computer sciences 0.2479207 0.009173072

6                               Otherb 0.5059506 0.011422802

7 Physical sciences and earth sciences 0.3117802 0.018656152

8       Psychology and social sciences 0.5832883 0.011162981

        50%       75%       90%       98%      100%

1 0.6816616 0.6865097 0.6927421 0.6930276 0.6930990

2 0.2220837 0.2279097 0.2304950 0.2320948 0.2324947

3 0.5058046 0.5084105 0.5096879 0.5120727 0.5126689

4 0.5453048 0.5518069 0.5544373 0.5564318 0.5569305

5 0.2466649 0.2529897 0.2614564 0.2629909 0.2633745

6 0.5105890 0.5143730 0.5165125 0.5207457 0.5218040

7 0.3140528 0.3249957 0.3339903 0.3353106 0.3356407

8 0.5849175 0.5890918 0.5946965 0.5973446 0.5980066

Ok, so the quantile data isn't super useful with the NSF data, but still.  Means and sigmas for each factor, and for the work data, pulling the quantile stuff was useful.  Then, hubris took hold, and I wondered if I use this to do linear fits to each factor as well.  The answer is no, not with ddply, because that reads and writes a data frame, and lm outputs a model object.  So you have to use dlply to save those models in a list, and then make a data frame from the stuff you care about in the list:

f = dlply(nsf, .(broad_field), lm, formula = R ~ year)

coeffs = ldply(f,coef)

coeffs$x2010 = (coeffs[,2] + 2010 * coeffs[,3]) * 100

coeffs$m = coeffs[,3] * 100

                           broad_field (Intercept)          year    x2010           m

1                            Education   -5.318890  0.0029835418 67.80287  0.29835418 

2                          Engineering   -7.627153  0.0039025598 21.69928  0.39025598 

3                  Humanities and arts   -1.145504  0.0008204401 50.35801  0.08204401 

4                        Life sciences   -9.438112  0.0049634604 53.84433  0.49634604 

5    Mathematics and computer sciences    1.055520 -0.0004017907 24.79207 -0.04017907 

6                               Otherb   -3.713211  0.0020990852 50.59506  0.20990852 

7 Physical sciences and earth sciences  -10.098766  0.0051793762 31.17802  0.51793762 

8       Psychology and social sciences   -4.194659  0.0023770883 58.32883  0.23770883 

And yes, x2010 is the same as the mean above, and m (percent change per year) is approximately 1/5 of the sigma above.  Taking a quick average of just the math, engineering, and physical sciences data for 2017 gives something like 28%.  The lack of improvement in math is a problem, and engineering starts at a deficit.  The depressing thing is that even with the highest improvement rate for the physical sciences, it's still something like 40 years until parity.

Also a plot:

library(ggplot2)

ggplot(nsf,aes(x=year,y=R,color=broad_field)) + geom_point() + geom_line() + geom_smooth(method=lm)

ggsave("/tmp/with_fits.png")

Yes, I'm sure there's a line thickness parameter I could have changed.  Whatever.

• More ggplot information.
• Hrm.
• That sandwich thing popped up again.
• Sits.
• I like this, because even though I have no clue what fandom this is from, I can pick up the entire story just from the pictures and the caption.
• Just some Wonder Woman.

Tuesday: My plan for today was to eat hot dogs, take a walk, eat more hot dogs, and watch some fireworks.

I did all those things, and did more editing for SciPy stuff.  I do not like how Snek handles Unicode, because that handling appears to consistently be "be as difficult about everything as possible."

I also drew a very bad picture of how different parts of a probability distribution are included in a calculation.

 And I battled more with my Magikarp, so I have more data to include in the plot from yesterday:

This is the "most useful" version.

This splits two data points into their true input JP and the 5% and 25% scaled values.  I had boosts on those two runs, so the scaled values are probably more correct.

I also need a better way to put code into blogs.  I've decided I'm going to copy my R code when I make things like this, so I can find them later.  It really is super easy:

j = data.frame(read.table("./jp",sep='\t',header=TRUE))
ggplot(j,aes(x=jp,y=height)) + geom_point() + stat_smooth(method='lm',formula=y ~ sqrt(x), col='red')
ggsave("/tmp/jp2.png")

And I took a walk:

The moon was out.

 I saw four geckos and one dog that threw up.  I was happier for one of those.

This didn't quite capture all the pinks and purples.

And then I came home and had my code insulted.  If the input data wasn't so bad, the code would do a better job on it.  More webpages should be designed with a structured data concept.

Then there were fireworks:

Boom!

Ka-ka-boom!

Ba-fszzzzzzzz.

Boom boom boom siziziziziz.

Boom ka boom!

Boo boom!

Boom boom boom boom!

Boom ka boom boom!

Boom boom boom boom boom boom boom boom!

Fazzzzzzzzzzzzzzzzzzz.

• Bee facts.
• Squirrel.
• Squirtle.
• I'm sad that it's ending, but he has to do what he thinks is best.
• I don't have anything to add.

Monday: I should have just taken the four day weekend. I would have been about as productive.

However, last night, I was able to plot up and fit how the height a  Magikarp jumps scales with JP:

JP/1000, obviously.  But my "this feels like it's a square root" idea seems to be supported.  This is just one pass through the league I'm currently battling.

And today I finally decided that I wasn't going to find an easy way to reshape this table of PhD recipients, so I just wrote a perl script to do it.  Once it was in nice "field/year/numbers" format, I was able to get R to make a pretty plot:

So I guess the 24% value probably isn't too far off.

f = data.frame(read.table("./matched.dat",sep='\t',header=TRUE))

ggplot(data=f,aes(x=year,y=(female / (female + male)), color=field)) + geom_line(aes(group=field)) + geom_point()

• A good story.
• "Millions of identical human men".
• Cocoa.
• I am glad to see that Squirrel Girl and Koi Boy being friends is backed up by nature.

Tuesday: Damn it.

I didn't get my talk all finished today, I was going to pack tonight and didn't do that, and I have to get up stupid early tomorrow for a doctor appointment.  What did I do instead?

I did sentiment analysis of my archive of tweets, and then learned about the aggregate function to break it up into monthly averages.

Wait, not only are my tweets generally positive, it's been getting increasingly positive?  What?  That doesn't seem right.

I liked how the moon looked this evening.

• It's all the fault of this blog post.  And for R for being reasonably easy to impulse install packages.
• I've had this happen before.
• Secret birb.
• What?
• That's good.

Sunday: The last day off for two weeks.

So, of course, I spent a chunk of the evening queuing data to download at work, and launching some processing jobs.  This doesn't work well with my plan to "don't do work stuff as much as possible."

I discovered today that because I've been walking, I've apparently lost weight.  That means I now have pants that do not fit because they're too big.  Whoops.  I can't wait until we get to the world completely run by millennials, because at that point, everyone will stop being dumb and we'll all just wear whatever's comfortable, because who wants to wear shit that isn't soft and stretchy?

Today's actual main plan was to go to Target.  Sunday is always a bad day to go to Target, because everyone is there.

Then sushi, because I honestly couldn't decide what I wanted, and unagi kind of sounded good.

It was super good.  Today also marks the end of my reign of terror at Dole plantation:

What?

Someone finally killed Cloyster.

Either that or there's some check in place to make sure you don't hold a gym for longer than two weeks.  That'd be fine with me.  I did briefly think about driving out to take it back, but ended up feeling like sitting on the couch being lazy was a better option.

Final interesting thing:  Julie asked a question about truncating values in a vector that I couldn't answer because I don't know pandas.  "How would I do that in R?" I thought.  So I came up with a solution.  Then a better solution.  Then a better solution.  Then a solution that is basically the same, just slightly lazier.

# Read example data
df = data.frame(read.csv("/tmp/x.csv",header=TRUE))
df
  title values otherstuff
1     A      1        0.6
2     B      2        0.5
3     C     33        0.4
4     D      4        0.3
5     E      5        0.2
6     F     79        0.1

# Use which in index mode to find values above some threshold and assign a new value. Replace accepts the list of indices that match the query, and replaces the values they point to.
df$values2 = replace(df$values,which(df$values > 5, arr.ind = TRUE), 5) ; df
  title values otherstuff values2
1     A      1        0.6       1
2     B      2        0.5       2
3     C     33        0.4       5
4     D      4        0.3       4
5     E      5        0.2       5
6     F     79        0.1       5

# Define a logical selection above the threshold.  Replace accepts a list of Booleans to indicate which values to replace.
L = df$values > 5
df$values3 = replace(df$values,L,5); df
  title values otherstuff values2 values3
1     A      1        0.6       1       1
2     B      2        0.5       2       2
3     C     33        0.4       5       5
4     D      4        0.3       4       4
5     E      5        0.2       5       5
6     F     79        0.1       5       5

# Do that all at once, avoiding using an intermediate vector.
df$values4 = replace(df$values,df$values > 5,5); df
  title values otherstuff values2 values3 values4
1     A      1        0.6       1       1       1
2     B      2        0.5       2       2       2
3     C     33        0.4       5       5       5
4     D      4        0.3       4       4       4
5     E      5        0.2       5       5       5
6     F     79        0.1       5       5       5

# Do it with a with, why not, which is functionally the same, just avoids $ characters.
df$values5 = with(df,replace(values,values > 5,5)); df
  title values otherstuff values2 values3 values4 values5
1     A      1        0.6       1       1       1       1
2     B      2        0.5       2       2       2       2
3     C     33        0.4       5       5       5       5
4     D      4        0.3       4       4       4       4
5     E      5        0.2       5       5       5       5
6     F     79        0.1       5       5       5       5

I keep overthinking things in R.

• I don't know.  It seems like a good idea, but land supply in high population areas is an issue.  This would work best in low-density/low-cost-land areas, but that seems like this has the potential to add to blight if the region loses population.  Definitely a good idea for places with land and homelessness.
• Animals and flowers.
• I've been reading this series as part of my "take advantage of Amazon sale"/"accidentally look like I'm committing fraud" thing.  It's very good, but I wish they hadn't split the A/B story issues into two trades.  I get that they did it so each book would have a "complete" story, but I feel like I'm missing out on the back and forth between the two.

Thursday: Writing never feels productive.

Probably at least partially because I hate everything I write.  Oh well.  There's always tomorrow.  And maybe the weekend.  And part of next week.

Today's Pokemon was Mimikyu.

Twelve.

I'm sure this is equally easy elsewhere, but I really enjoyed that R is so simple:

zpt <- data.frame="" header="TRUE))</font" read.csv="" zpt.dat="">

summary(zpt)

levels(zpt$filter)

mycolors = c('yellow','green','blue','red','purple','violet')

with(zpt,plot(as.Date(dateobs),zpt_obs,col=mycolors[filter]))legend(as.Date("2016-08-01"), 27.25, legend=levels(zpt$filter), col=mycolors, pch=19)

dev.copy(png,"./date.png"); dev.off()

I'm sure there's a way to tunnel a mysql connection through ssh so I could read the data directly from the database, but I figured that would take too much time to sort out.  Or, I'd google it now, and there'd be literally piles of pages saying "yeah, it works basically exactly how you expect."  Most of R seems to be that way.

• I guess I'm buying a new book at some point.
• I feel old.
• Apparently the Wonder Woman movie is good.  That makes me so happy.
• It's good to hear things like this.
• And see things like this.
• I wish more comics understood this.
• And wrote her like this.
• "One Curlyboy with no friend."
• Watch out, kitty!