User:Sean OByrne
Sean O'Byrne
Email: <<MailTo(sbobyrne AT SPAMFREE tpg DOT com DOT au)>>
Bayesian statistics demonstration
Here is a little demonstration program, based upon one of the early examples in the textbook 'Data Analysis: A Bayesian Tutorial' by D.S. Sivia. It may be useful for anyone who, like me, is just starting to learn Bayesian statistics. It's an example involving determining whether a coin is biased towards heads or tails, and it shows how the degree of certainty about the result changes with the number of tests.
Where to sample a Gaussian curve to get the best fit
This is an interesting problem arising from my research. One part of my work involves measuring velocities using fluorescence signals. The Doppler shift of the fluorescence peak is measured by exciting fluorescence in a gas using a laser tuned to an energy level transition of nitric oxide in a flow. The fluorescence intensity has a Gaussian distribution as wavelength is changed. The speed of the flow is related to the position of the peak of this Gaussian distribution.
As this is an experiment, every measurement will have an uncertainty associated with it, and several such experiments are required to determine the position of the Gaussian centre. Typically, the domain of possible wavelengths is split into relatively equally-spaced regions, and samples are taken within that region. This is not optimal, because a number of the points will be outside the curve, where the signal says little or nothing about the peak location.
This code is my attempt to make an automated algorithm to pick the best place to sample the Gaussian such that it will give the most information about the peak location. The background and algorithm behind the program is described in this conference paper:
I am not particularly happy with the code I wrote for this. It works, but it could be a lot better written. In particular, I found it difficult to interface with the nlls nonlinear least-squares routine, and had to use a few globals, which became lots of globals, combined with for loops etc. I'd appreciate anyone taking the effort to improve the code and make it more J-like. I have a lot to learn before my code could be considered elegant.
Here is the code...
and here is someone else's nonlinear least-squares routine that I used. File:Nlls j602 2 jan 2010.ijs
The code can be run by typing n GaussFill, where n is the number of additional points you want to adaptively sample in the range _10 to 10. The algorithm starts with values at _10,_5,0,5,10. If it's working OK, the points chosen will tend to cluster around the sides of the Gaussian. In its default form, the Gaussian has 5% noise added. The verb UniFit models the 'dumb' case of uniformly distributed points for comparison with the adaptive solver,
While I think I solved my problem, I can't help but think that there must be better and more elegant ways of solving this problem. It's really a very simple problem to pose, but not so simple to solve. There is an alternative way of doing this, called Bayesian adaptive exploration, but I don't know if I could implement it based upon the description in the paper presenting it. Once again, if you have a better way of doing this than my primitive method, I'd be very interested to hear about it.