Possible solution to a major problem.

Allowing me to do this (left: original; right: modified reconstruction).

Given image i, such that i = a + n + c, where a is the astronomical signal, n is the random noise, and c is the detector specific corruption, define a new image i' = (a - a') + n + c, where we've removed an estimate of the astronomical signal a'.  This has a Fourier transform F(i') = (A - A') + N + C, using the fact that Fourier transforms are linear.  Assuming a' is a sufficiently suitable proxy for a, and that n is small (specifically that ||N|| << ||F(i')||), this leaves F(i') ~ C.  Constructing a Fourier mask that removes C results in t ~ (a - a') + n + c.  The difference of i' and t provides a clean estimate of c, which can then be removed from i: o = i - (i' - t).

I need to test that photometric properties are retained by this transformation, and I'd like to speed it up some and sort out some of the edge case issues (visible around the blank bands above), but I think this solves the bulk of the issues with this problem.

The Supreme Court

Specifically, killing them.

Hang on there, Secret Service.  We're doing it statistically, so it's just an math problem.  No one is actually getting hurt, and I mostly did it to see how soon it's likely that there will be a vacancy due to death.  Someone could resign at any time, and that would make a vacancy, but that's not something that can really be predicted.  Death is a bit tricky, too, but due to the fact that despite the best attempts of some people, we still have a government that collects piles of data, so we can use the cohort life expectancy data from the Social Security Administration.

From that data, we're going to use the "High-cost" columns, and use the "at age 65" data.  The logic behind this is based on two facts.  First, the majority of the difference between the "at birth" and "at age 65" expectancies is due to infant mortality and things like that.  If you survive to 65, you didn't die as a baby.  The supreme court is universally comprised of people who did not die as a baby.  The choice of the "High-cost" is that being a supreme court justice is a pretty good job, so they're likely to have good health care.  The final assumption I have to make is that people born before 1940 have life expectancies similar to those born in 1940.  The table only goes back so far, and I don't really want to extrapolate it.

The next assumption is that the probability of surviving a year is based on the logistic function, largely due to a lack of better options.  The form I chose was P(age) = 1 / (1 + exp(-(expectancy - age))).  This sets the 50% probability at the expectancy age.  I could scale it, but I don't have any good idea of how to do that, so I'm just letting that be 1.0.

Using this table of birth dates and ages, we can make the following plot:

This shows the probability that that justice is alive.  I'm sorry that this claims Ruth Bader Ginsburg is 52% dead.

Ok, using these probabilities, we can plot up the number of justices statistically alive by simply summing these curves:

This suggests that (barring the Senate being jerks), Obama is likely to appoint two justices.  His successor falls into a generational gap between the Breyer and Thomas, and isn't likely to appoint anyone, even in an eight year term.  It's fairly linear over the next five terms, with about one each term.  This assumes that all newly appointed justices are younger than the current youngest, and that we're only concerned with the current set.

So again, barring the Senate terrible, we should have a significantly more liberal court within the next four years.

Batman and Robin are happy about that fact, too.

Because fundamentally misunderstanding physics annoys me

I just read an article where some moron is writing that people should build their own anti-satellite system to blind spy satellites using laser pointers.

Thinkgeek sells this laser pointer with a beam divergence of 1.2 milli-radians.  

Ok, fine, this page lists a few others, and it looks like the best they can come up with is a 0.5 mrad.  Geosynchronous orbit is 35786 kilometers up.  Using the power of tan(x) ~ x for small x (and noting that 5e-4 radians is really small), this gives a beam diameter (assuming a point source at the laser) of 17.893 kilometers.  This seems to work out to the irradiance at that satellite is something like a factor of 3.1234e-15 smaller than it is at the source.   In conclusion,

Thank you, Lois.

More Fourier transform stuff.

Consider this image:

This is just a sine wave tipped at 45 degrees.  Now, look at the power spectrum:

Ok, the two spots make sense, as those correspond to the point in the (U,V) plane that has that frequency, along with it's reflection (the negative frequency).  But what's up with all that other crap?

FFTs make the assumption that the signal is periodic, and since it only has the image to work with, the default thing to do is to calculate a circular FFT, where the left and right edges wrap around, and the top and bottom do as well (making a video game torus).  However, since the sine signal doesn't nicely match up at the edges (most obvious at the top and bottom), this introduces a lot of extra Fourier components that correspond to those + shapes.

