Distribution of Actual Births Around Estimated Due Date

Given that we're now in overtime here for our second child, I have become very interested in the question of what our "due date" meant in the first place.

The plotted on the left are data from a study of Canadian births, 1972-1986 (Arbuckle & Sherman 1989). I then started trying to aggregate some data for my own plots.

The first thing I noticed required a little attention was the "fence post" problem relating to recording gestational age. If a study A records births in the 39th week, where exactly is that on the x axis? Well, if we assume that they started counting with 1, then the 39th week is actually 38.5 +/- 0.5 weeks (they are counting fence). But if study B records births from week 39-40, the implication is that they started at 0 (they are counting fence posts), so 39-40 is 39.5 +/- 0.5. That was a little tricky.

Then we have studies that bin over different time intervals. If another study C records births at 37-41 weeks, how do we relate that to A and B? The intelligent way to plot this would be to use the probability density of going into labor, which divides out by the length of time over which the observation was made. So we should be careful to do that.

And then there's just the general problem that a lot of studies list percentages for births in each time bin, but don't list total populations or error bars, so we don't know what the errors in their measurements were. Ugh. So I crossed my fingers and hoped they followed good practices with their significant figures: I assigned an error equal to the last significant digit they list.


I got most of my data from this semi-thorough compilation of census data and going to some of the original sources. The data aren't great, but they are adequate (see plot to left). I fit to the data an increasing exponential tail to the left, plus a normal distribution. The fit isn't great (despite some claims in the literature of it being normally distributed), but it captures enough of the overall distribution.


The width of the normal distribution was 1.5 weeks, centered at 39.5 weeks. This seems consistent with several sources that suggest ~10% of pregnancies would go into the 42nd week if they were allowed to. It also suggests that the "due date" means the mean of the normal portion of the distribution. Yet fully 1/2 of pregnancies will go beyond the expected due date, and 1/6th will go past 41 weeks, according to this coarse fit.

Of course, none of this accounts for biases that we know exist. The growing prevalence of inductions and C-sections move births earlier artificially. Although the statistical significance may questionable owing to systematic biases (self-selection for uncomplicated pregnancies, etc.), it appears that the recorded midwife births go later than the aggregate (presumably hospital-dominated) births at the 2-sigma level. This may potentially indicate that without intervention biases, the distribution of birth dates around the expected due date could be broader and weighted toward later dates.

Who Dominates Health Care Costs?

There's nothing like a little bout with MRSA to make one pay a little more attention to the state of health care legislation. Two nights in the ER are definitely making me thankful for health insurance. Knowing that it wasn't going to cost me an arm and a leg to get antibiotics through an IV (in fact, it probably saved me the leg) definitely helped me to seek care early, rather than waiting for the infection to get truly life-threatening. And that probably saved in health care costs in the long run.

I've heard it argued many times by the other side that universal health care will drive up the cost of health care for everyone, because so-called "healthy people" will be paying, through their premiums, for the bills of the "unhealthy". Ignoring that:

1. the above is a tautological statement about what insurance is
2. people routinely go from the "healthy" group to the "unhealthy" group and back again
3. we should maybe feel a moral obligation to care for the unhealthy

Yeah, ignoring that, I wanted to know if the underlying assumption was, in fact, true. Who dominates health care costs? Is it the small number of extremely sick people? Or is it the larger number of moderately sick people?

To answer that question, I went searching for the population distribution of health care costs. I found the following publication: Variations in Lifetime Healthcare Costs across a Population (Forget et al. 2008). To the left are reproduced Figs. 4 and 5.

Given all the hype, I was somewhat underwhelmed to see that these curves depict (with the exception of an excess at the lowest cost bin) a gamma distribution. This isn't surprising, because a gamma distribution is supposed to represent the sum of a bunch of exponentially-distributed random variables.

