A moment with the Universe.


A sad grey twilight watches me from outside the window. Faint traces of the sun’s orange being given up to the quiet waiting blue some way up from the horizon.

How long has it been since I opened these little attics of my mind? How long since I last deemed they were of value? How long since I let the universe touch me, filling a moment? How long since I let it hold me and stare into my eyes, unhurried with questions?

Aren’t there really two ways of seeing everything? One says that the only things that matter are those detached of irrational feeling, while the other says that only how you feel about the universe is real. In the end, we are left to choose. How long has it been since I last chose the second?

Darkness now trickles closer to the horizon, chasing the last orange around the planet, leaving the city draped in beads of harsh sodium.

Isn’t it sad how it gets harder every day to see the face of the universe and contemplate its incredible existence, even for a fleeting moment? Isn’t it ironic that what is most difficult to see is the greatest encompassing truth? Is it not childish to clutter one’s fleeting time with small playthings, that have only to do with marginally extending that time? Isn’t it sad how many hurdles must be overcome, in violation of the code of compliance, to sit down with the universe for a moment of uninterrupted privacy?

Isn’t it sad how, despite knowing all this, one must chase this conditioned lunacy day after day until time runs out?

How will I be any different? It is an increasingly escaping possibility that I shall.

For now, just a moment will have to do.

Until the next time, universe. Hope to see you again. I hope you can find it in your heart to forgive me.

Thoughts from a drunken party

Are there any utilitarian side-effects to a conscious attempt towards the dissolution of the self? It might seem an odd question to ask, for the very task of this self-learned dissolution is only taken up by those who have reached a state of mind that has conclusively disengaged from any utilitarian goals.

I suppose I need to elaborate my definition of utilitarian in this context. The dissolution of the self and ego are lofty tasks that are hard and long to achieve. Are there other observable effects that arrive earlier and easier on the way that lend towards the final goal and towards happiness and peace in general? That is what I meant by the question.

I realized tonight on the way home from a drunken and disappointing party that there is at least one; I am sure there are more.

I realized that we are too eager to judge people we do not know too well yet, but have the slightest bit of information that we can base any judgement on (a bit of this was going on in both directions at the party). The poverty of information about people we have just met enables us not to hit any contradictions as we construct an image of them that lends heavy credence towards this judgement that we wish to paste on them.

Why are we so eager to judge people in the first place? The answer is obvious: it makes us feel better about our own selves. This is the underlying incentive and motivation behind judgement of any kind really.

Now comes in the utilitarian part of working to dissolve the ego. If you have taken the goal seriously to any degree and have given it some thought, you realize that working to elevate yourself above the perceived moral or intellectual characters of others is damaging to this goal, for it only serves to build the ego that you have realized is doing you no good. Not only is it damaging to this particular goal, it is a misdirected ambition to pursue, for it can only be fueled by things such as information-deprived quick judgements of people. The more you get to know and like the other person, the more they are on your side and they become your friends, an extension of your self and the things you stand for, and then you start using only their good qualities to judge other people and their friends against. It is an entirely opportunistic and self-deceptive process once you start to think about it.

Therefore, a willingness to erode the significance and importance we attach to the ego brings with it a necessary inhibition from being quick judges of other people. One is more accepting, forgiving and patient with people as they grow to understand them as a human being. In fact, I believe it makes it easier to imagine oneself to be the other person, which is a thought that if attended to with any sincerity even for a few seconds, makes harsh judgements nearly impossible henceforth.

And of course, such a quality would bring happiness, companionship, trust and peace in many ways, long before the dissolution of the self is achieved. I think there are many fruits worthy of picking on the way to this final goal; in fact the journey might scarcely be worth it without these collateral pickings.

Calculating the Lyapunov Exponent of a Time Series (with python code)

(In a later post I discuss a cleaner way to calculate the Lyapunov exponent for maps and particularly the logistic map, along with Mathematica code.)
I found this method during my Masters while recreating the results of an interesting paper on how some standard tests for chaos fail to distinguish chaos from stochasticity (Stochastic neural network model for spontaneous bursting in hippocampal slices, Biswal and Dasgupta, 2002).
The Lyapunov exponent is a measure of sensitive dependence on initial conditions, i.e. how quickly two nearby states diverge.
Now consider two points in the time-series, ti and tj, whose values are very close. That means the system reached near the same state at the ith and jth iterations. Now consider the two sequences  ti,  ti+1,  ti+2 … and  tj,  tj+1,  tj+2 … We wish to know how these two sequences diverge from each other. For this, consider the distance between the two sequences after k steps: d(k) = | ti+ktj+k |. (This is for a 1D time series. For higher dimensions, you can define this to be the Euclidean distance and modify the code accordingly.) If the system is chaotic, d(k) will initially rise exponentially with k. For this, one can plot ln d(k) vs k and apply a linear fit. The slope will be an estimate for the Lyapunov exponent.
(Since the system is bounded, the two nearby states will not diverge indefinitely though. Their exponential divergence will stop after some length. We must fit the straight line only within this region.)
Now, this was for a single pair of initial states. The Lyapunov exponent is an average of this divergence exponent over all nearby initial pairs. So for this, define d(k)>, where is averaging over all starting pairs  ti,  tj, such that the initial distance d(0) = | t– tj | is less than some fixed small value. The program finds all such initial pairs, calculates d(k)>, plots it against k, and the slope of the initial linear part gives us the Lyapunov exponent.
Python Code
The following code takes a text file with the time series, ‘timeseries.txt’, as the argument. The text file must contain only the time series values in a single column, no serial numbers or any other text before or after. It asks for the starting diameter within which to limit the initial pairs. It displays how many such pairs it is finding in the time series, so you can vary the diameter based on this.
It outputs a text file, ‘lyapunov.txt’ with two columns, k and d(k)>, which you can then plot and fit in the correct region by visual inspection.
from math import log