The solution to this problem is to taper the image down so that it is effectively periodic.  The easiest way to do this is to just multiply by a Gaussian so that the edges drop down to zero at the edges, like this:

Ok, what's this look like in (U,V)-space?

The two dots we expect.  There's some lower level noise, but this has the same colormap as the previous spectrum.  The dots aren't perfect points, which I think tells us something about the size of the Gaussian we used.

This is important, because the science stuff I've been playing with has odd edge effects that impose a strong sinc^2 signal onto the power spectrum.  This makes it difficult to isolate bad frequencies.  Remembering this trick solves that problem, and should make cleaning those images much easier.

Another reason why Fourier transforms are magical.

Let's pretend we have an image, and that image is just a bunch of dots:

Here are three, with different numbers and positions for the dots (randomly distributed), except for the centered single dot case.

I've chosen Gaussians for the shape of the dots, because those are simple to work with.  Now, take the 2D Fourier transform of these three images.  What does the power spectrum look like?

Those curly-cue shapes are ~20 orders of magnitude smaller than the peak.  That's what I get for plotting logarithms.

To within the floating point accuracy of my computer, these power spectra are all identical.  No matter where the dots are on the image, the power spectrum is identical.  How does that work?  Basically, the power spectrum just says "the image is this big, and the things that show up on it look like this."  All of the position information is stored in the phase information between the real and imaginary components:

The single dot has a constant phase image.  The other two have phase images that look like random noise, but are actually encoding the position information for each dot.  Next interesting thing:  what happens when we increase the size of the dots:

Bigger dots have a smaller sized dot in the power spectrum, as wider things are dominated by smaller frequency terms (also the fact you can directly calculate that F(Gaussian(sigma)) = Gaussian(1/sigma)).

Measures of variability

Actually: estimates of the dispersion of a distribution.  Everyone knows the standard deviation

sigma = sqrt(1/N * sum(x_i - mean))

That's pretty much the one everyone uses.  

Another strategy is to use the interquartile distance:

IQD = x_75% - x_25%

If you look at this for a Gaussian distribution (which is the assumption the standard deviation makes), you can use the cumulative distribution function for the normal distribution to see that you can get a "sigma" from this where

sigma_IQD ~= 0.7413 * IQD

This measure is somewhat more robust against outliers, as you can pile a bunch of points outside the inner half of the distribution, and they won't change the IQD much.  However, if you have a bias where all your outliers are on one side (say, all greater than the median of the distribution you care about), then you can only have 25% of all your points being outliers.

One way to fix this is to use the median absolute deviation:

MAD = median( abs(x_i - median(x_i)))

It's not too difficult to directly prove that for a symmetric distribution (Gaussian, or a Gaussian + unbiased outlier distribution)

MAD_symmetric = 0.5 * IQD_symmetric

However, if the outliers are biased, MAD is able to accept up to 50% ourliers before it really breaks down.  As above, you can create an effective sigma:

sigma_MAD ~= 1.4826 * MAD

Finally, you can histogram all the data, and do a Gaussian fit to that histogram to determine the width.  This has the problem that you need to have a decent number of samples before the fit will converge well.  It also can be more computationally expensive, since you have to fit a function.  Ok, now let's look at a bunch of plots.  For these, I created 10000 random samples, with some fraction "outliers."  The real distribution is a Gaussian of mean 0 and stddev 1.0.  The outliers are drawn from a flat distribution from -5 to 5 (for the unbiased case) or from 0 to 10 (for the biased case).

Biased.

Unbiased.

From these, you can see that:

1. For the biased outlier distribution, IQD fails around 25% contamination.  MAD is better until around 50%.  Histogram fits are generally better, but fail catastrophically around 60%.
2. For the unbiased outliers, IQD and MAD are identical as suggested above.  Histograms are still better.
3. Standard deviation is universally bad when any outliers are present.

Biased.

Unbiased.

This is the same set of plots, now done for N = 100.  Histograms are now good only in the unbiased case with less than 50% contamination.  For the biased case, they're not better than the standard deviation above 20% contamination.

So that's largely why I like using a sigma based on the MAD.  It's pretty much universally better than the regular standard deviation when outliers are present.  Fitting Gaussians to histograms can do a better job, but that uses more complicated math, so it may not be useful if you don't have fitting code ready to go.

