Oct 15

Random samples in JS using R functions

For a JavaScript-based project I’m working on, I need to be able to sample from a variety of probability distributions. There are ways to call R from JavaScript, but they depend on the server running R. I can’t depend on that. I need a pure JS solution.

I found a handful of JS libraries that support sampling from distributions, but nothing that lets me use the R syntax I know and (mostly) love. Even more importantly, I would have to trust the quality of the sampling functions, or carefully read through each one and tweak as needed. So I decided to create my own JS library that:

  • Conforms to R function names and parameters – e.g. rnorm(50, 0, 1)
  • Uses the best entropy available to simulate randomness
  • Includes some non-standard distributions that I’ve been using (more on this below)

I’ve made this library public at Github and npm.

Not a JS developer? Just want to play with the library? I’ve setup a test page here.

Please keep in mind that this library is still in its infancy. I’d highly recommend you do your own testing on the output of any distribution you use. And of course let me know if you notice any issues.

In terms of additional distributions, these are marked “experimental” in the source code. They include the unreliable friend and its discrete cousin the FML, a frighteningly thick-tailed distribution I’ve been using to model processes that may never terminate.

Dec 14

Can pregnant women intuit the sex of their children?

“So let’s start with the fact that the study had only 100 people, which isn’t nearly enough to be able to make any determinations like this. That’s very small power. Secondly, it was already split into two groups, and the two groups by the way have absolutely zero scientific basis. There is no theory that says that if I want a girl or if I want a boy I’m going to be better able at determining whether my baby is in fact a girl or a boy.”

– Maria Konnikova, speaking on Mike Pesca’s podcast, The Gist.

Shown at top, above the quote by Konnikova, is a simulation of the study in question, under the assumption that the results were completely random (the null hypothesis). As usual, you’ll find my code in R at the bottom. The actual group of interest had just 48 women. Of those, 34 correctly guessed the sex of their gestating babies. The probability that you’d get such an extreme result by chance alone is represented by the light green tails. To be conservative, I’m making this a two-tailed test, and considering the areas of interest to be either that the women were very right, or very wrong.

The “power” Konnikova is referring to is the “power of the test.” Detecting small effects requires a large sample, detecting larger effects can be done with a much smaller sample. In general, the larger your sample size, the more power you have. If you want to understand the relationship between power and effect size, I’d recommend this lovely video on the power of the test.

As it turns out, Konnikova’s claims notwithstanding, study authors Victor Shamas and Amanda Dawson had plenty of power to detect what turns out to be a very large effect. Adding together the two green areas in the tails, their study has a p-value of about 0.005. This a full order of magnitude beyond the generally used threshold for statistical significance. Their study found strong evidence that women can guess the sex of their babies-to-be.

Is this finding really as strong as it seems? Perhaps the authors made some mistake in how they setup the experiment, or in how they analyzed the results.

Since apparently Konnikova failed not only to do statistical analysis, but also basic journalism, I decided to clean up on that front as well. I emailed Dr. Victor Shamas to ask how the study was performed. Taking his description at face value, it appears that the particular split of women into categories was based into the study design; this wasn’t a case of “p-value hacking”, as Konnikova claimed later on in the podcast.

Konnikova misses the entire point of this spit, which she says has “absolutely zero scientific basis.” The lack of an existing scientific framework to assimilate the results of the study is meaningless, since the point of the study was to provide evidence (or not) that that our scientific understanding lags behind what woman seem to intuitively know.

More broadly, the existence of causal relationships does not depend in any way on our ability to understand or describe (model) them, or on whether we happen to have an existing scientific framework to fit them in. I used to see this kind of insistence on having a known mechanism as a dumb argument made by smart people,  but I’m coming to see it in a much darker light. The more I learn about the history of science, the more clear it becomes that the primary impediment to the advancement of science isn’t the existence of rubes, it’s the supposedly smart, putatively scientific people who are unwilling to consider evidence that contradicts their worldview, their authority, or their self-image. We see this pattern over and over, perhaps most tragically in the unwillingness of doctors to wash their hands until germ theory was developed, despite evidence that hand washing led to a massive reduction in patient mortality when assisting with births or performing operations.

Despite the strength of Shamas and Dawson’s findings, I wouldn’t view their study as conclusive evidence of the ability to “intuit” the sex of your baby. Perhaps their findings were a fluke, perhaps some hidden factor corrupted the results (did the women get secret ultrasounds on the sly?). Like any reasonable scientist, Shamas wants to do another study to replicate the findings, and told me that has a specific follow-up in mind.

Code in R:

