Modern governments are gung ho to privatize just about everything they can lay their hands on these days. In the '80s, the Thatcher government was busy selling the utilities. Her successor sold the railways. Here in Canada we've begun to open up the phone services to competition and have sold the national airline. And in the states there are experiments to privatize the prisons, the police force and much of the military.
One thing that has remained stubbornly in the public domain is the currency. And by this, I don't mean the creation of money, which now too, is largely in the private sector. I mean the currency itself. If we are not happy with our dollar and what it buys and its value (which, after all, is nothing--it's just a piece of paper, or in these times a blip in some computer data bank somewhere), why shouldn't we exchange it for another one, or for that matter for a Canadian Franc, or Dragma or Schekel or Pound or Ob or Guinea or whatever it is that the backer decides to call it?
Lets face it, competition for the bank conglomerates is sorely needed. These days, because of fees, abysimally low interest rates and other factors, if your net worth isn't close to a cool million, you're better off keeping your money in your shoe. And as the WikiLeaks debacle proves (in addition to the fact that freedom of the press is a convenient myth) the banks now wield an excessive amount of power.
What we need is a new standard of value. If the "Occupy" Movement is really serious about achieving its goals, two things need to be done.
The first thing that we need to do is organize a general strike. Nobody does any work, and by work I mean any revenue-generating activity. I'm not saying we need to go hungry. You can do anything you want to, as long as it doesn't involve cash. So for instance, you can bring your neighbours food or fix their cars for them and so on.
It is now well understood how fiat currency created primarily through debt helps concentrate wealth in the hands of the few--primarily those who already have control over a considerable amount of wealth. The more the wealth is concentrated, the more control the global elite ultimately possess. This paper shows just how concentrated that power has now become. So the second thing the Occupy Movement needs--and this compliments the first--is a new standard of value. Not necessarily currency since any new currency will be in immanent danger of being usurped for the purposes of control by an elite few just like what has happened to the existing currency.
My own vision for a better society and a better system of distribution ("market") is simple. We just need a clearing-house or bulletin board for work that needs to be done and goods that need to be distributed, whether for ourselves or for our neighbours. Sort of like Craigslist except everything in the want-ads is free. If I am free and I possess the skills, then I can do one of the tasks. Likewise, if my neighbour is free, he can come do one of the tasks I've posted on the board. If I have some spare bread that I need to get rid of, I just post it on the board.
Here is a book that discusses strategies for alternative currency. I think it is one of the most important books written in recent times.
Thursday, May 16, 2013
Friday, May 10, 2013
Celebrations are in order
I stated in the previous post that all of the work from the past year-and-a-half was done for free. This is not entirely true. First, I did received 1500 Euros for a stale contract. Two of the papers on arxiv are offshoots of that work. But more importantly, I received a $10 donation quite recently. This is actually my first donation, and it's still sitting in my Paypal acount--or at least what's left of it after the fees. Once a few more donations start rolling in I intend to divert them to a dedicated account.
Why did I not announce it sooner? First of all, and I don't know if this says something about me or about society in general or perhaps both, but when I first received the e-mail, I assumed it was a scammer. I even filled out an electronic report to the RCMP about it! It turns out my Paypal account hadn't been verified. Second, after almost two years of dedicated work as an independent research scientist this is all I have to show for it, at least monetarily. I realize the necessity, passed down to us from all the great moral teachers, of gratitude. But I am somewhat embarassed about this.
Anyway, thanks very much to William Oquendo for your generous donation to the libAGF project!
Why did I not announce it sooner? First of all, and I don't know if this says something about me or about society in general or perhaps both, but when I first received the e-mail, I assumed it was a scammer. I even filled out an electronic report to the RCMP about it! It turns out my Paypal account hadn't been verified. Second, after almost two years of dedicated work as an independent research scientist this is all I have to show for it, at least monetarily. I realize the necessity, passed down to us from all the great moral teachers, of gratitude. But I am somewhat embarassed about this.
Anyway, thanks very much to William Oquendo for your generous donation to the libAGF project!
Jumping on the crowd-funding bandwagon--again
Recently I've become interested in "men's issues" because I believe that modern feminism has gone too far. This is particularly true in academia where intellectual freedom (assuming it ever existed at all) has practically disappeared. A Harvard president is asked to resign after he merely suggests that men and women may have different abilities, meanwhile, academic feminists can write the most blatant hate-speech against men and other, so-called, "privileged" groups without fear of repercussion.
One of the characters I've encountered through men's rights activists is Anita Sarkeesian. She is a feminist "vlogger" who submitted a project to Kickstarter.com for the purpose of examining female stereotypes in video games. She collected over $158 000.
That's a lot of money for something so trivial. Women are stereotyped in videogames? Get over yourself. As if men aren't stereotyped as well. It seems inevitable in a medium where the characters are composed of pixels that a certain amount of cartoonishness is going to creep in. Yes, your female avatar may have an unnaturally large chest and be rather scantily dressed, but so too, in all probability, will your male avatar. Big deal.
Anyway, this has inspired me to try my hand once again at collecting donations using crowdfunding. This time, however, I'm going to use generic sites, rather than science specific ones. I'm also going to focus on those that use a "keep-it-all" funding model.
In the past year-and-a-half alone, I have made over 250 code commits to free software projects, have submitted a total of five articles to the arxiv.org database, as well as dozens of contributions to Wikipedia. These activities represent advances in diverse fields including statistics, machine-learning, numerical computation, atmospheric physics, remote sensing and environmental science. Little of which has yet been paid for. If something as frivolous as tropes in video games can raise over six figures, surely real science deserves at least a small nugget.
So give me money. You owe it to me. Or at least somebody does...
One of the characters I've encountered through men's rights activists is Anita Sarkeesian. She is a feminist "vlogger" who submitted a project to Kickstarter.com for the purpose of examining female stereotypes in video games. She collected over $158 000.
That's a lot of money for something so trivial. Women are stereotyped in videogames? Get over yourself. As if men aren't stereotyped as well. It seems inevitable in a medium where the characters are composed of pixels that a certain amount of cartoonishness is going to creep in. Yes, your female avatar may have an unnaturally large chest and be rather scantily dressed, but so too, in all probability, will your male avatar. Big deal.
Anyway, this has inspired me to try my hand once again at collecting donations using crowdfunding. This time, however, I'm going to use generic sites, rather than science specific ones. I'm also going to focus on those that use a "keep-it-all" funding model.
In the past year-and-a-half alone, I have made over 250 code commits to free software projects, have submitted a total of five articles to the arxiv.org database, as well as dozens of contributions to Wikipedia. These activities represent advances in diverse fields including statistics, machine-learning, numerical computation, atmospheric physics, remote sensing and environmental science. Little of which has yet been paid for. If something as frivolous as tropes in video games can raise over six figures, surely real science deserves at least a small nugget.
So give me money. You owe it to me. Or at least somebody does...
Friday, April 12, 2013
A sparse matrix calculator
For the past few months I've been working quite hard on both the ctraj and "libpetey" software. I hope to have new releases of both within the next month.
One of the most significant additions is a sparse matrix calculator. This tool allows you to perform operations on sparse matrices much as one might do to simple floating point numbers on a pocket calculator or a calculator app. This addition is so significant, in fact, that I'm consider creating a whole new project to contain it.
Two of the "peteysoft" projects already output sparse matrices. The Gaspard-Rice simulator, gr_scat, can output a sparse matrix representing a discretization of the Laplace or Helmholtz equation which models the wave scattering behaviour of the system. Also, the ctraj tracer simulation outputs a series of sparse matrices that represent the tracer dynamics. To model a tracer, the initial tracer configuration is multiplied through with these matrices. Storing them as an array is both more space and time efficient than multiplying them through.
I have already used the matrices output from ctraj in a quite sophisticated analysis: by correlating the singular vectors with a series of sparse measurements (e.g. from a satellite), they can be used in a type of dynamical tracer reconstruction I call "PC proxy." I submitted a paper about PC proxy analysis to Atmospheric Environment--it can be found here. The genesis of the sparse calculator project comes about because the reviewers were interested in learning in more detail how the method works, what makes it better than existing methods and how better to parametrize it. Since I don't really understand the method all that well myself, I decided to write this software tool to conduct a series of numerical experiments.
Sparse matrics can be used to model a large class of partial differential equations. In the past, the process of matrix multiplication would be coded directly in Fortran. A similar method is used to compute the tangent model for generating singular vectors from numerical weather models. Each time a matrix multiplication is desired, to be fed back to the eigenvalue subroutines, the model is re-run from scratch. Porting them to a data structure instead is not only faster, it is far more versatile. Since the matrix is typically quite sparse, the extra storage is not particularly burdensome.
The basic idea behind the tool is that files generated from other programs are imported and used as variables. Because these files can become quite large, the matrices are only loaded into RAM when they are used in a calculation. Except for scalars, variable names always correspond to a file name.
Using the ctraj package, we can generate a series of sparse matrices which represent the tracer dynamics for a given range of dates, say May 1, 2000 to May 5, 2000, and a given vertical level, say 500 K. First, you will need to download the NetCDF files containing NCEP temperature and wind data: air.2000.nc, uwnd.2000.nc, vwnd.2000.nc. Next, extract and convert the desired velocity data using the nc2ds command:
> nc2ds -i 2000/5/1 -f 2000/5/6 -T 500 vfield_n.ds
> nc2ds -- -i 2000/5/1 -f 2000/5/6 -T 500 vfield_s.ds
Next, we run the tracer model on this data:
> ctraj_tracer vfield_s.ds vfield_n.ds tracer.dat
Now we can find the top 10 singular vectors of the tracer dynamics using sparse_calc and put them in a file called sv.vec:
> echo "sv.vec=svd(sparse_array(tracer.dat)*transpose(tracer.dat), 10)"\
| sparse_calc -A 50
Finally, plot and display the first singular vector:
> plot_frame.sh 0 sv.vec sv0.ps
> gv sv0.ps
Now lets say you want to do the same thing, but only for the first day's
worth of data:
> echo \
"sv_1day.vec=svd(sparse_array(tracer.dat(0:4))*transpose(tracer.dat(0:4)), 10)"\
| sparse_calc -A 50
The sparse matrix calculator will be available in the next release of libpetey , currently a part of the msci project. The tool can handle generalized sparse matrices, sparse matrix arrays, full matrices and simple floating point calculations. There are routines for matrix solution, inversion and eigenvalue calculation. In addition, there is a generous online help function as well as a man page to quickly get you up to speed using it.
One of the most significant additions is a sparse matrix calculator. This tool allows you to perform operations on sparse matrices much as one might do to simple floating point numbers on a pocket calculator or a calculator app. This addition is so significant, in fact, that I'm consider creating a whole new project to contain it.
Two of the "peteysoft" projects already output sparse matrices. The Gaspard-Rice simulator, gr_scat, can output a sparse matrix representing a discretization of the Laplace or Helmholtz equation which models the wave scattering behaviour of the system. Also, the ctraj tracer simulation outputs a series of sparse matrices that represent the tracer dynamics. To model a tracer, the initial tracer configuration is multiplied through with these matrices. Storing them as an array is both more space and time efficient than multiplying them through.
I have already used the matrices output from ctraj in a quite sophisticated analysis: by correlating the singular vectors with a series of sparse measurements (e.g. from a satellite), they can be used in a type of dynamical tracer reconstruction I call "PC proxy." I submitted a paper about PC proxy analysis to Atmospheric Environment--it can be found here. The genesis of the sparse calculator project comes about because the reviewers were interested in learning in more detail how the method works, what makes it better than existing methods and how better to parametrize it. Since I don't really understand the method all that well myself, I decided to write this software tool to conduct a series of numerical experiments.
Sparse matrics can be used to model a large class of partial differential equations. In the past, the process of matrix multiplication would be coded directly in Fortran. A similar method is used to compute the tangent model for generating singular vectors from numerical weather models. Each time a matrix multiplication is desired, to be fed back to the eigenvalue subroutines, the model is re-run from scratch. Porting them to a data structure instead is not only faster, it is far more versatile. Since the matrix is typically quite sparse, the extra storage is not particularly burdensome.
The basic idea behind the tool is that files generated from other programs are imported and used as variables. Because these files can become quite large, the matrices are only loaded into RAM when they are used in a calculation. Except for scalars, variable names always correspond to a file name.
Using the ctraj package, we can generate a series of sparse matrices which represent the tracer dynamics for a given range of dates, say May 1, 2000 to May 5, 2000, and a given vertical level, say 500 K. First, you will need to download the NetCDF files containing NCEP temperature and wind data: air.2000.nc, uwnd.2000.nc, vwnd.2000.nc. Next, extract and convert the desired velocity data using the nc2ds command:
> nc2ds -i 2000/5/1 -f 2000/5/6 -T 500 vfield_n.ds
> nc2ds -- -i 2000/5/1 -f 2000/5/6 -T 500 vfield_s.ds
Next, we run the tracer model on this data:
> ctraj_tracer vfield_s.ds vfield_n.ds tracer.dat
Now we can find the top 10 singular vectors of the tracer dynamics using sparse_calc and put them in a file called sv.vec:
> echo "sv.vec=svd(sparse_array(tracer.dat)*transpose(tracer.dat), 10)"\
| sparse_calc -A 50
Finally, plot and display the first singular vector:
> plot_frame.sh 0 sv.vec sv0.ps
> gv sv0.ps
Now lets say you want to do the same thing, but only for the first day's
worth of data:
> echo \
"sv_1day.vec=svd(sparse_array(tracer.dat(0:4))*transpose(tracer.dat(0:4)), 10)"\
| sparse_calc -A 50
The sparse matrix calculator will be available in the next release of libpetey , currently a part of the msci project. The tool can handle generalized sparse matrices, sparse matrix arrays, full matrices and simple floating point calculations. There are routines for matrix solution, inversion and eigenvalue calculation. In addition, there is a generous online help function as well as a man page to quickly get you up to speed using it.
Thursday, February 14, 2013
Additional thoughts on software design: keeping it simple
There are some aspects of software design that are quite difficult to teach and even more difficult to learn. Like elegance in design. I remember helping a colleague with a subroutine to randomly permute a list--just like shuffling a deck of cards. He had this enormous subroutine that generated random numbers then checked them against other numbers. I just shook my head. It's counter-intuitive--at least until you figure it out--but the best way I've found is to generate a list of random numbers, sort them, and then use the shifts so generated to permute the list. In the language we were working with (Interactive Data Language) it's a one-liner.
I don't remember having trouble coming up with this simply because I've worked with simulated random numbers for so long. I've also used sorting for a unique application: calculation of equivalent latitude. Check out the Wikipedia page (which I wrote). I think the algorithm for calculating equivalent latitude is the epitome of elegance and I've discovered a number of very similar tricks for quite different applications.
Not to say that I haven't created my share of bloated code. I remember once writing a piece of code to find local maxima in a gridded, two dimensional field. What I came up with involved dividing the data values into several ranges and using a domain-filling algorithm (not of my own design). Not very accurate and even less fast. A simple, gradient ascent algorithm is faster by far (I know because I eventually rewrote it) or what about just checking for a local maximum in the eight nearest-neighbours?? I recall having tried both at the time, but couldn't get either working. That's the other lesson: sometimes it takes patience to ferret out all the bugs in a piece of code. It also takes confidence in believing that a particular approach will work.
Another particularly embarrassing piece of bloatware I've written involved binning measurement locations from satellite data in the process of interpolating them to a regular grid. The program did work quite well, but I'm pretty sure the binning actually made the algorithm slower. I'm not certain of that since I haven't rewritten it yet. Meanwhile, the indices it generated were gigantic. (I'm still convinced that some type of binning can produce significant speed gains in nearest-neighbour type problems, but I haven't yet been able to come up with the goods.)
I don't remember having trouble coming up with this simply because I've worked with simulated random numbers for so long. I've also used sorting for a unique application: calculation of equivalent latitude. Check out the Wikipedia page (which I wrote). I think the algorithm for calculating equivalent latitude is the epitome of elegance and I've discovered a number of very similar tricks for quite different applications.
Not to say that I haven't created my share of bloated code. I remember once writing a piece of code to find local maxima in a gridded, two dimensional field. What I came up with involved dividing the data values into several ranges and using a domain-filling algorithm (not of my own design). Not very accurate and even less fast. A simple, gradient ascent algorithm is faster by far (I know because I eventually rewrote it) or what about just checking for a local maximum in the eight nearest-neighbours?? I recall having tried both at the time, but couldn't get either working. That's the other lesson: sometimes it takes patience to ferret out all the bugs in a piece of code. It also takes confidence in believing that a particular approach will work.
Another particularly embarrassing piece of bloatware I've written involved binning measurement locations from satellite data in the process of interpolating them to a regular grid. The program did work quite well, but I'm pretty sure the binning actually made the algorithm slower. I'm not certain of that since I haven't rewritten it yet. Meanwhile, the indices it generated were gigantic. (I'm still convinced that some type of binning can produce significant speed gains in nearest-neighbour type problems, but I haven't yet been able to come up with the goods.)
Thursday, February 7, 2013
Stuff that I'm working on
I started this article some time ago, so some of the stuff is already out of date. It's quite long as well, so be warned...
While they should make the software faster and more robust, they are rather routine and somewhat boring improvements. What follows are some more interesting ideas, in order of interest and complexity. Of course, what makes an idea interesting, often also makes it complex and difficult to implement.
Another issue is tolerance of the root. This is exacerbated by the problem just discussed: very low x-tolerances can't be used. If we are using y- tolerances, then for the error function, clearly they must be varied depending on the value of x. For a root-finding algorithm, the idea of "relative tolerance" doesn't have much meaning [3].
This has now been implemented in version 0.9.6.
In the absence of binning, the time efficiency of finding a set of k-nearest neighbours is on the order n, where n is the number of samples. With binning, we can do one of two things: we can either reduce the coefficient, or we can get it down to order klog n.
In the first case, the standard, rectangular binning method in which each dimension is evenly divided into several blocks, would apply. Ideally, the bins should be large enough that most of them contain several samples, but small enough that the number of samples in each bin is significantly smaller than the total. As the number of dimensions increase, the number of bins becomes unmanageably large, while the fraction occupied also gets smaller. I have designed binning routines of this type for application over the surface of the globe. Oddly, they don't provide any speed advantage over simply running through all the samples. Either there is too much book- keeping and overhead, or there are optimizations implemented in the processor, the compiler or both that already provide a similar advantage.
An interesting solution to the multiplication of bins at higher dimensions and the presence of empty bins is to discretize each of the coordinate samples based on the bin size. Each dimension is assigned a rank and the resulting list of integer numbers sorted as you would sort any other integer. The storage requirement is typically no more than double the original. If the bin size is set so that the average number of samples per bin is very low, on the order of one, then we can build up the number of bins to check in a concentrically expanding hypersphere until we have the required number of samples. While this idea looks good on paper, the required amount of book-keeping and overhead is far more than simply searching through each of the samples one-by-one, especially since finding a given bin is now a logn operation.
Both these examples are using a somewhat naive rectangular binning strategy. A more sophisticated technique would use an adaptive binning strategy that depends on the actual distribution of data. Something like a KD-tree. I still have yet to try these out--coding one up would be far from trivial, while I find I have a distinct aversion to black-boxes so I tend to use pre-compiled libraries rather sparingly. The KD-tree also suffers from an excess of book-keeping as even searching for a single nearest neighbour requires checking every adjacent bin. The number of adjacent bins growing, of course, with the number of dimensions and often lying on very different branches of the tree. Searching fork -nearest-neighbours or all
neighbours within a radius would be even more problematic.
e2=(P/n/hD)2+(h2/n grad2 P)2
calculated respectively, in the limit as the filter width becomes very small and in the limit as it becomes very large. Here, P is the actual density function, n is the number of samples, h is the bandwidth and D is the number of dimensions. The actual density function, P, is of course, not known. But there is no reason why we cannot fit the function through repeated estimates using different bandwidths. Since e=P-f where f is the estimated density we can rearrange the above equation as follows:
n2(f-P)=(1/hD-n2)P2/f + h4 (grad2 P)2/f
Since we have several different samples, we can construct a system of equations by substituting values for f and h. The constant term, P, on the left-hand-side can be eliminated by subtracting different equations in the system:
(1/hD-n2)(1/fi - 1/fj) P2 + (hi4/fi - hj4/fj) (grad2 P)2 = n2(fi - fj)
so in theory we need a minimum of three samples to solve for P and grad2P. In practice, the method should be more stable by taking more than the minimum number of samples and performing a least-squares fit.
I recently took a second look at this project--I've had the method implemented for some time but it didn't work. I've now got it returning reasonable results for the optimal filter width, at least when it works at all, though not if you take the P2 factor from the above equation. Part of the problem is that the matrix is very ill-conditioned--the two variables we are solving for are separated by many orders of magnitude likewise all the different elements of the matrix.
The reviewers were interested first, that the manuscript be fleshed out somewhat, and second, how better to determine the parameters of the method: how many principal components to use and how long of a lead time to use. I am attempting something more general here: to build a simple sparse matrix calculator so that I can perform a series of numerical experiments with the matrices generated by the tracer simulation. I'm hoping that will help me to better understand just how the method works: for instance, how is it possible that you can correlate measurements from before the matrix was integrated and the method will still give good results?
The sparse matrix calculator has already been implemented and has even passed a few very basic tests. I still need to test it with actual, practical matrices like those from the tracer simulation, and of course it needs many refinements and features to make it complete, but the basic idea is there and it's running. At the moment I'm working on refactoring the entire trajectory suite and I'd like to finish that first before continuing with this.
There are two major ways you can divide up the classes: in the first case so that all of the conditional probabilities can be solved for and in the second case so that only the conditional probability of the "winning" class is returned. In the first case, you need at least (n-1)/2 borders--that is, for each border you can write two equations relating the conditional probabilities plus one more for the entire system (they sum to unity.) In the second case, the method is hierarchical: divide the set of classes into two groups, divide each pair of groups into two more and so on. Consider the case of dividing up four classes:
x | o
-----
+ | *
The four classes are x, o, + and *. In the first case, we think of this diagram as representing two borders: one horizontal, and one vertical. In the second case, consider first a horizontal border, and then two vertical borders. It stands to reason that the two methods can also be combined.
[2] Donald Knuth (1998) The Art of Computer Programming, Volume 3: Sorting and Searching, 2nd edition Addison-Wesley.
[3] William Press, Saul Teukolsky, William Vetterling, Brian Flannery (1992). Numerical Recipes in C, 2nd edition Cambridge University Press, Cambridge.
[4] Peter Mills (2012). "Principal component proxy tracer analysis." arxiv:1202.2194
[5] Peter Mills (2006), "The influence of noise on a classical chaotic scatterer." Communications in Nonlinear Science and Numerical Simulation, 11 (8) 899-906.
New software release and planned improvements
The Peteysoft Foundation has released a new version of the libagf software and the libpetey scientific library required by it. Changes are fairly minor, but nonetheless important. The Runge Kutta integrator now has a thread-safe version. The function pointer passed to this routine is compatible with the GSL. New routines have been introduced including a symbol table class and a "smart" pointer class. More significant improvements are detailed in the next two sections, while many exciting planned new developments follow.k-least algorithm
The kleast subroutine selects the k smallest elements in an array. Past versions of this algorithm were based on a binary tree and have a theoretical average time efficiency of O(n log k), although tests suggest that the actual running time is slightly better, closer to n+k log k [1]. This is the theoretical best efficiency for a k-least algorithm that sorts the resulting output, necessarily true for a binary tree. According to Knuth, if the output is not sorted, the efficiency can be as good as O(n) [2]. This can be acheived through a modified quicksort algorithm, but not one based on a heap which again, sorts the output. A new, faster k-least algorithm based on a quicksort has been implemented and applied to the libagf library to pre-select the data and to find the k-nearest- neighbours.Solving for the filter width
The old algorithm for solving for the filter width while holding the total of the weights, W, was based on repeated squaring of the weights until the total passed the target value. The final result was found by exponential interpolation between the final value and the previous value. The purpose of this excessively clever method was to reduce the number of computed exponentials--the number is fixed at two per weight. Unfortunately, it was quite finicky and improper initialization or certain sample distributions (such as duplicate samples) would cause the method to hang. To address these issues, it was replaced with nothing more than a root-finding algorithm, although a fairly sophisticated one. This is the same as used for searching out the class borders: a third order interpolation method which brackets the roots. I call this root-finding method, "supernewton". [1]While they should make the software faster and more robust, they are rather routine and somewhat boring improvements. What follows are some more interesting ideas, in order of interest and complexity. Of course, what makes an idea interesting, often also makes it complex and difficult to implement.
The inverse error function and the superNewton algorithm
The inverse error function is quite useful since it tells you for a given integrated probability and Guassian statistics, how far from the centre of the distribution you need to be. So for instance, if you have n normally distributed random samples, you can use the inverse error function to figure out approximately where the most distant sample should lie. Unfortunately, most math libraries don't include a numerical implementation of it. I figured, no problem, I'll use my powerful numerical inverse method "superNewton" to implement my own. Rather than giving me a simple, efficient math function in record time, it revealed to me some rather striking limitations in the algorithm as it is currently implemented. Most importantly, it tends to get stuck always re-bracketing the the root on the same side. This could probably be fixed quite easily, but at the cost of some extra book-keeping. Simply track how many times the estimated root has landed on the same side of the actual root as well as the relative distance between the previous left and right brackets. If these values are too high, the estimated root can be shifted so that it's now on the opposite side of the actual root. I would suggest doubling the difference between the estimated root and the closest bracket.Another issue is tolerance of the root. This is exacerbated by the problem just discussed: very low x-tolerances can't be used. If we are using y- tolerances, then for the error function, clearly they must be varied depending on the value of x. For a root-finding algorithm, the idea of "relative tolerance" doesn't have much meaning [3].
Cross-validation of PDF estimates
Cross-validation or n-fold cross validation are simple and reliable ways of testing the accuracy of a supervised statistical classification algorithm. Something similar would be nice to have for probability density function (PDF) estimates as well. Unfortunately, unlike that for a classification problem, actual values of the PDF are not contained in training data for a probability estimation problem and the underlying PDF is frequently not known. This can be solved quite readily by adding an extra step. A new dataset with the equivalent distribution as the existing training data is first simulated using the same method as a metropolis Monte Carlo [3]. This dataset is not used for PDF estimation, but rather the original dataset is randomly split up, PDFs are calculated for each portion at the locations in the simulated dataset and then compared, e.g. through cross-correlation. Note that the parameters used (k for KNN, W for AGF) will need to scaled up in proportion (i.e. by n if this is considered a form of n-fold validation) when using the full dataset.This has now been implemented in version 0.9.6.
Binning
There are so many applications for an efficient, multi-dimensional binning algorithm, I don't even know where to begin. For instance, in the surgery step of contour advection with surgery, a nearest neighbour search must be performed for every pair of points to check if they are too close or not. Kernel estimators are more efficient if the samples used are restricted to a subset near the centre of the kernel function. Binning can also be used for collocation routines used for the purpose, for instance, of creating satellite inverse methods.In the absence of binning, the time efficiency of finding a set of k-nearest neighbours is on the order n, where n is the number of samples. With binning, we can do one of two things: we can either reduce the coefficient, or we can get it down to order klog n.
In the first case, the standard, rectangular binning method in which each dimension is evenly divided into several blocks, would apply. Ideally, the bins should be large enough that most of them contain several samples, but small enough that the number of samples in each bin is significantly smaller than the total. As the number of dimensions increase, the number of bins becomes unmanageably large, while the fraction occupied also gets smaller. I have designed binning routines of this type for application over the surface of the globe. Oddly, they don't provide any speed advantage over simply running through all the samples. Either there is too much book- keeping and overhead, or there are optimizations implemented in the processor, the compiler or both that already provide a similar advantage.
An interesting solution to the multiplication of bins at higher dimensions and the presence of empty bins is to discretize each of the coordinate samples based on the bin size. Each dimension is assigned a rank and the resulting list of integer numbers sorted as you would sort any other integer. The storage requirement is typically no more than double the original. If the bin size is set so that the average number of samples per bin is very low, on the order of one, then we can build up the number of bins to check in a concentrically expanding hypersphere until we have the required number of samples. While this idea looks good on paper, the required amount of book-keeping and overhead is far more than simply searching through each of the samples one-by-one, especially since finding a given bin is now a log
Both these examples are using a somewhat naive rectangular binning strategy. A more sophisticated technique would use an adaptive binning strategy that depends on the actual distribution of data. Something like a KD-tree. I still have yet to try these out--coding one up would be far from trivial, while I find I have a distinct aversion to black-boxes so I tend to use pre-compiled libraries rather sparingly. The KD-tree also suffers from an excess of book-keeping as even searching for a single nearest neighbour requires checking every adjacent bin. The number of adjacent bins growing, of course, with the number of dimensions and often lying on very different branches of the tree. Searching for
Optimal AGF
The current version of Adaptive Gaussian filtering is quite useful and has a reasonable justification. We can, however, write an expression for the error of kernel density estimator as a function of the bandwidth. There are two terms, a variance term and and a bias term:e2=(P/n/hD)2+(h2/n grad2 P)2
calculated respectively, in the limit as the filter width becomes very small and in the limit as it becomes very large. Here, P is the actual density function, n is the number of samples, h is the bandwidth and D is the number of dimensions. The actual density function, P, is of course, not known. But there is no reason why we cannot fit the function through repeated estimates using different bandwidths. Since e=P-f where f is the estimated density we can rearrange the above equation as follows:
n2(f-P)=(1/hD-n2)P2/f + h4 (grad2 P)2/f
Since we have several different samples, we can construct a system of equations by substituting values for f and h. The constant term, P, on the left-hand-side can be eliminated by subtracting different equations in the system:
(1/hD-n2)(1/fi - 1/fj) P2 + (hi4/fi - hj4/fj) (grad2 P)2 = n2(fi - fj)
so in theory we need a minimum of three samples to solve for P and grad2P. In practice, the method should be more stable by taking more than the minimum number of samples and performing a least-squares fit.
I recently took a second look at this project--I've had the method implemented for some time but it didn't work. I've now got it returning reasonable results for the optimal filter width, at least when it works at all, though not if you take the P2 factor from the above equation. Part of the problem is that the matrix is very ill-conditioned--the two variables we are solving for are separated by many orders of magnitude likewise all the different elements of the matrix.
Sparse matrix calculator
I have recently put out a paper that describes a powerful dynamical tracer reconstruction method called principal component proxy analysis [4]. The method is based on decomposing a matrix representation of the tracer dynamics using principal component analysis (PCA). Each PC is then regressed with a series of (typically sparse) tracer measurements and the entire field thus reconstructed. The manuscript was sent to Atmospheric Environment but was unfortunately rejected.The reviewers were interested first, that the manuscript be fleshed out somewhat, and second, how better to determine the parameters of the method: how many principal components to use and how long of a lead time to use. I am attempting something more general here: to build a simple sparse matrix calculator so that I can perform a series of numerical experiments with the matrices generated by the tracer simulation. I'm hoping that will help me to better understand just how the method works: for instance, how is it possible that you can correlate measurements from before the matrix was integrated and the method will still give good results?
The sparse matrix calculator has already been implemented and has even passed a few very basic tests. I still need to test it with actual, practical matrices like those from the tracer simulation, and of course it needs many refinements and features to make it complete, but the basic idea is there and it's running. At the moment I'm working on refactoring the entire trajectory suite and I'd like to finish that first before continuing with this.
Multi-class classification with AGF borders
If you have a discrimination border, this is only good for discriminating two classes. If you want to be more discriminating, you need more than one border. This is something I've been planning for a long time, but made very little progress on. The problem is: there is no one-size-fits-all solution. You can divide up the classes in many different ways and still get usable results and the best way to divide them may be highly problem-dependent. For instance: in the case of a discretized continuum variable, you'd probably want to make sure that all the borders only separate "adjacent" classes, that is, ranges that are next to each other in the discretization scheme. On the other hand, if there are four classes, but they can be reduced into two overlapping sets of two classes, then that right there tells you how to divide them up.There are two major ways you can divide up the classes: in the first case so that all of the conditional probabilities can be solved for and in the second case so that only the conditional probability of the "winning" class is returned. In the first case, you need at least (n-1)/2 borders--that is, for each border you can write two equations relating the conditional probabilities plus one more for the entire system (they sum to unity.) In the second case, the method is hierarchical: divide the set of classes into two groups, divide each pair of groups into two more and so on. Consider the case of dividing up four classes:
x | o
-----
+ | *
The four classes are x, o, + and *. In the first case, we think of this diagram as representing two borders: one horizontal, and one vertical. In the second case, consider first a horizontal border, and then two vertical borders. It stands to reason that the two methods can also be combined.
References
[1] Peter Mills (2011) "Efficient statistical classification of satellite measurements." International Journal of Remote Sensing, 32 (21): 6109-6132.[2] Donald Knuth (1998) The Art of Computer Programming, Volume 3: Sorting and Searching, 2nd edition Addison-Wesley.
[3] William Press, Saul Teukolsky, William Vetterling, Brian Flannery (1992). Numerical Recipes in C, 2nd edition Cambridge University Press, Cambridge.
[4] Peter Mills (2012). "Principal component proxy tracer analysis." arxiv:1202.2194
[5] Peter Mills (2006), "The influence of noise on a classical chaotic scatterer." Communications in Nonlinear Science and Numerical Simulation, 11 (8) 899-906.
Thursday, January 31, 2013
Bias in drug trials
I just came across this talk and I think it is very important. It deals with the tendency for journals and researchers to publish only positive results while either ignoring or rejecting negative results. The problem exists not just in clinical medicine, but in every area of science. Science is a very competitive business and positive results often herald a "new discovery." By contrast, we tend to think that a negative result says very little, thus it's in scientists' own best interests to publish only positive results, assuming that the journal will even accept the negative ones.
From my own experience, I have seen scientists not only picking the best exemplar from a series of experiments and publishing only that one, but also essentially cherry-picking the data by using some creative averaging. I also know that I am not immune from this myself.
From my own experience, I have seen scientists not only picking the best exemplar from a series of experiments and publishing only that one, but also essentially cherry-picking the data by using some creative averaging. I also know that I am not immune from this myself.
Subscribe to:
Posts (Atom)