def d(series,i,j):
    return abs(series[i]-series[j])

f=open('timeseries.txt', 'r')
series=[float(i) for i in f.read().split()]
eps=input('Initial diameter bound: ')
dlist=[[] for i in range(N)]
n=0 #number of nearby pairs found
for i in range(N):
    for j in range(i+1,N):
        if d(series,i,j) < eps:
            print n
            for k in range(min(N-i,N-j)):
for i in range(len(dlist)):
    if len(dlist[i]):
        print>>f, i, sum(dlist[i])/len(dlist[i])

The following is the plot and fit of the resulting data from a logistic map series with an appropriately chosen initial diameter.

Lyapunov Exponent of Logistic Map

I deliberately did not automate the plotting and fitting part, because a. it’s tedious and hard to write the code in a way that runs on most installations, and b. human eyes will do a much more reliable job of identifying where the linear portion ends.

R code for multivariate random-walk Metropolis sampling

I couldn’t find a simple R code for random-walk Metropolis sampling (the symmetric proposal version of Metropolis Hastings sampling) from a multivariate target distribution in arbitrary dimensions, so I wrote one. This is also my first R code.
It requires the package MASS to sample from the multivariate normal proposal distribution using the mvrnorm function. If you are using R on Windows, you can download the package zip for Windows from the link, and use Packages > Install package(s) from local zip files… from the GUI to install the package.
The reason I couldn’t write the code for a general Metropolis algorithm (i.e. for any arbitrary symmetric proposal distribution, not just normal) or a more general Metropolis-Hastings algorithm (with any arbitrary proposal distribution, symmetric or not) is that generating the proposal point would then require sampling from an arbitrary proposal distribution. This is only easy for a few standard distributions, but hard in general (which is the point of using such algorithms in the first place).

I. Function

The following is the function that does the Random Walk Metropolis-Hastings sampling when supplied with the required arguments. Notes about the arguments follow the code.

rwmetro <- function(target,N,x,VCOV,burnin=0)
    require(MASS)   #requires package MASS for normal sampling
    samples <- x
    for (i in 2:(burnin+N))
        prop <- mvrnorm(n = 1, x, VCOV)
        if (runif(1) < min(1, target(prop)/target(x))) 
            x <- prop
        samples <- rbind(samples,x)

II. Arguments

  1. target function
    The function defining the multivariate target distribution, written as a function of an n-vector, where n is the number of variables on which the distribution is defined. The different variables of your distribution must be written as x[1], x[2] etc.
    An example is the following, defining the function f(x,y)=exp(-5*abs(x^2+y^2-1)) in two dimensions:

    ring2D <- function(x)	# x is a vector
  2. N integer
    The final sample size (i.e., excluding the burn-in length).
  3. x numeric vector
    The starting point (vector) of your Metropolis-Hastings chain. This vector needs to be the same length as the dimension of the target distribution. Your target function must be able to accept this vector as an argument.
  4. VCOV numeric matrix
    The variance-covariance matrix for the multivariate normal that is used as the proposal distribution for random-walk Metropolis-Hastings. The length of this matrix also has to be the same as the dimension of the target distribution, i.e. the length of vectors acceptable by target and the length of x. You can vary the entries of this matrix and observe your results to see what works better for sampling your target distribution.
    The following line defines a matrix in two dimensions with .01 variance for each variable and no covariance between them.

    vcov2D <- .01*diag(2)
  5. burnin (optional) integer
    The ‘burn-in’ length for the chain. The number specified will be the number of initial samples chucked. If nothing is specified, it’s taken to be 0.

III. Output
numeric matrix
The output is a matrix where each row is a sample from your target distribution, excluding the initial burn-in samples. The number of rows is thus the sample size, and the number of columns is equal to the dimension of the target distribution. You can use this matrix however you want, to save, visualize or analyze your sample.
An example output of sampling a 2D distribution with sample size 5. Each row is an x,y sample.

[1,] 0.12961923 0.03708061
[2,] 0.10765053 -0.02798036
[3,] 0.01112930 -0.07255766
[4,] 0.06049625 -0.04546265
[5,] 0.1832164 -0.1244354

IV. Usage Example
Let’s take the target distribution we used as an example above, f(x,y)=exp(-5*abs(x^2+y^2-1)). This looks like a ring of radius 1 rising from the x-y plane and centered at the origin. Here is a gnuplot surface plot of the distribution (because I found it frustratingly hard to figure out a level plot in R):

