Hello all,
There is a lot happening in my life right now and I am spread quite thin among the various projects that I take an interest in. In particular, I am thrilled to publicly announce on this list that Continuum Analytics has received DARPA funding (to the tune of at least $3 million) for Blaze, Numba, and Bokeh which we are writing to take NumPy, SciPy, and visualization into the domain of very large data sets. This is part of the XDATA program, and I will be taking an active role in it. You can read more about Blaze here: http://blaze.pydata.org. You can read more about XDATA here: http://www.darpa.mil/Our_Work/I2O/Programs/XDATA.aspx
I personally think Blaze is the future of array-oriented computing in Python...
Passing the torch of NumPy and moving on to Blaze
 |
| Rough Hex Output & Smooth Surface Reconstruction for Dogleg |
There was a pretty cool amateur rocket shot recently that was an attempt to win the Carmack micro-prize. The rocket is called Qu8k, designed and built by Derek Deville and friends. One of the stipulations of the prize is collection of a GPS location by the onboard avionics at an altitude above 100kft (as long as the velocity is low, this should theoretically not require an unrestricted GPS).
Of course, since the other stipulation of the prize is a detailed report about the shot and the data collected, this gives us number crunching nerds a neat data set to play with. Derek posted a spreadsheet of the accelerometer data, along with a simple first order integration (twice) to get velocity and altitude. I tried an FFT-based method to compare against the first order approach in the spread-sheet, and a second-order trapezoidal rule integration. The python script to do the integration and make the two plots below is on github.
The errors in the numerical integration are not terrible (the plots of the trajectories using the different approaches are indistinguishable in the eye-ball norm).
One of the cool things about Python is the array masking capability. That makes implementing the temperature ratio component of the 1976 US Standard atmosphere (to estimate Mach number) 7 lines of code.
This is an extension to the FFT-based oil extraction model (see the Mobjectivist blog for more more details). The basic approach remains the same, each phase of the process is modeled by a convolution of the rate in the previous phase with an exponential decay. Now we are going to apply some of the ideas I presented in Uncertainty Quantification with Stochastic Collocation.
Here is the original Python function which applied the convolutions using an FFT,
and here is the function modified to use the complex step derivative calculation method:
Notice that we need to transform the real and imaginary parts of our solution separately so that the nice subtractive-cancellation-avoidance properties of the method are retained (see [1] for a more detailed discussion).
Now we have the result of the convolution in the real part of the return value and the sensitivity to changes in rate in the imaginary part (see Figure 1).
Now we proceed as shown before in the UQ post, except this time we have a separate probability density for the result at each time. Figure 2 shows the mean prediction as well as shading based on confidence intervals (at the 0.3, 0.6, and 0.9 level). The Python to generate this figure is shown below, it uses the alpha parameter (transparency) available in the matplotlib plotting package in a similar way to the vizualization approach shown in this post.
Compare the result in Figure 2 with this Monte Carlo analysis. It’s not quite apples-to-apples, because that Monte Carlo includes uncertain discovery (the input) as well. Also, this result is just a first-order collocation, so it assumes that the rate sensitivity is constant with variations in the result. This is not actually the case, and you can see if you look closely that the 0.9-interval actually dips into negative territory, which is unphysical. This is a good example of the sort of common-sense checking Hamming suggested in his “N+1” essay as a requirement for any successful computation: Are the known conservation laws obeyed by the result? [2]. Clearly we are not going to “overshoot” on extraction and begin pumping oil back into the ground.
On a somewhat related note, Hamming had another insightful thing to say about the importance of checking the correctness of a calculation: It is the experience of the author that a good theoretician can account for almost anything produced, right or wrong, or at least he can waste a lot of time worrying about whether it is right or wrong [2]. Don’t worry,verify!
[1] Cervino, L.I., Bewley, T.R., “On the extension of the complex-step derivative technique to pseudospectral algorithms,” Journal of Computational Physics, 187, 544-549, 2003.
[2] Hamming, R.W., “Numerical Methods for Scientists and Engineers,” 2nd ed., Dover Publications, 1986.
(If you are serious about the art of scientific computing, that Dover edition of Hamming’s book is the best investment you could make with a very few bucks.)
This post is about visualizing confidence intervals. Matplotlib has some really neat capabilities that are useful in that regard, but before we get into the pictures, here’s a short digression on statistics.
I like resampling-based statistics a lot, as an engineer their practicality and intuitiveness appeals to me. It has been shown that students learn better and enjoy statistics more when they are taught using resampling methods. Here’s a nice description of the motivation for resampling (ok, maybe it's a bit over the top):
For more than a century the inherent difficulty of formula-based inferential statistics has baffled scientists, induced errors in research, and caused million of students to hate the subject.Complexity is the disease. Resampling (drawing repeated samples from the given data, or population suggested by the data) is a proven cure. Bootstrap, permutation, and other computer-intensive procedures have revolutionized statistics. Resampling is now the method of choice for confidence limits, hypothesis tests, and other everyday inferential problems.
In place of the formidable formulas and mysterious tables of parametric and non-parametric tests based on complicated mathematics and arcane approximations, the basic resampling tools are simulations, created especially for the task at hand by practitioners who completely understand what they are doing and why they are doing it. Resampling lets you analyze most sorts of data, even those that cannot be analyzed with formulas.
One of the useful things to do is resampling of residuals. The assumptions underlying most models are that the magnitude and sign of the residuals are not a function of the independent variable. This is the hypothesis which we’ll base our resampling on. First we fit a model (a line in these examples) to some data, then calculate all the residuals (difference between the data and the model). Then we can apply a couple of different resampling approaches towards understanding the confidence intervals.
A permutation-based way of establishing confidence intervals is easily accomplished in Python using the itertools module’s permutation function. This is an exact method, but the number of permutations grows as the factorial of the sample size. The six-sample example shown in figure 1 has 6! = 720 possible permutations of the residuals. With only ten samples the number of permutations grows to 3628800 (over three million). The point of this post is visualization, so plotting three million lines may not be worth our while.
Figure 1 shows the least-squares-fit line in solid black, with the lines fit using the permuted residuals slightly transparent. Here’s the Python to accomplish that:
We’ve used the alpha argument to set the level of transparency in each of the permutations. Where more of the lines overlap, the color is darker. This gives a somewhat intuitive display of the ’weight’ of the variation in the model fit (confidence).
As mentioned above, doing exact permutation-based statistics becomes computationally prohibitive rather quickly, so we have to go to the random-sampling based approximations. The bootstrap is one such approach.
Rather than generating all possible permutations of the residuals, the bootstrap method involves drawing random samples (with replacement) from the residuals. This is done in Python easily enough by using the random_integer function to index into our array of residuals.
Figure 2 shows the method applied to the same six data points as with the permutation.
Since the sampling in the bootstrap approach is with replacement we can get sets of residuals that have repeated values. This tends to introduce a bit of bias in the direction of that repeated residual. You can see that behaviour exhibited in Figure 2 by the spreading of the lines around the middle of the graph, in contrast to the tight intersection in the middle of Figure 1.
Of course the reason the bootstrap method is useful is we often have large samples that are impractical to do with permutations. The Python to generate the 12-sample example in Figure 3 is shown below.
These sorts of graphs probably won’t replace the standard sorts of confidence intervals (see the graph in this post for example), but it’s a kind of neat way of looking at things, and a good demo of some of the cool stuff you can do really easily in Python with Matplotlib and Scipy.
The explicit recognition of its [the FFT’s] importance is one of the few significant advances in numerical
analysis in the past few decades. Much has been learned about various aspects of computation, but no other
recent discovery has so profoundly affected so many different fields as has the fast Fourier transform.
–R.W. Hamming [1]
This tutorial will demonstrate Gaussian convolution / deconvolution and Abel inversion of something resembling microwave interferometry data. The tutorial uses Scipy [2], but the concepts (as well as most of the function names and even the underlying FFT libraries) transfer directly to other environments (Matlab, Octave, etc). When in doubt, Google or Hutchinson [3] are your friends.
The Microwave Interferometer beam is assumed to have a Gaussian distribution of intensity. The equation for a Gaussian distribution is
 | (1) |
This means that the measurement at each spatial position is a “mixture” of effects from neighboring areas along the beam’s axis. To understand this, we’ll do some convolutions with Gaussians. First off, we need an array of points from a Gaussian, but we need them in a particular order. Here’s a Python function which does what we need:
The reason for reordering is because we’ll be using FFTs to do the convolutions, and the FFT assumes periodic data. Figure 1 shows the results for several values of the standard deviation.
Once we have the proper distributions defined, the convolution just requires an FFT of the theoretical curve, and the Gaussians, point-wise multiplication in the transformed-domain, followed by an inverse transform. The Python to do this is shown below:
Figure 2 shows a theoretical phase shift curve for a calibration cylinder along with convolutions of that same curve with Gaussians of various standard deviations.
Convolution is easy enough, but what we have in the case of an interferometry measurement is the convolved signal, and we’d like the ’true’ measurement. We can get this by deconvolution. This is simply division in the frequency space as opposed to multiplication:
Figure 3 shows the results for a couple test σ’s. Unfortunately, deconvolution is numerically ill-posed and very sensitive to noise, it tends to amplify any noise present. Unless we are very lucky and guess the exact σ and our measurement has no noise, we end up with ’wiggles’.
We can fix things by first ’smoothing’ with a small Gaussian, and then performing some test deconvolutions.
The deconvolution that uses σ > σtrue still displays overshoots even with the smoothing, it also becomes too ’sharp’. The calibration measurements can be used to find the optimal smoothing and sharpening σ for the particular interferometer / frequencies being used.
Many of the same tools presented in Section 1 can be extended to the calculation of the Abel inversion. In microwave interferometry we have a ’chord-averaged’ measurement, but we would like to understand the radial distribution of the quantity (phase shift / plasma density).
 | (2) |
where F is our chord-averaged measurement, and f(r) is the radial distribution that we would like to recover. To calculate the inverse transform we first need an estimate of the derivative ∂F∕∂y. Like deconvolution, numerical differentiation tends to amplify noise (numerical or otherwise). To avoid this we can convolve the derivative of a Gaussian with the function we wish to differentiate and then perform an inverse transform to recover a ’smoothed’ numerical derivative. The derivative of the Gaussian is
 | (3) |
A Python function that gives a properly re-ordered first derivative is
In the same way that we used the Gaussian for smoothing we can use its derivative for differentiation. Here’s the Python to do that:
Figure 5 shows the results. The wiggles that we don’t manage to smooth in the initial sharpening stage become amplified by differentiation.
Once we have a numerical approximation for the derivative, calculating the integral is straight-forward:
The one minor ’gotcha’ is dealing with the integrable singularity in 2. We just copy the value of the nearest neighbor over the Inf, the integral will approach the correct value in the limit of small Δx (lots of points).
The analytical inverse Abel transform of the theoretical phase shift curve can be easily calculated (Maxima [4] script shown in section A.2). The theoretical phase shift has the form
 | (4) |
which gives the derivative as
 | (5) |
and finally, the inverse Abel transform as
 | (6) |
Often, despite our best efforts, the data are not centered. There is a method to estimate the optimal shift required to center the data using the FFT [5]. Figure 7 shows one of the “measurements” shifted by a small amount.
The method of Jaffe et al. [5] is to maximize
![N∕2-1 ∑ [ -2πmk Δfn ]2 Re(e Lnh) -N ∕2](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgey9LM9xfdMaMFFvvbB04nDHSHO_-ExftZuSHvzZsFrCXgXiBtRokrFGpJWrcVxGiYDun0bgzrIN5JFgE6NqfwC5YrlTe1jB0mbBTlza_ON4ncYOieht3pjtMcpejqQI6FKsI5LHl_pLM/s1600-h/abel6x.png) | (7) |
Where Lnh is the Fourier transform of the data, and mkΔ is a shift. Since the data is real, the real part of its transform is the even part of the signal. The shift which maximizes the magnitude of the even part is the one which makes the data most even (centered). A graphical analysis of the data can provide a range of shifts to try. The Python to perform this maximization is shown below (this relies on being able to index arrays with arrays of integers, which can be done in Matlab as well).
Figure 8 shows that the objective function (equation 7) reaches a local maximum at the shift which was originally applied (shown in Figure 7).
The techniques presented in section 1 can be combined to fit a calibration model to interferometry measurements of a calibration cylinder (or any other object for which you have the theoretical phase-shift curve). The purpose of a calibration model is to provide estimates of parameters like the standard deviation of the Gaussian beam for use in deconvolution and inversion of subsequent test measurements.
There steps in fitting the model are
The Python functions to accomplish this are shown below. One error function uses a sum of squared errors and the other uses a maximum error.
The error in the fitted model following this approach is shown in Figure 9. Both the least squares error and the min-max error methods recover the correct deconvolution σ.
Jaffe et al. [5] also suggest “symmetrizing” before calculating the Abel transform by removing the imaginary component of the transformed data. They also estimate the variance of the inverse Abel transform based on both the portion of the data removed by filtering (smoothing) and symmetrizing.
Rather than calculating the derivative necessary for the inverse Abel transform with a convolution, a pseudo-spectral method could be used. The trade-off being that very accurate derivatives could be attained at the expense of greatly amplifying any noise. A more practical variation, that is not susceptible to noise, would be to calculate the integral with a pseudo-spectral approach rather than the composite trapezoidal rule. This would give a more accurate integral, and bring the numerical and analytical curves in Figure 6 into closer agreement.
The complete Python script to do all of the calculations and generate all of the plots in this tutorial is available.
The Maxima batch file used to calculate the inverse Abel transform of the theoretical phase-shift curve is available.
[1] Hamming, R.W., “Numerical Methods for Scientists and Engineers
[2] Jones, E., Oliphant, T., Peterson, P., “Scipy: Open source scientific tools for Python,” 2001–, http://www.scipy.org
[3] Hutchinson, I.H., “Principles of Plasma Diagnostics
[4] Schelter, W., et al., “Maxima: a computer algebra system,” http://maxima.sourceforge.net
[5] Jaffe, S.M., Larjo, J., Henberg, R., “Abel Inversion Using the Fast Fourier Transform,” 10th International Symposium on Plasma Chemistry, Bochum, (1991).