To find the contribution of people in each cost bin to the total health care cost of the population, we simply need to multiply the population of that bin (drawn from a gamma function) by the mean health care cost of that bin (a linearly increasing function). Setting the mode of the gamma distributon to $90k for females, and tweaking the k and theta parameters (I'll chi-by-eye it at k=4.5, theta=1.0) we get the following distributions of fractional population (black) and fractional total health care cost (red), as a function of lifetime healthcare cost:

So who dominates health care costs? Those just slightly above the mode, which is to say, the large number of people who are just a little sicker than most. And that really could be any of us, folks.

Graph-SLAM

After a trying, but ultimately successful month spent extracting the family from Puerto Rico and re-embedding us in Berkeley, I'm just starting to get back on top of things enough to think about posting...

I had lunch yesterday with a good friend of mine, Pierre, who co-founded a company that specializes in sensory and mapping systems such as those that are used to create Google's "Street View". I was impressed to learn about their system for combining data from GPS, LIDAR, car odometers, and IMUs to create a consistent picture of how a vehicle is located and oriented in space as a function of time. They've spent a lot of time calibrating their systems, and use some sophisticated MCMC post-processing methods for deriving the actual trajectory of a vehicle.

Although the antennas in the PAPER array (that's the low-frequency interferometer I'm working on), are much less mobile than a car, there was considerable overlap between the problem Pierre has been working to solve and the calibration problem I am facing the requires positioning antennas and celesital sources as a function of time in the face of ionospheric distortion, variable gains, etc. Pierre pointed me to Graph-SLAM as a formal description of the problem that we are trying to solve, and suggested that Kalman Filtering with RTS Smoothing was a powerful technique for converging to the optimal solution (with covariance information) in linear time.

The Voynich Manuscript: Bootstrapping Language

The internet is a powerful and dangerous thing. It all started when I read the latest xkcd, which warned me that visiting cracked.com was dangerous. Then I read about the "6 phenomena that science can't explain" (which was a very dramatic title for some underwhelming mysteries), and next thing I knew, I was reading about the Voynich Manuscript. I learned about cryptography, glossolalia, the Manchu language, among other things. Then I took a look at the manuscript and before I knew it, I had a transcribed version of the manuscript in electronic form using the European Voynich Alphabet. And it just went downhill from there.

To summarize, the Voynich Manuscript (hereafter VMS) is a handwritten text with some illustrations some 500 years old. It uses glyphs no one knows how to read, it is not clear if it corresponds to a known language, it may or may not be encrypted, and little progress has been made in deciphering any of it, despite the fact that some bright people have tried. So I decided to have a crack at it.

The reason I got interested is because of the similaries to SETI. Arecibo, back in 1974, transmitted a message off into space that had been designed to be decrypted. We might receive a message like that some day. Or we might intercept something much like the VMS--a bunch of data in a language that we have no prior knowledge of--and we may be finding ourself trying to figure out how to bootstrap a language. That is to say, to learn the grammar and semantics of a language from a static example, without outside help.

Is this possible? For grammar, I'm pretty sure of it. I can imagine an algorithm (maybe Maximum Entropy Modeling and Bayesian learning applied to grouping and parsing) that uses correlations in the appearances of language elements (starting with letters and building up) and correlations in the behaviors of these elements relative to one another to build a model for parsing a language. For the VMS, I used something similar to this (not the MEM and Bayesian part) to show that spaces, line breaks, and paragraph breaks show similar grouping correlations relative to other VMS letters, and so can probably be considered one grammatical element of whitespace. That's a pretty simple thing to deduce, but it was actually something I was worried about in getting started with the VMS.

Sematics is another issue. Once upon a time, I would have had an optomistic answer to bootstrapping sematics from text that was not written for that purpose. However, after watching my child mysteriously acquire language, illustrating how hard-wired the human brain is for learning language from another human, and how much it relies on shared experience and feedback, I'm less sure.

I would be interested to know if there is a field of mathematics that studies sematics and the properties that a self-contained system needs to have to be able to generate sematical relationships. The Arecibo message relied on a shared physical environment to try to bootstrap sematics. I wonder if it would be enough to describe the rules of the grammar of a language in the language itself. That one, once the reader had deduced the relationships between elements, you would have a shared knowledge of that subject that might enable a reader to correlate the structure of the descriptions with the grammatical structure and thereby establish the first sematical relationships.

Anyway, after preliminary analysis of the VMS, I'm pretty sure that it's not random gibberish (there are correlations between elements on levels ranging from letters to words), and if it's encrypted, it's a weak form of encryption that preserves these correlations. My pet theory, extended from the glossalalia idea, is that this is actually plaintext in a natural language with an invented set of symbols, but that the natural language might be the accidental or intentional creation of a savant or scholar.

The Need for Speed

On the drive from San Juan to Arecibo this morning, I got to wondering about where my average driving speed fell in the distribution of drivers here in Puerto Rico. In the states, I felt like I was a pretty average driver, but here en la isla, the distribution of driving speeds is different. There are a lot of fast drivers, too be sure, but there is also a subpopulation of drivers whose speed is significantly (~10 mph) below the speed limit. This may be because relative to the US, PR is economically depressed and so more old cars are on the road, or as a reaction to the more erratic driving habits there seem to be here, but anyway, I definitely pass more people than pass me now.

So in an effort to discover where my driving speed fell relative to others (and in and effort to alleviate the boredom of driving 1.5 hrs alone), I started counting how many cars I passed and how many passed me as I was going 65 mph (the speed limit). Out of 55 pass events, only 9 involved me getting passed. To make this a tractable problem in my head, I decided to assume that driving speeds were normally (gaussian) distributed about a mean--even though this contradicts my anecdotal evidence above. Using this approximation, my first instinct was to say that 1/6 of the cars on the road were faster than me, and since ~2/3 of samples are within +/- 1 sigma of the mean in a gaussian distribution, 1/6 of the samples would be above +1 sigma. So I approximated that I was a 1 sigma driver.

But then it occurred to me that I needed to control for a significant sample bias. This is because the test I was doing wasn't randomly selecting cars and comparing my speed to them. Cars were far more likely to get selected if the difference between their speed and my speed was large. A car going the same speed as me would never pass me, and I would never pass it. But I would assuredly pass almost every car on the road that was going 10 mph as I went 65. The "road distance" that I sampled for different velocities is proportional to abs(v-v0), where v0 is my velocity. The effect this had on my samples was to underweight speeds close to my own and overweight the wings of the gaussian distribution I had assumed as my model. If I drove exactly the mean velocity, this effect would not be terribly important--if the model were correct, it would still be the case that as many cars passed me as I passed. But as my velocity moves away from the mean velocity, the "normal drivers" who are only going a little faster than me get undersampled, so I only see the drivers who are tearing around like a bat out of hell. At the lower end, I still see the real slow-pokes on the road, but I start seeing people who are going a bit faster than that, of which there are a lot more. The effect of this sample bias, it seems, would be to make it seem that I'm a farther outlier in my driving speed than I actually am.

So now, to figure out where I fall in the (normal) distribution of driving speeds, I need to know exactly what the mean driving speed and what sigma is, so that I can compensate for the abs(v-v0)
sampling factor. That means I need to figure out 2 numbers, but unfortunately, I only measured 1 number (that 1/6 of the passes while driving were me being passed) so I won't be able to properly constrain this problem. However, I should be able to figure out the mean on my drive home by finding the speed at which as many people pass me as I pass. For now, let's say this is 55 mph. Then all I need to do is find the sigma for which a gaussian distribution around 55 mph downweighted by abs(v-65 mph) has 1/6 of the area lying above 65 mph. I just solved that numerically on my computer, and it's saying that the best-fit sigma is ~22 mph. So that puts me at about +1/2 sigma. That seems reasonable.

An interesting next step (which I'm not going to do right now since I need to get to work) would be to translate the sample error in my pass measurements into an error in the determination of sigma, and then the error in my driving speed percentile.

What are the chances?

I cringe every time that I hear this phrase. I heard it most recently when my friend had her car stolen from San Juan. It was recovered in a semi-drivable state in Bayamon. She invested a couple thousand dollars to get rid of the "semi", only to have the car re-stolen a couple of months later. This time, when the car was recovered in Dorado, there was no "semi" to be had, so she's currently trying to sell it for pieces. While the police were fairly understanding (though less than helpful) the first time her car got stolen, the second time was occasion for all sorts of raised eyebrows and skepticism. And in exasperation my friend uttered the phrase in question.

"What are the chances" is a Pandora's box of bad statistics. Statistics is about hedging your bets given incomplete information, but this phrase is always uttered after the fact, when we have (relatively) complete information: it happened. So unless you plan on repeating the experiment, the chances are one. It happened.

As an example, let's take the famous Goat/Car Puzzle. There are 3 doors; one has a car behind it and the other two have goats. After you pick a door, the game-show host opens one of the other doors and reveals a goat. You are then offered the option of switching your choice to the other door. If you play this game repeatedly, you'll win more often if you switch your choice. But the instance just played out, the car was behind one of the doors, and if that was the door you picked, your chances of getting the car were 1. If you didn't, your chances were 0.

You might object: "What were the chances beforehand, when I didn't know where the goat and car were?". But to do that, you need to make some assumptions. You need to assume that at each playing of the game, the cars and goats are randomly assigned and/or you randomly pick doors. Otherwise, it might be the case that the car is always behind door #1, and you always pick door #2. Your chances of success wouldn't be so good in this case. You might have decent prior knowledge of how cars, goats, and doors are picked in this example, but for everyday occurrences, we usually have much more limited prior knowledge. Are cars randomly stolen, or are certain brands targeted? Are certain areas targeted? People often assume that these processes are random, but they rarely are. With limited priors on these events, the question "what are the chances" can't be answered with any certainty and any answers given should be taken with a great big shaker of salt.

Furthermore, people have selective attention. We ignore whole heaps of ordinary outcomes and only pay attention to ones that strike us as interesting. As a friend of mine once said: "Low probability stuff happens pretty regularly because stuff is happening all the time." Even if the processes involved are random, unlikely outcomes are to be expected if the processes are repeated often enough. People tend to ignore the ordinary outcomes, exclaim at the extraordinary ones, and then assume that something deeper is afoot. In my friend's case, the police started wondering if she was being personally targeted or if she was really bad at locking her car. But even if we assume a random model of car thefts, some number unlikely outcomes doesn't automatically imply that our random model is wrong.

Finally, we also need to keep in mind that in complex systems like real life, there may be a huge number of possible outcomes. But something has to happen. When you roll a die, each number only has a 1/6 chance of coming up. Would you roll a die once and then exclaim: "Wow, it came up six! What are the chances?" In real life, there might be millions of outcomes, each with one-in-a-million chance of coming true, but the fact that one of them happens shouldn't be surprising.

Fighting against all of the pitfalls inherent in asking "What are the chances?", I've developed a reflexive response: "What are the chances?" One.

Compressed Sensing and Wiener Filtering

Today I'm trying to expand my understanding of how we can best remove contaminant signals from the data we take with the Precision Array for Probing the Epoch of Reionization (PAPER). There is a specific problem I want to make sure we can solve for PAPER. Foregrounds to our signal, particularly synchrotron radiation, are expected to be very smooth with frequency. The idea put forth by the MWA and LOFAR groups is that by observing the same spatial harmonics at multiple frequencies, we should be able to remove such smooth components to suppress them relative to the cosmic reionization signal we are looking for. However, generating overlapping coverage of spatial harmonics as a function of frequency is expensive. My intuition is that since foregrounds do not have a spatial structure that changes dramatically with frequency, we shouldn't need to sample a given spatial harmonic very finely in frequency to get the suppression we want. This would allow us to spread our antennas out a little more and get measurements of the sky at a variety of spatial modes.

In many ways, our problem is analogous to what was done with the Cosmic Microwave Background (CMB). For foreground removal in CMB work, Tegmark and Efstathiou (1996) begin with an assumption that foregrounds can be described as the product of a spatial term and a spectral frequency term. This allows them to construct Wiener filters that use the internal degrees of freedom of their data, together with a model of their foreground and a weighting factor based on the noisiness of their data, to construct a filter for removing that foreground. For the most part, this is standard Wiener filtering, except they have to be careful about what they do to their power spectrum, so they apply a normalization factor to correct for a deficiency in Wiener filters. Tegmark (1998) goes on to generalize this technique for foregrounds that vary slowly with frequency. I'm in the process of wading through these papers, but they seem to be directly applicable to what we are doing, and seem to confirm my suspicions that synchrotron emission should be well-enough behaved to require only sparse frequency coverage of a wavemode in order to be suppressed.

Another tactic that I am investigating is that of compressed sensing which I was alerted to in talks by Scaife and Schwardt at the SKA Imaging Workshop in Socorro this last April. The landmark paper on this principle seems to be Donoho (2006), where it is shown that the compressibility of a signal (being sparse for some choice of coordinates) is a sufficient regularization criterion to faithfully reconstruct signals using a small number of samples. In a way, this technique has an element of Occam's Razor in it--it tries to find a solution, in some optimal basis, that needs the fewest non-zero numbers to agree with the measured data. At least, that's my take on it without having finished the paper.

The relevance of compressed sensing to image deconvolution is explored in Wiaux et al (2009), and it seems to be powerful. I'm excited by this deconvolution approach because it meshes well with the intuitive approach I've been taking to deconvolution, which was to use wavelets and a Markov Chain Monte Carlo optimizer to find the model with the fewest number of components that reproduces our data to within the noise. Compressed sensing seems to be exactly this idea, but is agnostic about the basis chosen, instead of mandating one like wavelets. Anyway, this technique may also be relevant to our foreground removal problem because we might be able to use it to construct the minimal foreground model implied by our data. For synchrotron emission, which should have smoothly varying spatial structure with frequency, I envision that this could construct a maximally smooth model that would allow us to use sparse frequency coverage to remove the foreground emission to the extent that it is possible to do so.

More on MCMC in Python

Markov-Chain Monte Carlo (MCMC) seems to be a promising technique for the calibration/imaging problem that we are facing with our experiment the Precision Array for Probing the Epoch of Reionization (PAPER). Yesterday, in addition to taking a crash-course in MCMC, I also started playing with PyMC, which implements, among other things, MCMC using Metropolis-Hastings chains. A first shot at a simple fitter using PyMC went something like this:

import pymc, numpy as n, pylab as p
from pymc import deterministic

x = n.arange(-10., 10, .01)

def actual_func(a, b, c): return a*x**2 + b*x + c

sig = .1
tau = 1/sig**2
noise = n.random.normal(0, scale=sig, size=x.size)
mdata = actual_func(1., 2., 3.)
mdata += noise

a_fit = pymc.Uniform('a', 0., 2.)
b_fit = pymc.Uniform('b', 1., 3.)
c_fit = pymc.Uniform('c', 2., 4.)
tau_fit = pymc.Uniform('tau', tau/3, 3*tau)

@deterministic
def func(a=a_fit, b=b_fit, c=c_fit): return actual_func(a, b, c)

mdata_var = pymc.Normal('mdata', mu=func, tau=tau_fit,
    value=mdata, observed=True)

mdl = pymc.Model([a_fit, b_fit, c_fit, d_fit, mdata_var, tau_fit])
mc = pymc.MCMC(mdl)
mc.sample(iter=1000,burn=250)  

a_hist, a_edges = n.histogram(a_fit.trace(), 40, (0,2))
a_bins = (a_edges[:-1] + a_edges[1:])/2
p.plot(a_bins, a_hist)

p.show()

I chose this example over the ones that came with PyMC because it was much closer to the kind of calibration problems we will be trying to solve with PAPER. An interesting point that came up early on was the reliance of MCMC on an estimate of noise levels in the data. I remember this from the Maximum Entropy Method (MEM) for deconvolution that I coded up in AIPY. An interesting technique, used here, is to actually include the noise level (tau_fit) as a variable that is fit for. This way you can specify a starting guess for noise, along with a confidence interval, and not have to pull a magic number out of a hat. In this example, the fitter does a good job of accurately determining noise levels. I think what happens in this code is that the current guess for the noise level is used as part of the calculation that determines the next state to jump to in the Markov chain, and that new state make include a revision of the noise level. This clearly might be instable for more complex systems, so I imagine some amount of care must be exercised in leaving noise as a free parameter.

Image Deconvolution with Markov-Chain Monte Carlo

I've decided to morph this blog to be more about short updates pertaining to current ideas I'm thinking about, rather than the long-winded philosophical rants I've posted so far. Hopefully this might keep me more engaged as a blogger and maybe even help me keep better track of the things I'm working on. Starting up this blog again, by the way, is a shameless procrastination technique, since my dissertation is due in about 1 month, and all my writing energy should really be focused on that...

After attending an SKA Imaging Workshop in Socorro, NM a couple of weeks ago, I've developed an interest in Bayesian statistics and Markov-Chain Monte Carlo (MCMC) techniques as they pertain to interferometric imaging. Having never taken a stats course, I'm scrambling a little to absorb the vocabulary I need to understand papers written on the subject. Fortunately, in this era of wikipedia, getting up to speed isn't that hard. After reading wiki articles on MCMC, Markov Chains, and the Metropolis-Hastings algorithm, I dived into EVLA Memo 102, which talks about a first shot at using MCMC for image deconvolution.

The Maximum Entropy Method (MEM) is a classic deconvolution technique (one I've already reimplemented for AIPY), but I'd like to go a bit further down this road. According to the standard implementation (which I gleaned from reading Cornwell & Evans (1984) and Sault (1990)) this algorithm uses a steepest descent minimization technique based on the assumption of a nearly diagonal pixel covariance matrix (i.e. the convolution kernel is approximately a single pixel). While this is an effective computation-saving assumption, I found that for the data I was working with, this assumption lead to the fit diverging when I started imaging at finer resolutions.

I think MCMC, by not taking the steepest decent, might be able to employ the diagonality assumption more robustly. I also think it's high time that deconvolution algorithms make better use of priors. The spatial uniformity prior in MEM makes it powerful for deconvolving extended emission, while the brightest-pixel selection technique in CLEAN makes it effective for deconvolving point sources. There's no reason we can't build a prior explicitly for a deconvolution algorithm that tells it to prefer single strong point sources over many weaker point sources, but also tells it that when all else is equal, entropy should be maximized.