Kernel Density Estimators

Since high school science class, I've been making graphs that show one variable as a function of another - force as a function of extension, redshift as a function of distance, intensity as a function of wavelength, et cetera. But around that time I was also writing code - probably generating some kind of chaos or fractal - that needed to plot the distribution of data points in one dimension. Seems simple, right? Well, the standard solution is to produce a histogram - divide up the range into "bins", count how many fall into each bin, and plot that. But choosing the right number of bins is something of an art: too many and your plot is hopelessly spiky, too few and you can't see any features your data might have. So I wondered if there was something better. It turns out that there is (though it's not always clear when it actually is better); a tool called a kernel density estimator.



First, the problems with histograms. I generated a collection of 400(*) photon arrival phases corresponding to two equal peaks. Above is plotted a histogram of their arrival times. Below is a histogram of their arrival times shifted by half a bin. To me, at least, it's not at all obvious that the two histograms are of the same distribution, let alone the same sample.



Since the signal is periodic, one natural thing to try is to work with the Fourier series. It's not too hard to construct some coefficients of the Fourier series of the photon arrival times. If I used all (infinitely many) coefficients I'd get a collection of spikes, one at each photon arrival time, which isn't too useful. But if I discard all Fourier coefficients after a certain point, I smooth out those spikes into a more reasonable shape. Here's what I get using ten harmonics:



The two peaks look more similar now, which makes sense since the Fourier coefficients don't really care about the starting phase. There are those nasty ripples, though - in fact, they go negative, which is rather distressing for something that's supposed to be a probability distribution. (Incidentally, there's some Fourier magic that goes into that error band on the sinusoid but I don't want to get into it right now.)

One way to get rid of those ripples is to think of this as a problem in digital filtering, where they are an example of the "ringing" that can occur in digital filters with too-sharp cutoffs. In that context, the usual solution is to taper the coefficients to zero smoothly, and that's what I will do. But there is a statistical way to think about what's going on.

First of all, tapering or cutting off the Fourier coefficients can be thought of as multiplying the Fourier coefficients by another set of Fourier coefficients. By a wonderful theorem of Fourier analysis, this amounts to convolving the set of photon arrival times by a kernel. That is, we take the forest of delta-functions representing the photon arrival times, and replace each delta function with a copy of some smoother function, and add them all together. This process, when carried out on the real line, is known to statisticians as constructing a kernel density estimator (the "kernel" is the smoother function, and it is used to estimate a probability density).

Simply chopping off the Fourier coefficients corresponds to the kernel sin((n+1/2)x)/sin(x/2), the periodic "sinc" function. This has a peak at zero but oscillates from positive to negative, so it's not too surprising that we got ripples. So to get rid of the ripples, we just choose a kernel that is positive and whose Fourier coefficients we know. There is a natural (if somewhat obscure) candidate: the von Mises probability distribution. This is a circular analog of a Gaussian, both in appearance and in a technical sense: the Gaussian distribution is the distribution with maximum entropy for specified mean and standard deviation. The von Mises distribution has maximum entropy for specified "circular mean" (which includes information about spread as well as location). In any case, it's a positive smooth periodic distribution whose Fourier coefficients can be computed in terms of Bessel functions using your favourite scientific computing tool. So using it instead gives this:



This looks pretty good; it shows the distribution pretty cleanly, with no ringing. There's only one fly in the ointment: while I don't need to specify a starting phase, I do need to come up with a parameter - analogous to the number of harmonics - that determines the width of the kernel. If I choose too wide a kernel, it flattens out all the features in my data; if I choose too narrow a kernel I get something horrible like this:



So in a way I haven't solved the original problem with histograms that motivated me. Fortunately, one of the selling points of kernel density estimators is that the statistical literature is full of papers on how programs can automatically choose the degree of smoothing. None of their techniques (that I have found) are appropriate for periodic data, but I have my own ideas about that (to be written up in future).

(*) 400 is not actually a very small number of photons; Fermi often gets something like one photon a week from a source.

Ghetto apple crumble

I made a surprisingly successful apple crumble today. Apple crumble, as you may know, is basically apples with sugar and cinnamon, with a topping of butter, flour, oatmeal, and sugar, all in about equal quantities. Almost idiot-proof. So why "surprisingly successful"? Well, my oven died halfway through (in fact after preheating but before the crumble went in). It turns out you can microwave an apple crumble and it's fine. Who knew?

Recipe below.

Filling:

• Five apples, peeled and diced. Granny Smith are good for a bit of tartness to balance the sweet crumble. (Peeling is technically optional but the peel makes for an awkward texture.)
• 5 mL cinnamon.
• 100 mL sugar.
• a little butter.

Crumble:

• 125 mL brown sugar (white will do fine).
• 125 mL rolled oats.
• 125 mL flour.
• 100 mL butter.
• a mL or two of vanilla.
  

Put the apples in a deep baking dish, sprinkle the sugar and cinnamon on top and shake a bit; scrape a little butter on top.

In a bowl, add the crumble ingredients. In the likely case that the butter is hard, microwave it for a few seconds to soften it up. Mix the ingredients together; you should get a crumbly mass. Pack it in on top of the filling.

Microwave on 50% for fifteen to twenty minutes; when it's ready the top layer will have sunk somewhat because the apples will have softened. As it cools the crumble on top will harden somewhat.

Serve hot or cold. Particularly good with a scoop of vanilla ice cream.

Alchemy

The alchemists' dream (one of them, anyway) was always to make gold. We now know there are very good reasons they couldn't: since gold is an element, making it from anything that doesn't contain gold(*) requires you to change the nuclei of the atoms involved, while all the alchemists had access to was the electron shells around the atom(**). So their efforts were basically hopeless. Now, though, we do have the ability to manipulate nuclei, and in fact we do so on industrial scales. So could we make gold? In fact, let's be more ambitious: could I make gold in my basement? The answer is, surprisingly, yes.


First of all, the easiest nuclear reaction to go for is to transmute mercury into gold. Mercury is commercially available (though somewhat encumbered by very sensible environmental concerns) and not actually very expensive - on the order of $18 per kg (***). Gold, by comparison is more like $40000 per kg. So there's room for some profit here.

How could I make the nuclear reaction happen? This particular reaction needs so-called "fast neutrons", that is, neutrons that are still zipping around at the high energies typical of nuclear reactions, as opposed to neutrons that are bouncing around at energies consistent with room temperature. I could stick the gold in a "fast breeder reactor", but I don't actually have one in my basement, and they're kind of hard to build. I could use a particle accelerator to generate some neutrons (basically by bashing nuclei around until some neutrons fall off) but while I do have a particle accelerator in my basement, it takes one a lot more serious than I can reasonably build to get neutrons out. Nuclear fusion reactions give out neutrons, though, so all I'd have to do would be to build a fusion reactor in my basement. Improbable as it sounds, this actually is feasible, provided I'm not trying to get any energy out.

The trick is that there's a fusion reactor, called the Farnsworth-Hirsch fusor, that is surprisingly simple to build. It is actually something of a cross between a fusion reactor and a particle accelerator: I'd set up an electrical potential in a spherical configuration, accelerating deuterium ions towards the center, where they'd crash into other deuterium ions, occasionally hard enough for fusion to happen. This fusion would release a fast neutron.

To make gold, then, all I'd have to do would be to build a fusor, surround it with a blanket of mercury, run it for a while, and then extract the gold from the mercury. Simple, really.

Let's look at the economics, though.

Suppose we want to make a kilogram of gold, giving us $40000. We need about a kilogram of mercury, costing $18. We also need about 5 g of deuterium (assuming perfect efficiency), which would cost about $30. Finally, we need the power to run the fusor. That's not going to be cheap. An optimistic number for the best fusor ever built is about 10^12 neutrons per second from about 4 kW input. That amounts to 9*10^14 neutrons per kilowatt-hour. Assuming perfect efficiency again, we need about 3*10^24 neutrons for our kilogram of gold, or 3 terawatt-hours, about the world's total energy production for an hour and a half. At $0.10 per kilowatt hour (I live in the land of cheap hydroelectricity) that's three hundred billion dollars.

There's a somewhat more disturbing possibility, though. Gold is easily obtained; you can just buy it. But as a global society, we try very hard to make sure you can't easily get plutonium, particularly plutonium-239, since that is well-suited to building atomic bombs. (You can make bombs out of uranium too, but that requires you to separate the different isotopes, which are very nearly chemically identical. Plutonium, on the other hand, can easily be separated from uranium since it is a different element.) Uranium isn't too hard to come by, especially "depleted uranium" (uranium with most of the uranium-235 removed) - armies fire the stuff at each other, for example. And if you had lots of U-238, a fusor would let you make plutonium out of it. The cost would be high, hopefully prohibitively so, but you could do it without doing anything that would put you on the radar of the IAEA. Fortunately, the power use is so outrageous we don't really need to worry about it.

So, in short, I could make gold in my basement, but not any appreciable quantity, and not for any kind of sensible price.


(*) Since gold is a noble metal, there aren't many chemical compounds that contain gold; unlike, say, iron, gold is often found on Earth as lumps of raw gold. So while in principle alchemists could have started from some gold compound and gotten the original gold back, this would not have been a very interesting accomplishment.

(**) There are actually situations where you can affect the nucleus by manipulating the electron shells. For example, if an isotope decays by electron capture, you can drastically slow down its decay by stripping away all its electrons. But stripping away all the electrons from a reasonably heavy element is one of those things that's virtually impossible under terrestrial conditions but not too rare astrophysically. In any case this has no effect on stable isotopes.

(***) Canadian dollars and American dollars are equivalent to astronomical accuracy.

Curve Fitting part 4: Validating Bayesian Code

In my previous post, I wrote a tool to use Bayesian inference to ask whether a collection of photons represented a pulsed signal, and if so, what its parameters were. It gave plausible answers, but knowing my own fallibility, I really want some more careful test before trusting its output.

I previously talked about how to implement such a correctness test in a frequentist setting. But in a Bayesian setting, what one gets are probability distributions describing our knowledge about the hypotheses - how can a program test those?

After racking my brains for a while trying to find a way to do this, I asked about it on the scipy-user mailing list, and received a few useful suggestions, the most valuable of which was to send a polite email to Professor Tom Loredo. I did, and he replied immediately with a very helpful email, and some links to other people's work on the subject.

It turns out that, in the context of Bayesian fitting, the issue is called "calibration". The key realization is that if you select a model according to your prior distribution, generate a data set, and then do your fitting to generate posterior distributions, then your "true" parameters that you used to generate the data set look just as if they had been drawn from the posterior distribution.

This makes a certain kind of sense: after all, your posterior distribution is supposed to describe the distribution of models that might have generated your data set.

So how do you turn this into a test? After all, if you just get one sample from a distribution, it's pretty hard to say anything very definite, especially when the tails - the most unlikely extreme parameters - are not necessarily very well modeled. If you try to generate another sample, you pick a different model and must then generate a different data set, so you get a new point but a totally different posterior distribution. So now you have two points, drawn from two different distributions, and your situation has not necessarily improved.

There's a trick to let you work with a large sample: you transform them to all have the same distribution. You can do this because you know the posterior distribution you're working with; in my case I have its values sampled evenly. So I can construct the cumulative distribution function and use its inverse to get the posterior probability of a value less than the true model. This should be distributed uniformly no matter the shape of the posterior.

In fact, I'd rather use a slightly different trick: instead of looking at the probability of a value less than the true model, I'll use the probability of a value more extreme than the true model. Essentially I'll combine both tails. This has a more direct bearing on the question of confidence intervals, and still results in a uniform distribution as I repeat the process many times.

Once I have a long list of many values that should be distributed uniformly, there are any number of tests I can use; I'll use the handy but not particularly good Kolmogorov-Smirnov test. My unit testing code now looks like this:

def test_credible_interval():

    M = 100
    N = 50

    ps = []
    for i in range(M):
        f, p = np.random.rand(2)

        events = bayespf.generate(f, p, N)

        phases, fractions, r, P = bayespf.infer(events)

        frac_pdf = np.average(r,axis=0)
        fi = np.searchsorted(fractions, f)
        p = np.sum(frac_pdf[:fi])/np.sum(frac_pdf)
        p = 2*min(p, 1-p)

        ps.append(p)

    assert scipy.stats.kstest(ps,lambda x: x) > 0.01

Note that I don't bother testing the posterior for pulse phase. This is just laziness.

In a real Bayesian problem, there would usually be many more parameters, and I would probably not be able to evaluate the posterior on a grid like this. I'd probably be using some sort of Monte Carlo method, which would return a list of samples drawn from the posterior instead. But there are reasonable approaches in this setting too.


So far so good. But what about the probability that the signal is pulsed? In some sense hypothesis testing is just a special case of parameter fitting, but the discreteness of the "parameter" poses a problem: there's no way we can hope for a nice uniform distribution when only two values are possible. Professor Loredo very kindly sent me an unpublished note in which he addresses this problem, showing that if your code accepts a hypothesis whose posterior probability is P_crit, then a correct piece of code will be right - have chosen the correct hypothesis - with probability at least P_crit. Unfortunately P_crit is only a lower limit on the probability of getting things right; but there isn't really a way around this: suppose my code were working with a million photons a run. Then it would, in almost every case, give a probability very close to zero or one. There would be very, very few errors, no matter what value P_crit you set.

The fact that P_crit is only a lower limit does mean that this doesn't allow you to catch code that is too conservative: if your code makes fewer errors than its probabilities claim it should, this test has no way to tell that that's what's happening.

One must also choose P_crit carefully. In my example, if I choose P_crit=0.5, then code that simply chose randomly, or returned the same value all the time, would pass: after all, my priors specify equal probabilities for each model, and with P_crit=0.5 the code only needs to be right half the time. On the other hand, with P_crit really high, the code will almost never be wrong, although it only takes a few errors to fail, but it will almost never actually come to a conclusion, so you'll wait forever for evidence. There should be some fairly well-defined optimum value of P_crit, and I need to think more about what it is.

In any case, having selected P_crit and the other parameters, I can write a unit test for the hypothesis testing component:

def test_pulsed_probability():
    np.random.seed(0)
    p_accept = 0.75
    M = 50
    N = 30

    accepted = 0
    correct = 0
    total = 0
    while accepted<M:
        if np.random.rand()>0.5:
            f, p = np.random.rand(2)
            pulsed = True
        else:
            f, p = 0., 0.
            pulsed = False

        events = bayespf.generate(f, p, N)

        phases, fractions, r, P = bayespf.infer(events, n_phase=100, n_frac=101)

        if P>=p_accept:
            accepted += 1
            if pulsed:
                correct += 1
        if P<1-p_accept:
            accepted += 1
            if not pulsed:
                correct += 1
        total += 1
    assert 0.01<scipy.stats.binom(M,p_accept).cdf(correct)

As a final note, I'd like to thank Professor Loredo for his help, but note that any errors in the above code or description are entirely mine.

Curve Fitting part 3: Bayesian fitting

When you fit a curve to data, you would usually like to be able to use the result to make statements about the world, perhaps something like "there's a fifty percent chance the slope is between 1 and 2". But this is a bit peculiar from a philosophical point of view: if your data is a set of measurements of some real-world phenomenon, then it's a bit funny to talk about probabilities that the slope has certain values. The phenomenon has some fixed slope, so we can't somehow repeat the experiment many times and see how often the slope is between 1 and 2. But there is a way to make such a statement meaningful: Bayesian inference.

The basic idea is that you use probabilities to quantify the degree of uncertainty you have about the world; you are using a system of logic that uses probability and probability distributions rather than binary logic. This may sound fuzzy and esoteric, but Bayesian logic is used very successfully in, for example, automatic systems to evaluate whether email is spam.

When applied to data, Bayesian reasoning lets you make meaningful statements about probabilities of parameter values, at the cost of making some explicit assumptions going in, and also at the cost of some potentially substantial computation. I'll work through an example of fitting a light curve to a set of photon arrival times, using a Bayesian procedure.


First the problem setting: suppose we observe a pulsar, whose period we know exactly (perhaps from radio observations), with an X-ray or gamma-ray telescope. We see some (fairly small) number of photons, and we want to know whether the flux we see is modulated at the pulsar period. We "fold" the photon arrival times, recording the pulsar phase when each one arrives. So our data is a collection of some hundreds or thousands of numbers between zero and one.

The model we'll fit includes some fraction f of photons whose arrival times are distributed as a sinusoid with a peak at phase p; the remaining fraction (1-f) are independent of phase.

The key idea of Bayesian curve fitting is that if you have some collection of hypotheses Hi about the world, each having some probability P(Hi), and you make some observation, there's a simple procedure to update these probabilities to reflect your new knowledge.

From a philosophical point of view, it's a bit worrisome to have to supply hypotheses (called "priors") about what the world is like before we ever make any observations. It amounts to making assumptions in the absence of data. But in fact there are assumptions built into the usual "frequentist" methods of inference as well, and in Bayesian inference the ability be explicit about the hypotheses at least makes it clear what's going on.

What assumptions should we make for our pulsar? Well, there might or might not be pulsations, so we'll have two basic hypotheses: no pulsations and pulsations. Absent any information, we'll assume these are equally likely. Then, if there are pulsations, we need to specify prior distributions of phase and pulsed fraction. Since both these parameters are between zero and one, we'll just take a so-called "flat prior" that makes all values between zero and one equally likely.

Given these priors, we need to figure out how to use our observed photon arrival times to update the priors to give us "posteriors". The general formula is:

P(Hi|D) = P(D|Hi) P(Hi) / P(D)

That is, the probability of hypothesis i given the data, P(Hi|D), equals the probability of the data given Hi, P(D|Hi), times the prior probability of Hi, P(Hi), divided by the probability of the data given any hypothesis.

The first thing to note is that P(D|Hi) is just what we need to evaluate for a maximum-likelihood estimate: how likely data like what we observe is to arrive given some hypothesis. We only need to define it up to a constant, since it appears in both numerator and denominator. For our problem, the probability density for pulse arrival times is p(f,p,t) = f(1+cos(2 pi (t-p)))+(1-f). So P(D|Hi) is the product of p(f,p,ti) for all events ti.

How do we form P(D)? Well, since we have two hypotheses, H0 (no pulsations) and H1 (pulsations), we can write P(D) = P(D|H0)+P(D|H1). Further, H1 is actually a family of hypotheses depending on two parameters, so we need to integrate P(D|H1) over all possible values of the parameters.

If we apply the above formula, then, we should get two things: a posterior probability that there are any pulsations at all, and a posterior probability distribution for the two parameters.

Let's look at python code to implement this. First of all, we're going to need to be able to generate fake data sets:

def generate(fraction, phase, n):
 m = np.random.binomial(n, fraction)
 pulsed = np.random.rand(m)
 c = np.sin(2*np.pi*pulsed)>np.random.rand(m)
 pulsed[c] *= -1
 pulsed += 0.25+phase
 pulsed %= 1

 r = np.concatenate((pulsed, np.random.rand(n-m)))
 np.random.shuffle(r)
 return r

This routine generates the photons in two parts. First it decides randomly how many come from the pulsed component. Then the photons from the pulsed component are generated uniformly. To convert this to a sinusoidal distribution we select some of the photons in the lower part and move them to the upper part. We then add in some uniformly-distributed photons, and shuffle the two samples together.

Now we write a routine to evaluate the probability density function:

def pdf_data_given_model(fraction, phase, x):
 return fraction*(1+np.cos(2*np.pi*(x-phase)))+(1-fraction)

Note that in spite of appearances, this routine can act on an array of values at once; this is important since python's interpreted nature makes each line of python take quite a long time.

And now the fitting routine:

def infer(events, n_phase=200, n_frac=201):

 events = np.asarray(events)
 phases = np.linspace(0,1,n_phase,endpoint=False)
 fractions = np.linspace(0,1,n_frac)

 lpdf = np.zeros((n_phase, n_frac))
 for e in events:
     lpdf += np.log(pdf_data_given_model(fractions, phases[:,np.newaxis], e))

 # This weird-looking hack avoids exponentiating very large numbers
 mx = np.amax(lpdf)
 p = np.exp(lpdf - mx)/np.average(np.exp(lpdf-mx))

 S = np.average(np.exp(lpdf))

 return phases, fractions, p, (S/(S+1))

This uses one of the simplest approaches to calculating the distribution and its integral: just evaluate on a grid. Integration then becomes averaging. More sophisticated Bayesian problems usually involve high-dimensional integrals, and so a whole elaborate machinery has evolved for efficiently evaluating these (for example the python package pymc).

Finally, some wrappers to generate a fake data set, call the fitting routine, and plot and print the results:

if __name__=='__main__':
 import pylab as pl

 events = generate(0.2,0.5,200)
 phases, fractions, r, P = infer(events)
 print "Probability the signal is pulsed: %f" % P

 pl.subplot(211)
 pl.contourf(fractions, phases, r)
 pl.xlabel("Pulsed fraction")
 pl.ylabel("Phase")
 pl.xlim(0,1)
 pl.ylim(0,1)

 pl.subplot(212)
 p = np.average(r,axis=0)
 li, mi, ui = np.searchsorted(np.cumsum(p)/np.sum(p),
         [scipy.stats.norm.cdf(-1), 0.5, scipy.stats.norm.cdf(1)])
 pl.plot(fractions, p)

 pl.xlabel("Pulsed fraction")
 pl.ylabel("Probability")
 pl.axvline(fractions[li])
 pl.axvline(fractions[mi])
 pl.axvline(fractions[ui])
 print ("Pulsed fraction: %f [%f, %f]" %
         (fractions[mi], fractions[li], fractions[ui]))

 pl.show()

One key point here is that when I want to know the distribution of pulsed fraction but don't care about the phase, I integrate (i.e. average) the joint distribution along the phase direction.

This gives us the following plot:



And the following output:

Probability the signal is pulsed: 0.450240
Pulsed fraction: 0.210000 [0.100000, 0.315000]

So it looks like the fitting routine is working: even with relatively few photons and a small pulsed fraction, it has selected quite good best-fit values. The probability that the signal is actually pulsed seems a little low, but keep in mind that we have only two hundred photons, and only forty of these are actually pulsed (while a Poisson uncertainty on the number of photons would be something like 14). But giving plausible results is not really enough: I want to systematically test this routine for correctness. But that will be another post.

Testing statistical tests

Often one wants to ask something like "is there a periodic signal in this data?" or "is there a correlation between these two data sets?". Often there is some way to calculate a number that represents how much signal there is or how much correlation there is. But of course there is always noise and uncertainty in the data, and so it's possible that the apparent signal or correlation is actually just noise. So for such a number to be useful, it must also come with an estimate of how likely it would be to get such a value if there were no signal, only noise. Such an estimate is often called a p-value.

I should say, before I go on, that there are a number of things a p-value isn't: it's not (one minus) a probability that the signal is real. It's not a measure of the strength or importance of the signal. It's not even enough to tell you whether you should believe the signal is real: if you look at a hundred thousand data sets, you should expect to find many with p<0.01 even if there's nothing there. But this has been better discussed elsewhere, and it's not my topic today.

Today I'm addressing a specific issue: suppose that you have implemented some new way of quantifying the presence of periodic signals, or correlations, or whatever, and your implementation returns a p-value. How do you make sure that p-value is actually correct?


Of course, you can (and should!) go over your code carefully, with a critical eye. But evaluating the p-value often involves evaluating some arcane special function. How can you be sure that you implemented it correctly, or for that matter, that you got the math right in the first place? It turns out there's a nice way to test your code as a black box.

For our example test, let's take a well-known test: the chi-squared test. As I'm going to use it here, we have m measurements, each with measurement uncertainty one, and we want to know whether the underlying quantity is the same in all measurements. The chi-squared test works by constructing the quantity chi-squared, which is the sum of the squares of the differences of the measurements from the mean:

def chi_squared(values):
  return np.sum((values-np.mean(values))**2)

This quantity has a known distribution, naturally enough called the chi-squared distribution. It is parameterized by a number called the "degrees of freedom", which in our case is m-1, the number of measurements minus one because we subtracted the mean. For my p-value, I will just ask what the probability is that random noise would produce a chi-squared value larger than what we observed. (Note that we might well want to consider unusual any data set for which the chi-squared value is unnaturally small as well as ones where the chi-squared is unnaturally large, but that can only happen if we have overestimated the uncertainties on our measurements, so I'll leave it aside.)

def chi_squared_p(m,c):
  return scipy.stats.chi2(m-1).sf(c)

For a distribution in scipy.stats, the "sf" method evaluates the "survival function", that is, what fraction of the distribution's values are larger than the given value.

So now we have a test and a p-value calculation. How are we to check that we implemented it correctly? Well, let's pick an m and a p-value we're interested in, say p0 = 0.05. Let's also pick a number of repetitions, N. We will generate N fake data sets in which there is nothing but noise, and see how many times our statistic reports p<p0. This should be something close to p0*N. How close? Well, the statistic should behave like flipping an unfair coin with probability p0, so we can use the binomial distribution to find limits that should contain the number of false positives 98% of the time.

def test_p(N,p0,m):
  false_positives = 0
  for i in range(N):
      if chi_squared_p(m,chi_squared(np.random.randn(m)))<p0:
          false_positives += 1

  assert scipy.stats.binom(N,p0).cdf(false_positives)>0.005
  assert scipy.stats.binom(N,p0).sf(false_positives)>0.005

There we have a test; the more trials you let it run (N), the tighter constraints it puts on the p-value. Unfortunately, it fails 2% of the time even if everything's fine, and it's random, so if it fails, you can't restart running it in the debugger (since the next run will get different values). There's a way around these problems, too. The first problem can be avoided by running the test once, then if it fails, rerunning it. Then the test only reports a failure 0.04% of the time if all is well, but a real problem in the algorithm will probably still show up. The second problem you can solve by seeding the random number generator every time you run the test. In python, decorators provide a handy way to avoid having to write boilerplate code to do both these things for each test; I have two decorators, seed() and double_check, that do this. I'll omit their code here because they have some nastiness to work around limitations in nosetests (but you can find them, along with a more detailed example of the techniques in this post here).

This test is nice, but a little limited: it only tests one p value, p0. Since every time the statistical test runs it returns a p value, surely we can do better?

Well, since the p value is supposed to be the probability that a random data set will exceed the value obtained in a particular test, if we generate lots of noise samples, we should have the fraction whose value is less than p0 roughly equal to p0 for every p0 - in other words, the p-values returned should be uniformly distributed. So we can use a statistical test for uniform distribution to check whether they're coming out right! One such test, present in scipy, is the Kolmogorov-Smirnov test. That gives the following scheme:

def test_ks(N,m):
  pvalues = [chi_squared_p(m,chi_squared(np.random.randn(m))) for i in range(N)]

  D, p = scipy.stats.kstest(pvalues, lambda x: x)
  assert p>0.01

This does have the drawback that the K-S test is known to be most sensitive in the middle of the distribution, while what we care about is actually the tail. A minor modification can help, at the cost of some speed:

def test_ks_tail(N,m,p0):
  pvalues = []
  while len(pvalues)<N:
      pv = chi_squared_p(m,chi_squared(np.random.randn(m)))
      if pv<p0:
          pvalues.append(pv/p0)
  D, p = scipy.stats.kstest(pvalues, lambda x: x)
  assert p>0.01

This has the embarrassing problem that it's too good: it turns out that while this works for the chi-squared statistic I describe, for the K-S test itself, the p-value returned is actually rather approximate, so that this test tends to fail.

Finally, there's one essential thing to check: how stringent are these tests? If we return the wrong p-value, do they pick it up? A quick check can be done simply by modifying the chi-squred calculator to jigger the value a bit, then checking that the tests fail. With a thousand trials, the tests pass just fine if I return 1.01*chi-squared; I have to increase it to something like 1.05*chi-squared to start getting failures.

This brings us to the elephant in the room: the power of a statistical test. The p-value is really only half the answer; it tells us how likely we are to get an apparent signal when there's nothing there. It tells us nothing about whether we would actually see a signal if it were there; for that, you need a different quantity, the power of the test.

Cosmos

In honor of Carl Sagan day, I'd like to link to Cosmos:





Nearly thirty years later, and in spite of tremendous advances in astronomy, Cosmos is still a wonderful introduction. Rather than focus on the science, Carl Sagan does a great job of sharing the wonder of discovery.

I am kicking myself now: when I was much younger, I was visiting a cousin who works at Cornell, and he pointed to Carl Sagan's office and asked if I wanted to go up and say hello. Out of shyness, I guess, I demurred, but it would have been a fascinating visit.

Eureka?

The most exciting phrase to hear in science, the one that heralds new discoveries, is not 'Eureka!' (I found it!) but 'That's funny ...'

-Isaac Asimov

I recently found a very exciting new millisecond pulsar. But my first thought was not "Wow! A new millisecond pulsar!" but "Isn't that a suggestive bit of interference?"

To explain myself a bit further, I was looking at candidates from the drift-scan survey. These are all the periodic signals we picked up with the GBT, and they naturally include every cell phone call, car ignition, laptop computer, and worn-out electric blanket in the vicinity. Most are easy to distinguish from real pulsars, but some aren't. One characteristic pulsars generally have is they're noisy: after all, they're faint astronomical sources, so it would be very strange if they were so strong we couldn't see any noise in the observation. But in the plot showing 1023's signal, seen above, you can see, there's no noise. Turns out 1023 is just plain bright. I was not the first to misclassify it, though.

J1023 normally looks like a fairly ordinary 17th magnitude star. It's in star catalogs and photographic plates (enter RA 10 23 47.67 and Dec 00 38 41.2) going back as far as 1952, but it never got any attention in particular until 2002.

In 1998, the FIRST sky survey was carried out with the VLA. This is a rather different beast from a pulsar sky survey; it gives average brightnesses over several-minute integration times, so it's not going to have any luck detecting pulsations. But because the VLA is an interferometer, it's able to generate quite high-resolution images. This kind of survey, called a "continuum" sky survey (as opposed to a survey for spectral lines or a pulsar survey) is good for finding nebulae, radio stars, and galaxies. J1023, or to give it its full name, FIRST J102347.67+003841.2, caught the attention of people analyzing the FIRST data because it was variable: in the three observations just days apart, its radio brightness varied by a factor of at least three. Galaxies generally don't change so quickly, so they thought the source was worth following up.

Since the source had an optical counterpart, they chose to follow it up by taking optical spectra. These optical spectra, and those by another group, showed that it was blue and had double-peaked emission lines. A blue spectrum indicates very hot gas, emission lines indicate hot diffuse gas, and the fact that they were double-peaked is normally a sign that they were coming from an accretion disk: gas on one side of the disk is moving rapidly towards us, so its emission line is shifted towards the blue end of the spectrum, while gas on the other side is moving rapidly away, so its emission line is shifted towards the red end. Since the hottest gas is the nearest the center and the fastest moving, it produces the strongest emission, and you get a double-peaked spectrum. Other observers also looked at the source with high-speed photometry (just looking at how bright the object was) and found it was flickering, which is normal for accreting systems, as turbulent knots in the disk come and go. So people looked at all this data and classified J1023 as an unusual "cataclysmic variable", that is, an accreting white dwarf.

In 2002, observations showed that the emission lines were gone, the spectrum had gone back to the Sun-like colors observed in 1999 and earlier, and the flickering had tailed off. All that was left was a Sun-like star that varied in brightness by 0.4 magnitudes in a predictable fashion as it travelled around its orbit. This return to quiescence is normal for cataclysmic variables: they go through active phases and passive phases. But John Thorstensen and Eve Armstrong decided to try to come up with a model that explained the light curve (the brightnesses and colors).

When you have a binary system like this in which the stars are close together, the companion is heated by the white dwarf. When the hot side is facing us, the star is brighter and bluer than when the cool side is facing us. So it shouldn't be too hard to look at the light curve and figure out how bright the white dwarf is.

Well, as always, it's not as easy as it sounds, but by dint of great effort, Thorstensen and Armstrong managed to come up with a model that fit. But there was a problem: it needed a massive bright white dwarf. So massive, in fact, that it couldn't be a white dwarf, and so bright that we should have seen it in the spectrum, as a blue bump. There are many possible explanations for this sort of thing; software bugs, calibration problems, and so on. But they did a careful analysis, didn't find any of those, and decided to go out on a limb and suggest that the prevailing opinion was wrong: J1023 didn't contain a white dwarf at all, but a neutron star.

There the question stood for a few years. Homer et al. took X-ray observations of J1023, finding it was bright in X-rays, which supported Thorstensen and Armstrong but wasn't definitive. Thorstensen kept an eye on the system, taking a spectrum now and then to see if it was doing anything new. It wasn't. Then we came along.

When I found the pulsar, we weren't very sure of its position. We were using the low frequency of 350 MHz so that our beam would be big enough to see any piece of sky for two minutes as the Earth turned. But that big beam means that when we find a pulsar, all we know is a very rough position. Nevertheless, when we found a bright fast millisecond pulsar, we knew it was going to be interesting, so we requested a follow-up observation with the GBT.

The night before we were supposed to take the data, Jim Condon of the NRAO emailed us to point out that J1023 was in our beam. Ingrid Stairs, one of my collaborators on the survey, did some reading and found Thorstensen and Armstrong's paper. Suddenly there was the possibility that our millisecond pulsar might actually have been accreting in 2001. This was a big deal - Duncan Lorimer bet us all a drink that it wasn't the same source (covering his bets, I think - he was hoping it was the same source as much as the rest of us), and we were all excited. So we definitely wanted to find out whether it was the same source as soon as possible.

Unfortunately, the observation we had planned was another 350 MHz observation, just to see whether the source was real and to start to build up a timing solution for it. So as Scott Ransom and I prepared to run the observation, we argued: I wanted to look at 350 MHz first, to make sure there was something there at all, and Scott wanted to take a 2 GHz observation pointed at the position of J1023, so that the much smaller beam would tell us whether the source was really J1023 or not. In the end we compromised: we started at 350 MHz, and the pulsar came booming in right away. So we retracted the prime focus arm and switched receivers and pointed the GBT at J1023, and sure enough, there it was, loud and clear right at the position of J1023. We took the rest of the observation at 2 GHz, and immediately began requesting follow-up observations.

We initially planned to follow up with the GBT, since it was what we were used to, but Paulo Freire emailed us and asked us to please propose for time with Arecibo: there was not much else for Arecibo to look at at the time of day 1023 is visible there, and Arecibo's funding is threatened and it could really use a splashy discovery. (We were keeping the discovery quiet at this point, but Paulo was sharing an office with one of my collaborators, so there was no keeping it from him.) With Arecibo, this already-bright pulsar comes in beautifully, and we get nice clean timing.

Timing the pulsar, we quickly came up with a model of the orbit of the pulsar, and sure enough it agreed with Thorstensen and Armstrong's orbit. In fact, not only did the orbital period come out the same, if we extended our solution back to 2005, we got the same orbital phase as they did: over those three years, we were able to account for every single turn of the companion around the pulsar. Needless to say, this removed all doubt that our system was actually J1023.

The follow-up observations also revealed some peculiar phenomena, like plasma floating around the system and orbital period variations (very small, needless to say), but the essence of it is there: J1023 is a system with a neutron star and a companion that had an accretion disk for about a year around 2001 but is now a millisecond pulsar. The paper has just been published in Science, and is available on arxiv.org.

Unobtainium capacitors

Continuing my train of thought from the other day, what are the theoretical limits for a capacitor? Just how much energy can you store in one?

If you try to analyze it in terms of voltages and capacitances, it becomes very tricky: you can make thinner plates closer together without increasing the weight or volume of the capacitor at all. The key to making the problem tractable is to realize that the energy is stored in the electric field. The goal in trying to store as much energy in a capacitor as possible is to get the electric field as strong as possible over as much volume as you can manage.

Since this will be an order-of-magnitude estimate anyway, let's assume that the electric field completely fills the capacitor. That leaves the maximum strength as the last parameter. Since we're looking for theoretical limits, let's assume you're limited by the tendency of the electric field to ionize the atoms in whatever you make this capacitor out of. So an electric field of about a volt per angstrom is just about all we can hope for. This gives us an energy density of about 400 Joules per cubic centimeter, or 0.1 watt-hours per cubic centimeter. This really isn't that much, and in fact according to Wikipedia there is a company claiming its capacitors do substantially better than this. Granted, it's a patent application, but where did I go wrong?

Well, there are certainly things I ignored - dielectric constant, for example - but I think the basic issue is thati in real capacitors the limit really is the strength of electric field your materials can tolerate. Since the energy density depends on the square of this number, clever engineers who can come up with extraordinarily tough materials have a fair amount of room to beat my crude estimate. But it seems clear to me that that room is far from unlimited: capacitors will continue to improve as energy storage devices, particularly in aspects I haven't touched on like leakage, but don't expect one to power your laser pistol any time soon.

A Missing Link

We've been pretty sure that low-mass X-ray binaries turn into millisecond pulsars for a while now. But no one had ever seen the one turn into the other. Thanks to the drift-scan survey, now we have.

The system, J1023, is currently behaving just like a normal millisecond pulsar: it spins very regularly, and we see radio pulses every cycle. But in 2001, before anyone had any idea there might be a pulsar there, optical observations showed that it was flickering, blue, brighter than usual, and had double-peaked emission lines. All of these are signs of hot gas flowing in the system, and the double-peaked lines in particular are a pretty clear indication that the system had an accretion disk. Unfortunately nobody knew to point an X-ray telescope at J1023 while it was doing this, but it seems pretty clear in retrospect that it was an X-ray binary at the time.

In 2002, observations showed the emission lines were gone, the brightness and color were back to 1999 levels, and the flickering was tailing off to nothing. So the optical story stands today: no sign of an accretion disk. But there is clearly a radio pulsar there, and we're able to make some extremely good measurements of the system because of it. We know that the pulsar is 7.2 times as massive as its companion, for example, and if the pulsar has a typical neutron star mass of 1.4 solar masses, the system is 1.3 kiloparsecs away (about 4200 light-years) and we're seeing it at an angle of 46 degrees.

The image above is a pair of artists' renditions of it (well, sort of, I did most of the work using the software binsim and I'm hardly an artist, though Joeri van Leeuwen improved them significantly). They assume the pulsar's mass is 1.4 solar masses, and show the system as we think it was in 2000 and now. The hot disk of matter, present in 2000 but absent now, produced the optical emission that was observed, but (we think) it was blown away when a drop in the accretion rate allowed the radio pulsar to turn on, producing not just the beam of radio waves pictured in the image but also a powerful wind.

Since we know the system went through one roughly two-year active phase, it seems entirely possible that it will do so again within the next few years. If that happens, we'll be able to watch the formation of an accretion disk in a system where we have very good measurements of the orbit and system geometry from pulsar timing. That's never been seen before, and will be very exciting.


While I was the one who found this source, everyone involved in the drift-scan survey did essential work in making it possible for me to find it, and I worked with many other people to carry out the follow-up observations, so let me thank all my collaborators: Ingrid H. Stairs, Scott M. Ransom, Victoria M. Kaspi, Vladislav I. Kondratiev, Duncan R. Lorimer, Maura A. McLaughlin, Jason Boyles, Jason W. T. Hessels, Ryan Lynch, Joeri van Leeuwen, Mallory S. E. Roberts, Frederick Jenet, David J. Champion, Rachel Rosen, Brad N. Barlow, Bart H. Dunlap, and Ronald A. Remillard.
The paper describing the discovery has just been published in Science, and can also be read on arxiv.org.

Dar Williams

I only really started to listen to music in the 90s, so my idea of what music should be is bleeps and clicks, with human voices only in the form of ironic movie samples. (This is of course precisely the range of music that can be made by geeky shut-ins with piles of computers and electronics. Literally, in some cases; I saw this one live.)

But as the result of a hard disk crash (sadly not recorded), I've been listening to a lot of music from other sources lately, and I find I'm really liking Dar Williams. She's got a nice voice, clever lyrics, and how can you dislike someone who writes a song about the Milgram experiment?

Where do millisecond pulsars come from?

Most pulsars have periods around a second, and are spinning down very gradually. A few are clearly young (still hot, still in a supernova remnant, spinning down rapidly) and are faster, with periods as short as several tens of milliseconds. But there are other pulsars that are even faster - periods down to 1.4 milliseconds, that is, spinning some seven hundred times a second - and yet appear to be very old: cold, spinning down very gradually, very weak magnetic field. These are called millisecond pulsars, and it was a puzzle where they could possibly come from.

The key clue to the puzzle was their binarity. Most stars are found in systems of several stars orbiting each other. But pulsars are usually solitary. This is mostly because of their violent births: to make a pulsar, the star has to go supernova, and such a violent explosion tends to break up binary systems. It doesn't always break them up, though, and so some binary pulsars are known. If you look at the millisecond pulsars, though, most of them are actually in binary systems, unlike the normal pulsars (see the graph at right, based on the ATNF pulsar database). So we think the presence of a companion is key to making a millisecond pulsar.

The story, as we understand it, is this: A pulsar forms in a supernova. It is either in a binary system which survives, or it captures a companion. It lives out its life as a pulsar, spinning down gradually past the point where it is visible as a pulsar. The system stays like this for a long time. But eventually, one of two things happens: the companion starts to swell up into a red giant, or the orbit shrinks. In either case the system reaches a point where some of the matter at the surface of the companion is attracted more strongly by the neutron star than by the companion. This matter then falls down onto the pulsar.

Pulsars are not much more massive than the Sun, but they are much much smaller. So when matter falls onto one, tremendous amounts of energy are released. But remember that the two stars are in orbit, so the system is rotating. If you take a piece of matter from the companion, it will carry some angular momentum with it. To make the obligatory analogy, just as when a figure skater pulls in her arms, she speeds up, when you take matter from the companion and bring it in to the neutron star, the matter begins to rotate more rapidly. When this matter falls on the star, the star is spun up a little. I am glossing over numerous important details here, but the point is, when you start transferring mass to a neutron star, it is possible to spin it up.

So, we think that millisecond pulsars are old pulsars that have spun up, or "recycled", by accretion of mass from a companion. Observationally, we see systems where mass is being transferred: they're very bright X-ray sources. In a few cases we can actually tell how fast the neutron star is rotating, and sure enough, its period is down in the millisecond range. But these systems don't produce radio pulsations, presumably because all that matter falling in either blocks the radio waves or shorts out the emission mechanism (which needs a near-vacuum in the magnetosphere). So to make millisecond pulsars you need somehow to turn off the accretion and clear out all the matter, so that the radio pulsations can emerge. This transition hasn't been seen before, and the theorists have some difficulty explaining the population of objects that we see - while we see both millisecond pulsars and their hypothetical accreting progenitors, none of the progenitors seems to be positioned to turn into anything like the millisecond pulsars we actually see. So that last step, accretion turning off and radio pulsations starting, remains something of a mystery.

Energy storage

Dan of Dan's Data just posted an article about science-fiction-y batteries. He discusses a battery based on matter and antimatter but dismisses it as inconvenient (in that the energy comes out as a shower of high-energy radiation). Instead he suggests a battery based on just cramming a lot of electrons - several grams of electrons! - into a AAA-sized package, and finds that this actually stores vastly more energy. There's an error there, but it's an interesting one.


First of all let's do his calculation a little more carefully. Let's suppose we're taking 11.5 grams of electrons - roughly two billion coulombs - and cramming them into a AAA-sized package. For simplicity, let's assume the package is a conductive sphere of radius 1 cm. This is actually a capacitor, with capacitance about 1 picofarad (incidentally, 1 pf is a fairly normal size for a capacitor, though they don't usually choose a conductive-sphere form factor). So the energy stored in this capacitor is (1/2) Q^2/C, or a staggering 2*10^30 Joules, the sun's total luminosity for a couple of hours.

(All calculations are done with the charming frink, though of course the errors are my own.)

There's a problem with this, though. Actually, there are several, but let's tackle the most fundamental first: if I convert 2*10^30 J to mass using E=mc^2, I get about twenty billion tonnes. This isn't just a meaningless figure: the energy stored in the capacitor really does weigh this much. If it seems weird for energy to have mass, well, I have to agree, but you can actually see it by looking at the periodic table: a helium atom is basically four hydrogen atoms after you stick two of the electrons to two of the protons and squash all the nucleons together. But each hydrogen atom has a mass of 1.007825 u, while helium has a mass of 4.002602 u, rather than the 4.0313 u that you would expect. The difference in mass is the mass of the energy in the nucleus. Since in helium, that energy is negative (pulling the nucleons together rather than pushing them apart), the mass is also negative, and helium weighs less than the hydrogen that makes it. But it's the same phenomenon.

One interesting question is, where is all that mass? It turns out that there is a nice answer. The little conductive sphere is surrounded by a strong electric field. In fact Gauss' law lets us work out just how strong it is at the surface: 2*10^23 V/m. It turns out that electric fields store energy. The stronger the electric field, the higher the energy density it stores. If you add up all the energy stored in the electric field of this charged sphere, over all space, you get exactly the energy you had to put in to charge it up. In other words, the energy is stored in the electric field, strongest at the sphere but extending out to infinity. And when you have an energy density, you have a mass density. Dividing by c^2, you get a mass density for the electric field of 2*10^15 g/cm^3. All that mass is due to the electric field, mostly close to the sphere.

Incidentally, there's a neat little calculation here. We know that electrons are weird quantum-mechanical beasts whose position and velocity are unavoidably ill-determined, and in fact as ar as we can tell an electron is a point particle. But let's ignore that for a moment and try to figure out how big an electron is. We don't have much to go on, because they don't sit still, and they interact with other things without touching them. But there is one place to look: they have a mass. Suppose we pretend that an electron is a tiny conductive sphere, with no mass of its own, but with charge smoothly distributed over its surface. Then the electric field will have a mass of its own, just as with our little conductive sphere. What if we just declare that the mass of the electron is entirely due to this electric field? Well, then we can calculate a radius! The value you get is 1.4*10^-15 m, which is tiny enough (about a hundred thousand times smaller than an atom) to not be totally unreasonable. You get a slightly different answer if you assume the electron is non-conductive, but leaving aside that issue, this size is the classical electron radius: 2.8*10^-15 m. Astonishingly enough, this number comes up in various physical phenomena, for example Thomson scattering of light by an electron.

So we can't, even with some pretty impressive unobtainium, cram 11.5 grams of electrons into a AAA-sized space. But if we decide that 11.5 grams should be the total weight of the battery, including energy, how much energy can we store? Well, it turns out, just a little less than the matter-antimatter battery. Since the matter-antimatter battery turns its contents entirely into energy, there's no way to beat its storage efficiency. But the above calculation shows that the charged metal sphere is nearly as efficient. We just charge it up to a slightly less outrageous voltage of 4*10^13 V. Now we've got something as good as the antimatter battery, but that provides us with handy electrons! What more could you ask?

Well, for one thing, you need to get rid of the electrons somehow once you've used them. Otherwise your laptop will be constantly giving out static shocks. Try and drive an electric car with them and you'll be shedding kilowatts of static shocks. Not ideal. Also, speaking of static, the electric field of this metal sphere will extend out to infinity, and will be quite strong even at some distance. At the surface of the sphere, the electric field will be 4*10^15 V/m, which is 4*10^5 volts per angstrom. So the voltage difference across an atom will be on the order of a million volts, compared to the 13.6 volts it would take to pull an electron off a hydrogen atom. So this little sphere will shred all atoms that come near it. Worse, that field extends out to infinity, so that at a meter away it's still 40 volts per angstrom. This is going to be a real pain. Isn't there some way to shield it?

Well, yes. You just put a spherical shell around it with the opposite charge. Or you abandon the sphere and just set up two parallel plates with opposite charges. Now the fields of the two plates cancel out (nearly), everywhere except between the plates.

Actually, this kind of energy storage device exists, and you probably have thousands in your house: it's a capacitor. Now, granted, it doesn't store the kind of energy density Dan is talking about, since it's not made of unobtainium by aliens, but an appropriate bank of capacitors can store a very great deal of energy and deliver it very quickly.

But it doesn't sound so impressive to ask the aliens for a capacitor.

[Update: I talk more about what one can hope for from a capacitor here.]

Finding pulsars part 2

To find new pulsars, you obviously need to look at the sky with a big radio telescope. But once you've recorded this massive amount of data, you're only half done. Searching through the results, it turns out, is a massive computational task, even with modern computer clusters.

In principle, searching for pulsars in the data should be simple: they're nice periodic signals, so you just take a fast Fourier transform and the periodic signals show up as peaks in the power spectrum. In practice it's a lot more complicated than that.

The first problem you run into is interstellar dispersion. The space between stars is not empty: there is something like one hydrogen atom per cubic centimeter, on average. It is mostly hot, and at such a low density that means the hydrogen is in the form of plasma, loose electrons and protons. Plasma has a refractive index: it slows radio waves propagating through it very slightly. The absolute amount is small, maybe a few milliseconds, but it depends on the observing frequency. So if your pulsar has a one-millisecond pulse, then you might find that at the bottom of your observing band it's delayed by several milliseconds compared to the top of your observing band. This smears out a nice sharp pulse into a smeary blob that's much harder to detect. So when we look for pulsars, we need to divide the observing band into many channels and record the power in each channel separately. When you multiply a thousand channels by ten thousand power measurements a second, you get an awful lot of data.

To search for pulsars in this huge stream of data, you have to choose a dispersion measure (DM) and shift the channels accordingly, then add. Now you can do your Fourier transform and look for peaks.

The second problem that arises is acceleration. Some of the most interesting pulsars are not solitary stars. Instead, they are in orbit with a companion. Usually the pulsar is the more massive of the two, but it does not stand still; both stars orbit their common center of mass. Since these binary systems are often very close together and often have quite short periods - hours to days - the pulsar can be moving quite quickly, which produces a Doppler shift. A constant Doppler shift just changes the position of the peak, which is not a problem, but if the pulsar is accelerating - and moving in a circle involves acceleration - the Doppler shift can change. If the system is close enough, the Doppler shift can change substantially over the course of an observation, so that instead of the pulsar appearing at a single frequency in the Fourier transform, it gets spread out. Now, this may sound rare, and it is, but not only are these some of the most interesting systems to find, they are also the systems most likely to have been missed in previous surveys. So it's important to look for these accelerated systems, which you can do by looking not just for single peaks but for the characteristic shape of an accelerated peak. This is effective but takes substantial computing time.

Multiply together all the DM trials, Fourier frequencies, and possible accelerations, and you find that processing a single two-minute beam takes some thirty hours on a modern computer. Work out how many two-minute beams there are in the 1491 observing hours we carried out for the drift-scan survey, and you begin to see why we haven't yet processed all the data.

Once you've processed the data, though, you unfortunately do not get a nice tidy list of pulsars. There are rather a lot of processes that produce periodic radio signals, so what you get is a monster list of cell phone calls, car ignitions, telescope birdies, electric blankets, and a very few new pulsars. Winnowing through the huge mass of periodic signals is unfortunately less a question of clever Fourier techniques and more akin to spam filtering: many are obviously junk (they don't show interstellar dispersion, say), some are less-obvious junk, and a few are hard to tell from real pulsars even by humans.

My main contribution to the survey has been to run data through McGill's computer cluster and to sift through the results. I've done a certain amount of sifting by hand, but I've also written a few "spam filters" to try to reduce the number of obviously-not-pulsar results I have to look at. The real excitement, though, has been finding new pulsars and following them up.

Finding pulsars

How do we find new pulsars? Well, there are a few places we know pulsars are likely to be - globular clusters, supernova remnants, X-ray binaries - but for the most part we just need to look at lots of sky. This isn't really feasible in X-rays - which you can only observe from satellites, on which time is precious - and apart from the exciting new gamma-ray results from Fermi, the only frequency range at which you can reasonably hope to find pulsars is the radio. So we more or less need to run a big radio survey that looks at lots of sky, one little patch at a time. But how do you get time on a multimillion-dollar telescope to look at sky where there's probably nothing?

Sometimes you can make a good science case for doing a survey. The most successful pulsar survey, the Parkes Multibeam survey, did this: they had a new multibeam receiver that would drastically improve the rate at which they could survey sky, so it was clear that time spent looking for pulsars would pay off. And it did, handsomely; they discovered over a thousand new pulsars.

More often, though, there are more pressing projects fighting for telescope time. So pulsar surveys need to come up with some other way to get time. The GBT 350 MHz drift scan survey took a clever approach: we were able to use time the telescope wasn't able to be used for anything else.

The GBT is an amazing machine, but with seven million kilograms running back and forth constantly, the azimuth track recently started to wear out. So the summer of 2007 was spent replacing the track. During this time, obviously the telescope couldn't move in azimuth without risking falling off the end of the track, so it couldn't follow sources as the Earth turned. But the receivers were all working fine, and the telescope could move in elevation just fine. So the people behind the drift scan survey pointed out that if they could simply use the receivers during the repairs, recording the sky as it rolled past, they could get useful data during the summer. This made good sense to the telescope access committee, so they were given the time.

The final arrangement was that the workers would work four twelve-hour days a week, during which time no data would be taken (arc welders are not conducive to good radio reception anyway). The rest of the time was devoted to the survey, which would be allowed to move the telescope in elevation twice a day. The survey used a relatively low frequency, 350 MHz, so that the beam would be big enough for any given piece of sky to spend about two minutes in view of the telescope.

Once the survey was set up, all we needed to do (I say we because it was at about this point that I got involved) was move the dish once a day, start the data recorders, and deal with the data. Dealing with the data was more involved than it sounds, because the data came in at a rate of about 90 GB/hour. So leaving analysis aside for the moment, we filled a couple of 800 GB disks a day, and somebody had to move all the data onto them and box them up for shipping to the various universities that would analyze the data. So when I went down to the telescope to observe as part of the survey, what I actually did was run a couple of scripts a day and shuffle disks.

In the end, the survey covered roughly a third of the celestial sphere, producing 134 TB of data. We're still analyzing the data - maybe only a quarter of it has actually been processed at this point, a year and a half later - but we've already found a number of new pulsars, some of them fascinating. So I think there are worse things the GBT could have done with those two months of track repair.

Safety

I certainly understand the value of workplace safety. But really, two pages of alarming safety warnings for a java tone generator applet?

Also, among the very numerous warnings:

All sounds produced by DFG modules are licensed by Digital Recordings for your personal, environmentally-responsible and non-commercial use only.

"All sounds" is just sine waves, so I'm not exactly sure they can claim copyright on them, let alone use copyright restrictions to prevent me from damaging the environment (with sine waves?!)

Falling free

There was a rather neat paper on arxiv.org today: The Flow Of Granular Matter Under Reduced-Gravity Conditions, by Hofmeister et al. They point out that while we understand the movement of, say, sand fairly well under terrestrial gravity, objects like the Moon and Mars are covered with fine powdery regolith, and we don't really understand how it moves under the lower gravity there. This includes, for example, the steepest stable slope, obviously important for understanding the dunes on Mars. So they did some experiments in reduced gravity - which is much harder than it sounds.

Gravity is tricky stuff. It's such a weak force, you can't manpiulate it by moving sources around, like you can with electricity and magnetism. You basically need to use whole planets to make any appreciable amount of gravity. But you can fake it with accelerations - and in fact it turns out that "fake" isn't the right word, since accelerations are indistinguishable from gravity. So if you want to do some experiment in higher than normal gravity, you can put it in a centrifuge. But if you want less than normal gravity, that's not so easy.

There are basically four ways to get reduced gravity. The first is to fake it: many experiments (wind tunnel tests for example) are done on scale models. The scaling requires cleverness, because (for example) the surface-to-volume ratio of a toy airplane is very different from that of a real airplane. So you also need to scale things like air speed and air density appropriately. So if you're studying some gravitational phenomenon, fluid flow perhaps, you may be able to build a model scaled so that the needed gravity for the model is one Earth gravity. This only works for systems with pretty simple physics, though, for which you understand the scaling.

The second way to get reduced gravity is to do your experiment on the ISS. Unfortunately, it costs hundreds of thousands of dollars per kilogram to get anything up there, so this is rather a last resort.

There is also the famous "Vomit Comet", which flies along parabolic trajectories to provide about half a minute of microgravity to those aboard. This is great for training astronauts (and wealthy space nerds) but it's very difficult to keep residual accelerations small enough for precision experiments.

This brings us to the system Hofmeister et al. used: drop towers. In principle they're very simple: you put your experiment in a box, then drop it from a tall tower. For as long as it takes to fall, you have microgravity. If you want to double the time, you can even fling it up from the bottom. Of course, the devil is in the details; for example you want to evacuate the tower so there aren't residual accelerations from air drag, and you need to plan on the capsule being pushed a few centimeters to the side by the Coriolis force. Plus, of course, you need to stop it non-destructively when it reaches the bottom. But all these problems are tractable, and in Bremen there is a drop tower that can give up to nine seconds of free-fall.

Hofmeister et al. used the Bremen drop tower for their experiment. Perhaps perversely, in order to get the reduced but non-zero gravity they needed, they put their experiments in a centrifuge inside the capsule. They hooked up a high-speed camera, and were able to track each grain of sand as it slid down the slopes. They found that the steepest stable angles and velocity patterns of the flows were rather poorly described by existing theory. So it looks like the theorists have some work ahead of them before they can understand the Martian dune fields.

Telescope of the month: the GBT

The Robert C. Byrd telescope at Green Bank in West Virginia, called by almost everyone the GBT, is a radio telescope. With an elliptical dish 100 meters by 110 meters, it's the biggest fully-steerable telescope in the world - much to the chagrin of the people in the tourist center at Effelsberg, whose telescope is 100 meters but circular.


As a single-dish radio telescope, the GBT really only puts one pixel on the sky at a time, so making an image is a slow and laborious process. People do it, but where the GBT really shines is for looking at other aspects than spatial resolution. It really shines when looking at spectral lines - which lets you map out moving gas along the line of sight - and for pulsar observations.

The GBT is also used in conjunction with Arecibo for radar experiments: Arecibo has a big radar transmitter that can be used to bounce radio waves off objects in the Solar System. But since transmitting and receiving at the same station can be a problem, often projects will use the GBT as a receiver for these radar transmissions.

As a machine, the GBT is quite astonsihing: seven million kilograms of moving mass. The people at the visitor center like to say that it's the largest moving object on land, which is not quite true, but the counterexamples - the Saturn V aboard its transporter, and a colossal dragline excavator so big it walks instead of rolls - put the scale in context.

One of the most valuable features of the GBT is its location in the National Radio Quiet Zone. Not only is there no cell phone or pager service, the NRAO has a converted ice cream truck they drive around searching for sources of radio noise (for example in one story they tracked down one noise source to a local family's malfunctioning electric blanket; the family didn't want to part with it so the observatory bought them a new one). This absence of local interference makes a tremendous difference when looking for new pulsars (among many other tasks), and sure enough surveys at the GBT have been very productive of new pulsars.