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.