Using Python to Model a Levy Path

I have been on a bit of a project to broaden my mind and revist things which I had studied years ago using Coursera. As part of this, I took a course on Statistical Mechanics that was offered by Professor Werner Krauth at the Ecole Normale et Superieure in Paris. (Taught in english fortunately for me, as I wouldn't have even thought of trying to deal with it in french - my high school french is completely lost to me). Even in this case I was struggling a bit, and only got to the middle of the course before I decided that I was too far behind to catch up. I saved all the homework exercises and lecture videos, and resolved that I would come back and work through them as soon as I had the chance. The course involves describing physical systems by modeling them as systems of statistical behavior, rather than the deterministic models one might think of constructing using freshman physics. Strangely, the results produced by these models are identical in many circumstances to those produced by deterministic models. The course was taught using python as an environment to construct the models, which was a good choice for producing simple code which one could run and modify on the spot, but it put me in a little bit of pain moving back and forth between R, octave, (an opensource MATLAB clone), and python with the various courses I was taking. I found myself continually wondering, now what was the function called that did that again? And what were the arguments?...evidently not those...,hm.  I very much wanted to settle on one particular language and not jump back and forth, though I still haven't figured out which one of them that would be. The course began with an introduction to Markov chain methods and simulations, something which I had vaguely heard of, but never really known what people were talking about, and then moved into modelling simplified physical systems using Markov methods. I got to the first section of Quantum Statistical Mechanics and then decided that I was not going to be able to keep up with the material, so now I have been working through it more gradually. The section I am working on now involves Levy paths. These are something else I vaguely remember hearing about without having any idea what they were. They evidently appear in models of phenomena of all sorts, from quantum mechanics to fiancial analysis.


 This course described them as being a simpler way of obtain essentially the same output as the Feynman path integral, and, at least in the demo program, it really isn't different from a random walk with each change of direction being chosen from a normal distribution. In fact, a levy path can be constructed this way, though one must correct back to the final x end position and move each x slice back proportionally as well.



# Add a random normal with mean of x-previous.
Upsilon.append(random.gauss(Upsilon[-1], sigma))


 # Subtract the distance between upsilon_final and xend
 x = [0.0] + [Upsilon[k] + (xend - Upsilon[-1]) * \ k / float(N) for k in range(1, N + 1)]



One can see that the last position selected becomes the mean for the distribution from which the next position is selected. In the class we were provided with a Markov implementation of the path integral first before learning about Levy paths. The implementation involved generating proposed x positions from a random uniform distribution, calculating new weight and old weight, and then throwing away everything that did not fit into a random normal distribution. This, of course, makes no particular sense as an algorithm, but it is exactly how a Markov method works. One does not know anything about the probabilities of the solution space. One simply walks randomly around the solution space, and calculates the values of randomly chosen points. The wrapping algorithm that I was supposed to implement for this exercise caused me some confusion as well. One was supposed to implement some sort of list wrap like:


L = L[3:] + L[:3]


The first time I worked through the exercise, I did not understand that the previously generated value was being used to generate the current value, so things were just being scattered on a normal distribution, rather than gradually tending toward the value of xend. I also did not see how large versus small beta affected the calculation until I did a graph of the internal variables of levy_harmonic_path(). I still am not sure I am gaining anything by wrapping the list. I had thought it would be necessary for the sort of Markov method implementation that would be needed to simulate a Bose-Einstein condensation. One would just have the program grind round a million times and it would never be able to leave the list, but I am not really sure things work that way with this particular function - when you call it, it generates a new list every time, and overhead seems quite high for no particular gain. On the other hand, maybe I am just not quite understanding what they have in mind. Now, let's have a look at the inside of the function that is supposed to construct a Levy path and see if any sense can be made of it. The function is called like:

levy_harmonic_path(xstart, xend, dtau, N)

 xstart and xend are the initial and final positions of the particle of course. I found myself considerably confused by the lecture's insistence that the particle begin and end at the same position. The code given did not have any such requirement, but it is probably needed for a Markov simulation which always seems to run inside a boundary-less box with particles leaving from one side reappearing on the other. Dtau is beta divided by N, so it represents the temperature loss with each slice of the simulation. N is the number of slices we have divided the simulation into. In order to calculate a possible x position at a given slice we need to have the previous x position and the final x position. These will be used to calculate the mean and standard deviation of a normal distribution from which a random selection will be made for the new x position on the current slice. We are only moving in one dimension here, but we have a second dimension because we are moving from high temperature, or low beta, to low temperature, or high beta. The position of the particle becomes more and more uncertain as beta gets larger, but it has a definite position at xstart and xend, at least as far as this code is written. (This might be another reason why the code must be altered to wrap, particularly if one were to observe a Bose-Einstein condensation). This particular function models the behavior of a particle in a potential well, which means that a force is acting on it to keep it near x = 0. Now, when considered statistically, a force means that one is more likely to be found in one region than another, rather than some sort of causal factor producing a particular trajectory. This may seem obvious, but I had to think a bit before it occured to me.

Ups1 = 1.0 / np.tanh(dtau) + 1.0 / np.tanh(dtau_prime)

Ups2 = x[k - 1] / np.sinh(dtau) + xend / np.sinh(dtau_prime)

mean = Ups2/Ups1

sd = 1.0/np.sqrt(Ups1)

This is how the function calculates the mean and standard deviation of the distribution from which to make a random selection for the x position of the current slice.

With dtau_prime = (N-k) * dtau, so this is the total tau for all the slices remaining until we reach the final value of beta. (Each slice changes by dtau). Now, the first thing which occurred to me was that I had no real idea what the functions in the denominator looked like, so we can have a look at hyperbolic sine and hyperbolic tangent:

Hyperbolic sine seems to increase exponentially, and so flattens out everything the larger dtau and dtau_prime are, but hyperbolic tangent approaches an asymptote of 1, so as dtau becomes larger 1/tanh(dtau) will tend toward 1. dtau depends entirely on how large beta is, while dtau_prime decreases with each slice, so the formula has the effect of pulling us on a track toward xend, at least with small beta. When beta is large, Ups1 will tend toward 1 and Ups2 will tend toward 0, so we will make a random selection of x from the standard normal distribution, mean zero, standard deviation 1.

 BETA1.0 With small beta, the course is determined. One starts at the initial position and ends at the final position with only small deviations, but with large beta things are different.

BETA20.0 Now things look different. The mean for selection is being calculated by Ups1/Ups2 with this generally tending to 0/0. The position is very uncertain and the initial and final positions have little effect of the position of the particle at any given slice.

Trying out some techniques from Coursera Machine Learning Class

I finished the Machine Learning course from Andrew Ng on Coursera recently. The last topic he covered in the course was on recommender systems. There was a technique he described that involved multiplying a matrix of user preferences by a matrix of movie scores to fill in data on movies which had not been rated, and to allow the system to make recommendations for movies the users had not seen.  It immediately occured to me that this should be a general method one could use to fill in missing data in many sorts of situations.  I decided that I would need to investigate the technique further after the class ended.  I would try to find some relatively simple data file and see how one would go about applying the technique to it.  During a lecture, professor Andrew Ng mentioned that the technique he was describing was known as Low Rank Matrix Factorization.  I did a bit of web searching and found more information under "Low Rank Matrix Reconstruction."

When I had some free time, I went straight to the UC Irvine Machine Learning Lab website to see if I could find some relatively simple, real world data to try the technique on.  I found a data file on cardiac arrhythmia which seemed about right for a test: not too big, a fair amount of missing data, mostly numeric, and possibly low rank.  I know nothing about electro-cardiogram procedures, but thanks to Wikipedia, I discovered that they involved recording data through twelve leads.  The data file itself has 279 attributes recorded for 452 individuals.  In principle, all these 279 attributes should be created by some mixture of the twelve variables of the different leads, so this data should definitely be a low rank matrix, and possible to reconstruct using one method or another.

First I needed to do some data munging on the file. Beginning with code I found at:

http://www.togaware.com/datamining/survivor/Cardiac_Arrhythmia.html


I wrote some code in R to label all the variables, remove non-numeric columns, split the sample into training and test sets, and save the results into Rdata files.

Setup.R

First, I tried running svd() on the dataset with all NA's removed. This reduced the dataset to only 68 complete rows, but it seemed likely to give me some idea of how many features there would be in an optimal basis.

Eyeballing things, it seems as though 10 to 20 features would be a reasonable reduced basis, at least from the complete data examples.  This is good, since it does not obviously conflict with the number of leads used in the apparatus, which is the reason to settle on this number of parameters. It seems that we should be able to apply the method of minimizing parameters to this data to recover the NA values.

So on to the code which should let us try this in R. First, a file of which specifies the cost and gradient functions for optim(), then a script to run them and see whether this approach recovers missing values which we have zeroed out ourselves for a test.

CollaborativeFiltering.R


Test.R

The functions in CollaborativeFiltering.R are a basic translation of what I wrote in the Coursera Machine Learning class from MATLAB, or octave, to R.

Regrettably, none of this worked quite as I had thought. First, I had a bit of a struggle with optim() just understanding how one should use it with two entire matrices as a parameter. I decided that I should just flatten everything to a vector, just as I had in octave. Then, much more confusing, I found that optim() was zeroing out the entire second matrix. Since the process requires that one multiply matrix X * Theta' to regenerate the data matrix, this would not work at all. I found that when I fed the data in without any sort of normalizing, (which was required at least by the course lectures), I got back two matrices, but when I touch the data in any way, (either by dividing by the max of each feature or subtracting the average), the second matrix came back filled with zeros. This was quite useless for a method that required multiplying one matrix by another.

Since I wasn't having success with R, I decided to return to the original octave code, and see whether the method worked on this data. Unfortunately, it did not work there any better. The code ran and produced a matrix of the proper dimensions of cardiac readings by patients, but the numbers were not even the same order of magnitude at many points in the regenerated matrix.

Some things which occur to me:
The movie score matrix was very simple. It was composed of integers between 1 and 5 with many missing data points. Next, the problem that was being solved wasn't actually restoring data, it was yet another classification problem. There are two classes which a movie can belong to in this problem, these being interesting, and uninteresting. We don't need to accurately construct an unknown score, we only need to sort this particular movie into one class or the other. I hadn't been aware that I was actually solving a classification problem while I was working through the code for class. In fact, classification problems seem to be the only sort of problems Machine Learning concerns itself with at least by the course material for this class. This again is something that I had not noticed until I started working on this data set. Possibly this method could be improved by adding another variable like the lambda to push down on the mean squared error section of the code, but it would probably be best to try another approach. I read of a method that involved reconstructing a matrix by minimizing the sum of the singular values returned by svd(). This seems more likely. I might try this, as well as trying a neural net and K-means on this data set if I have some more free time.

A Jewel at the Heart of Quantum Physics - Some Thoughts

Browsing through my daily dose of slashdot I found a link I just had to follow, A Jewel at the Heart of Quantum Physics. This article described a mathematical structure that enabled one to solve calculations of particle interactions that had formerly been completely beyond calculation, in some cases on a sheet of paper. In order to arrive at this simplification the theory requires abandoning two core ideas:

Locality is the notion that particles can interact only from adjoining positions in space and time. And unitarity holds that the probabilities of all possible outcomes of a quantum mechanical interaction must add up to one. The concepts are the central pillars of quantum field theory in its original form, but in certain situations involving gravity, both break down, suggesting neither is a fundamental aspect of nature.
In keeping with this idea, the new geometric approach to particle interactions removes locality and unitarity from its starting assumptions. The amplituhedron is not built out of space-time and probabilities; these properties merely arise as consequences of the jewel’s geometry. The usual picture of space and time, and particles moving around in them, is a construct.

After reading the article, a number of questions occur to me which I do not have the mathematical chops to even begin to answer.

1. Is this the only structure possible, or are there  a large, or even infinite number of similar structures, each of which could form the basis of a potential universe?
2. Is this structure hidden in the quantum wave equation when you look at it  in the right way? It seems it really should be.
3. I am imagining something like a standing wave pattern that might have formed by addition with each phase transition of the Big Bang. Does this work at all?
4. Is this structure stable and eternal essentially, or is it temporary?

Money - value

Money is a strange thing, completely abstract and mathematical in certain regard, yet held to be the most tangible of all things. Lately, people have been making the rounds on NPR, etc., expressing the thought that our ability to rapidly and accurately assess value is increasing, and will increase greatly in the future. This will make it possible to end "second chances," and rapidly sort the productive from the unproductive. Value and money are actually rather trickier concepts than this, and their relationship rather more complex. Financial panics would not occur if our ability to assess value were truly as advanced as such claims would have us believe, since a financial panic is a market perception that monetary valuation and emotional and tangible valuation are seriously out of alignment. Money is a very strange instrument. Consider the case of Long Term Capital. Here, a hedge fund constructed a mathematically formula which enabled them to sell the risk out of their portfolio, while keeping the profit, thus apparently being able to invest with no risk at all. No real asset that one could use money to purchase behaves in such a manner. A farm, or factory, or piece of equipment all have risk of loss or damage attached inescapably to any investment. I would argue that a market panic is an essential process which forces market valuations back into conformity with emotional and tangible valuations.

Human "errors" with risk estimation and chance

Economists and statisticians are always telling us that humans fundamentally misjudge risk, and do not follow optimal strategies in all manner of situations.  I would say that we more likely following strategies that are optimized for evolutionary survival in the sort of environments which animals must move around in.  People who haven't taken a course in probability theory seem to assume that every example of chance pulls from sort of pool - tossing a coin uses up some pool of possible "heads" results, so that after a run, there must not be any "heads" left in the pool, and we must start getting "tails."  The sort of probability which people intuitively model functions much more like a deck of cards than a coin toss or a roll of the dice.  One might wonder why this might be, imagining oneself as some sort of animal, a small monkey for instance, since this would be where we might make evolutionary contact with our methods of estimating risk.  A monkey is in a tree, looking for ripe fruit.  When it finds a piece that is very nice, it has drawn a card from the ripe fruit deck.  Pretty soon, there are no cards left in the ripe fruit deck, and only unripe fruit remains.  The monkey could climb to the next tree, but it knows there might be pythons there.  There are no pythons in this tree, while in the next tree, there is a deck of python cards, with no clear knowledge of what it contains.  If you draw a python card, then more than likely, you are dead.  If you escape, then you have turned over the python card, and can avoid it.  The pool of fatal risk has been reduced by some amount, at least in the near term.  Everything the monkey is doing stays in fixed pools, in some limited interval at least, with the risk of error being quite possibly fatal.  This is the sort of risk evolution has designed us to process and estimate I would say.  Pools are closed, with each draw reducing the pool.  Cost of risk is very high.  Eventually, one must seek new resources, but risk should always be minimized.  The monkey is playing a card game rather than rolling dice.  But this game with predators will look different than the game with one's own species,

Putting Blackbox Models in Charge

Something I read in the Economist this week struck me, (April 20, 2013 edition - pg. 30 "That swooning feeling").  This seems worth quoting completely.  It is speaking about various theories about the repeated apparent dip in the economy each spring.  "One theory is that models interpreted the economy's plunge in late and 2008 and early 2009 as partly seasonal, and responded by nudging up subsequent winter figures and nudging down summer data to compensate."  Really!  Are the models that the Bureau of Labor Statistics this autonomous and opaque that no one who is in charge of using them can even figure out whether such a thing is occurring?  Their testing solution for this evidently works like this - "But the federal Bureau of Labor Statistics has found that the pattern persists even if the job numbers are seasonally adjusted without those recession months."  Could they really have so little comprehension of the model they are using, that they can't figure out whether such an effect is likely, and the only way of testing is to remove the data and run it again?  The stock market swings up and down, hiring and lay-off decisions are made, investments made or abandon on this number, and evidently no one has any idea how it actually comes into existence.  It would be better to publish the raw employment surveys and let people construct their own analysis than to use a black box of this sort.

Inflation and Wealth Concentration vs. Asset Bubbles

The Fed has been pumping money out into the economy for some years now, leading to a continual gnawing fear in some circles that inflation is about to skyrocket.  I would say that this is overlooking a transformation which has occurred over the course of the last couple of decades.  First, I would make some hypotheses about money.  I would say that money wants to multiply, and that money wants to concentrate.  These are a bit anthropomorphic, but money is concretized human desire, so a bit of anthropomophic speculation is perfectly reasonable.  Next, money dispersed over a broad population is what sustains demand, and therefore makes it possible for prices to rise - inflation in other words.  When money is concentrated beyond a certain degree, inflation will no longer be possible, instead asset bubbles will predominate.  Money will always attempt to concentrate and multiply, so left to its natural tendencies, it will flow to a smaller and small portion of the population.  This small proportion of the population cannot sustain demand of a broad consumer market, so inflation cannot increase.  Formerly, it was necessary to hire people and pay them an increasing wage in order to make more money, but automation and financial instruments have decoupled this connection to a great degree.  Paying money to labor is dispersion, and the return is most likely logarithmic, while investing in obscure financial instruments is concentration, and the return is exponential.  It is easy to see what the trend will be.  Now, the question would be, how to put some numbers on this and make some sort of convincing case?