Recovery

Let's assume a disaster strikes.  You have an incredibly efficient logistics team, who have access to large ships, air freight, and helicopters capable of dropping supplies into any location needed.

Assume you pack everything in a 20' shipping container.  This gives you 33.1 m^3 of volume, and a net load of 28,200kg.  What do we need to pack into this?

1. Water.  1 m^3 of water weighs 1000 kg, and contains 1000 liters.  Woo metric system.  ($6.20 / 12L)
2. Food. Let's go with MREs, and this page, so a case is 12 MRE/case, which weighs 20lbs and is 16"x11.5"x10".  Metricize, and we find that 12MRE = 9.07185 kg and 0.0301522 m^3. ($60/case)
3. Blankets. Let's get some of those space blankets.  A case of these contains 200, 0.2lb each, case size 17.5"x9"x10". ($1.35ea)  Let's do regular blankets too.  6"x6"x15", 3.6lb, $17.
4. Other supplies.  Let's say we want to reserve 25% of the container for medical supplies/generators/communication equipment/etc.

Ok, now solve.  One person needs about 2 liters of water a day under regular circumstances.  Let's bump this to 10L/day because you're rebuilding civilizationStaten Island.  One person also requires 3 meals a day, and should probably have one space blanket and two regular blankets.  Let's bump that to four regular blankets, and assume the excess represents extra clothing, etc.  This gives us the equations:

Usage = [W, M, S, B] = [10 * D * N, 3 * D * N, 1 * N, 4 * N] 

Weight (kg)  /unit   = [1        , 0.75599  , 0.090718, 1.6329]

Volume (m^3) /unit   = [0.001    , 0.0025127, 1.2905e-4,0.008849]

W = k1 * W_w + k2 * W_m + k3 * W_s + k4 * W_b

V = k1 * V_w + k2 * V_m + k3 * V_s + k4 * V_b

N_people = k4 / 4 = k3 => 4 * k3 = k4

         = k2 / (3D) = k1 / (10D) => k2 = (3/10) k1

W = D * N * (10 * W_w + 3 * W_m) + N * (W_s + 4 * W_b)

V = [...]

[12.267970 6.622318;0.017538 0.035525] * [N * D ; N]

Ok, I come up with plausible solutions:

N  D  Weight     Volume   W / Wc W / Wm   V / Vc V / Vm

58 30 21382.3622 32.5766 [1.0110 0.7582] [1.3122 0.9842]

60 29 21383.6069 32.6476 [1.0110 0.7583] [1.3151 0.9863]

62 28 21335.7796 32.6485 [1.0088 0.7566] [1.3151 0.9864]

64 27 21238.8805 32.5793 [1.0042 0.7532] [1.3124 0.9843]

65 27 21570.7380 33.0883 [1.0199 0.7649] [1.3329 0.9996]

67 26 21412.4990 32.9314 [1.0124 0.7593] [1.3265 0.9949]

69 25 21205.1882 32.7043 [1.0026 0.7520] [1.3174 0.9880]

72 24 21243.8591 32.8635 [1.0044 0.7533] [1.3238 0.9929]

75 23 21208.9221 32.9174 [1.0028 0.7521] [1.3260 0.9945]

Where Wc and Vc are the 75% critical limits, and Wm and Vm are the absolute container limits.

How much is this going to cost?  Let's take the maximal people solution.  This means we need:

75 * 23 * 10 L of water      @ $6.20 / 12 ea = $8912.5

75 * 23 * 3  MRE             @ $60 / 12 ea   = $25875

75           space blanket   @ $1.35 / 1 ea  = $101.25

75 * 4       regular blanket @ $17   / 1 ea  = $5100

1            shipping container @ $3000 / 1ea= $3000

============================================== $42988.75

So, for something in the $50000 range, you can support 75 people for three weeks.  Ok, so we can't solve all of Staten Island's problems with these, as supporting its population of 470467, you'll need 6300 containers, at a cost of about a third of a billion dollars.  Still, this seems like the kind of thing the government could just stockpile in a few places, put them on a train towards where it looks like disaster might be happening, and drop them into place as needed.  Having a dedicated recovery unit that is ready to go anywhere at anytime, the initial recovery impact is likely to be more dramatic, providing a more efficient (and PR popular) recovery.