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.

Distribution widths.

Copying comments from a source file is totally the same as writing a blog post, right? Has anyone really considered the full-width N-max of a gaussian before, and then extended that parameter to non gaussian data?  Here's how they would start, if they did:

/* Discussion on FWHM for Gaussian lines:
   Given a gaussian distribution defined by PDF = 1 / sqrt(2 * pi) * exp(-0.5 * u^2),
   then you can define a fractional max location by:
   PDF / PDF_MAX = exp(-0.5 * u^2) / exp(-0.5 * 0^2);
   PDF / PDF_MAX = exp(-0.5 * u^2)
   sqrt(-2.0 * log(R)) = u
   So, for a given ratio of the maximum R, the half-width R-max is located at u.
   This leads to the full-width R-max location being at
   FWRM = 2.0 * sqrt(-2.0 * log(R))

   Since we've used u, we can add the sigma dependence back in by noting that u = x / s:
   FWRM = 2.0 * sqrt(-2.0 * log(R)) * sigma
 
   The standard value is the FWHM, or R = 0.5:
   FWHM = 2.0 * sqrt(-2.0 * log(0.5)) = 2.0 * sqrt(2.0 * log(2.0))

   Do this for quarter max, R = 0.25:
   FWQM = 2.0 * sqrt(-2.0 * log(0.25)) = 2.0 * sqrt(2.0 * log(4))
        = 2.0 * sqrt(2.0 * 2.0 * log(2.0)) = sqrt(2) * FWHM
   And for 3/4 max, R = 0.75:
   FW3QM = 2.0 * sqrt(-2.0 * log(0.75)) = 2.0 * sqrt(2.0 * log(4/3))
         = 2.0 * sqrt(2.0 * (log(4) - log(3)))
         = 2.0 * sqrt(2.0 * (2.0 * log(2.0) - log(3)))
         =~ 0.64423 * FWHM

   If a distribution is "wingier" than a Gaussian, then FWQM will be inflated, such that
   FWQM_wingy > FWQM_Gaussian

*/

I seem to always have trouble on my return flights from Austin.

I seem to always have trouble on my return flights from Austin.  I've been stuck in Denver and Phoenix, and I was reasonably sure that I was going to end up stuck in San Francisco as well.

I arrived at AUS with plenty of time to spare, and went through the "Premiere Line" for security, which means I get to skip over all the waiting and go directly to the driver's license check.  Yay, frequent flier benefits!

I then got pizza:

which put me behind the worst fast food orderer I've seen in a long time.  First he didn't know what pizza he wanted. Then, after choosing and moving on to pay, he decided he wanted a second slice, and asked if it could just be put into the same box.  Since a box holds one slice, that couldn't happen, so he accepted the second box.

Anyway, by the time I got to the gate, I noticed an important fact: there was no plane connected to the jetbridge.  This continued for almost an hour, until the flight came in late, leading to the second problem: there was no crew to fly the now-available plane to SFO.  I've never been able to understand how airlines can manage their logistics so poorly.  You know you have a flight at airport X at time T.  That flight requires a set of N crew members, who need to be replaced after some number of hours.  Therefore, you need some number of crew sets k = ceil(dT / H).  Not having a crew I think has been the problem in both my Denver and Phoenix flight issues.

Upon arriving at SFO, I discovered the actual reason for the delay for the flight: the airport only had one open runway, and all traffic had to use that.  I didn't know anything about this, suggesting that United didn't do a very good job of notifying travelers.  The one benefit of this was that all flights were delayed, so leaving Austin an hour late didn't decrease my layover very much.  I went to the gate, and saw what I thought was my name noting that I'd been given an upgrade to first class, as seat 3-A was listed for someone with my same first initial and first three characters of my last name.

So, using this set of data for last name frequencies, and this data for first names (assuming first initial fractions are reasonably constant over time, so I don't need to convolve an age cohort function against time variable first names), suggests that my first name is more likely than random (0.068126 compared to 1/26 = 0.038462), but that my truncated last name isn't that common (0.0024997 compared to SMI = 0.010087).  Therefore, randomly selecting from the US population gives me a upgrade name that has a probability of 1.7030e-04.  Now, confining this to a plane makes it a birthday problem, and this seat map suggests a United 757-200 has 188 seats.  This leads to a probability that someone else on the plane has the same upgrade name as me of 0.031512 (Edit: For comparison, Roger Smith has a probability of 0.10684, significantly worse than me.) This is higher than I would have thought, but since it's a birthday problem, that's pretty much always the case.  Still, this would give me a ~97% confidence that I had the upgrade to first class, meaning I would be able to count on a meal on the flight.

I had some emergency chocolate that I'd packed before leaving home, and I discovered it while sitting there.  It had a raspberry filling.

Tl;dr: I didn't get the upgrade, this flight fell into the 3.15%.

In any case, I made it home without much trouble (despite all that), and was able to get to see a really nice sunset over the Pacific Ocean:

I'm sure the clouds have a specific name too. Altocumulus?  Is that right?

Seriously, it's things like this that make it hard for me to sleep