Target Distribution in Gnuplot
Let’s generate a sample of size 40,000 from this distribution with the starting point (0,0) and without any burn-in length, and with the variance-covariance matrix we defined before. This is done by calling the function with the correct arguments:

ringsample2D <- rwmetro(ring2D,40000,c(0,0), vcov2D)

This assumes that you have already defined the target ring2D and the matrix vcov2D in the way explained, preceding this function call.
ringsample now contains a random sample from the distribution.
With the following one-line R code, I made a scatter plot of the sample from the ringsample matrix.

plot(ringsample2D[,1], ringsample2D[,2], xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main='Metropolis-Hastings Sample',xlab='x', ylab='y', pch='.')

The following is the plot:
Metropolis-Hastings Sample, scatter plot

Putting all the code together in sequence, this is what the full code for defining the arguments, drawing the sample, and making the plot for this example looks like:

# Define arguments
ring2D <- function(x)    # x is a vector

vcov2D <- .01*diag(2)

#Define the sampling function
rwmetro <- function(target,N,x,VCOV,burnin=0)
    require(MASS)   #requires package MASS for normal sampling
    samples <- x
    for (i in 2:(burnin+N))
        prop <- mvrnorm(n = 1, x, VCOV)
        if (runif(1) < min(1, target(prop)/target(x))) 
            x <- prop
        samples <- rbind(samples,x)

# Call the sampling function with the arguments
ringsample2D <- rwmetro(ring2D,40000,c(0,0),vcov2D)

# Use the sample
plot(ringsample2D[,1], ringsample2D[,2], xlim=c(-1.5,1.5),ylim=c(-1.5,1.5), main='Metropolis-Hastings Sample',xlab='x', ylab='y', pch='.')

The following is another sample, this time for a 3-variable function, f(x,y,z)=exp(-5*abs(x^2+y^2+z^2-1)). This is just the 3D version of the previous function, and its density is mostly concentrated on a spherical shell of radius 1 and of some thickness, centered at (0,0,0). Notice how we never have to spare a thought about normalizing these functions to get the sample, which is one of the advantages of the Metropolis-Hastings algorithm.
The target function is now:

ring3D <- function(x)	# x is a vector

and the variance-covariance matrix is, say, .01 down the diagonal again:

vcov3D <- .01*diag(3)

Let’s go for a sample of 20000 points after a burn-in of 20000, starting from (0,0,0). The function is then called as:

ringsample3D <- rwmetro(ring3D,20000,c(0,0,0),vcov3D,20000)

ringsample3D is now a 3×1000 matrix, and I use the following to get a rotatable 3D plot of it (requires the rgl package):

plot3d(ringsample3D[,1], ringsample3D[,2], ringsample3D[,3], xlab='x',ylab='y', zlab='z',col='red', size=1)

Metropolis Sample

I would have preferred a spinning 3D gif, but I’m too tired to work that out right now.
Share the knowledge.

Partners meet halfway: a simple correlation study of an undergrad lab class

Last semester I taught two classes of an Introductory Physics (Electrodynamics) for Engineers lab course at the University of Texas at Austin. The first day that my students came in, most of them did not know each other. They just came in and sat at tables of two, and this pretty much became their permanent seats in the room for the rest of the course (that’s human habit). They were supposed to work with their partner on the experiments in the first part of the class, but individually on their worksheets in the second part. So choosing their seats also fixed the partners. Thus, it is a good assumption to say that the pairing had been pretty random, at least as far as academic ability of the students was concerned.

I noticed though that both partners in a group tended to move towards a common average in their performance over the course. If one partner was diligent, it would motivate (and sometimes force) the other to work harder, while the sloppiness and uncaring attitude of a partner would also negatively affect the performance of the other. I was happy to notice though that the ‘uplifting’ effect was usually more dominant than the ‘down-dragging’. In any case, this caused the performance of partners to pull closer, irrespective of their starting academic ability.

At the end of the semester, after I had the final scores that I handed out to them, I decided to do a simple correlation test. The following are two scatter plots of pairs of student scores for each of my classes, each pair ordered as (better score, worse score) (which is why all points are below the bottom-left to top-right diagonal of the plot). The linear fits through them have also been shown. Although I had suspected that there would be a significant correlation between partners’ scores, I was surprised at how high the correlation coefficients turned out to be. (Note that if you pair an even number of numbers randomly, on the average you should expect a linear correlation coefficient of zero.)





I had felt this trend even by the middle of the semester, and I had this crazy idea of using this to try and pull up the performance of the students who were not doing too good. I thought of sorting the students by their score so far into the course, and re-assigning partners so as to pair the best student with the worst, the second-best to the second-worst and so on. However, I know that habit is a deadly force, and if I suggest my students in the middle of the course that they have to say goodbye to their partner with whom they already have some steady chemistry going and have new partners again, ones that I have chosen for them, all hell would break loose. So I had to rest it in my head as a theoretical result with possible practical applications, to be taken up by a braver soul than me. But hey, here’s the data, so you can quote it if you want to try it and someone doesn’t take the idea too well.