trials = 100000
results = rep(0,trials)
for(i in 1:trials) {
	results[i] = sum(sample(c(0,1),48,replace=T))
extremes = length(results[results<=14]) + length(results[results>=34]) 
dat <- data.frame( x=results, above=((results <= 14) | (results >= 34)))

Sep 14

Labor day distribution fun

Pinned, entropy augmented, digitally normal distribution, of no particular work-related use and thus perfectly suitable for today. Code in R:

iters = 1000
sd = 2
precision = 20
results = rep(0,iters)
for(i in 1:iters) {
	x = floor(rnorm(20,5,sd) %% 10)
	results[i] = paste(c('.',x),sep="",collapse="")
results = as.numeric(results)

Sep 13

The week in stats (Sept. 9th edition)

Bayesian Evolution

Jul 13

The hat trick

In his book Quantum Computing Since Democritus, Scott Aaronson poses the following question:

Suppose that you’re at a party where every guest is given a hat as they walk in. Each hat has either a pineapple or a watermelon on top, picked at random with equal probability. The guests don’t get to see the fruit on their own hats, but they can see all of the other hats. At no point in the evening can they communicate about what’s on their heads. At midnight, each person predicts the fruit on their own hat, simultaneously. If more than 50% of the guests get the correct answer, they’re given new Tesla cars. If less than 50% of the guests get it right, they’re given anxious goats to take care of. What strategy (if any) can they use to maximize their chances of winning the cars?

Answer: there is no strategy that works.

Kidding! Of course there’s a strategy, as you can tell by the length of this post. Did you come up with any ideas? At first glance, it seems like the problem has no solution. If you can’t communicate with the other party goers, how can you find out any information about the fruit on your own head? Since each person was independently given a pineapple or a watermelon with equal probability, what they have on their heads tells you nothing about what you have on your head, right?

My own initial strategy, after considerable (but not enough!) thought, was to bet on regression to the mean. Suppose you see 7 pineapples and 2 watermelons. The process of handing out hats is more likely to generate a pineapple/watermelon ratio of 7 to 3 than 8 to 2 (it’s most likely to generate an equal number of each type, with every step away from a 5/5 ratio less and less likely). Thus, I figured it would be best to vote that my own hat moved the group closer towards the mean. Following my strategy, we all ended up with goats. What did I do wrong?

The key to solving this problem is to realize that the initial process for handing out hats is irrelevant. All that matters is that, from the perspective of a given person, they are a random sampling of 1 from a distribution that is known to have either 7 pineapples and 3 watermelons, or 8 pineapples and 2 watermelons. Thus, each person knows that the probability a randomly sampled guest will have a pineapple on their head is somewhere between 70% and 80%. More precisely, it’s either 70% or 80%. In any case, so long as every person votes for themselves being in the majority, then the majority of guests will be voting that they are in the majority.

I simulated this strategy using parties of different sizes, all of them odd (to avoid the issue of having and equal number of each hat type). Here’s the plot, with each point representing the mean winning percentage with 500 trials for each group size. As always, you can find my code at the end of the post.

As you can tell from the chart, once we have 11 or more guests, it’s highly likely that we all win Teslas.

One way to look at this problem is through the lens of the anthropic principle. That is, we need to take into account how what we observe gives us information about ourselves, irrespective of the original process that made each of our hats what they are. What matters is that from the perspective of each party goer, their view comprises a random sampling from the particular, finite distribution of pineapples and watermelons that was set in stone once everyone had entered the room. In other words, even if the original probably of getting a pineapple was 99%, if you see more watermelons than pineapples, that’s what you should vote for.

This problem, by the way, is related to Condorcet’s Jury Theory (featured on the most recent episode of Erik Seligman’s Math Mutation podcast). Condorcet showed, using the properties of the binomial distribution, that if each juror has a better than 50% chance of voting in accordance with the true nature of the defendant, then the more jurors you add, the more likely the majority vote will be correct. And vice versa. Condorcet assumed independence, which we don’t have because our strategy ensures that every person will vote the same way, so long as the difference between types of hats is more than 2.

# Code by Matt Asher for StatisticsBlog.com
# Feel free to modify and redistribute, but please keep this header
iters = 500
numbPeople = seq(1, 41, 2)
wins = rep(0, length(numbPeople))
cntr = 1
for(n in numbPeople) {
	for(i in 1:iters) {
		goodGuesses = 0
		hats = sample(c(-1,1), n, replace = T)
		disc = sum(hats)
		for(h in 1:n) {
			personHas = hats[h]
			# Cast a vote based on what this person sees
			personSees = disc - personHas
			# In case of a tie, the person chooses randomly.
			if(personSees == 0) {
				personSees = sample(c(-1,1),1)
			personBelievesHeHas = sign(personSees)
			if(personBelievesHeHas == personHas) {
				goodGuesses = goodGuesses + 1
		if(goodGuesses > .5) {
			# We win the cars, wooo-hooo!
			wins[cntr] = wins[cntr] + 1
	cntr = cntr + 1
winningPercents = wins/iters
plot(numbPeople, winningPercents, col="blue", pch=20, xlab="Number of people", ylab="Probability that the majority votes correctly")