Sampling a distribution

I was listening to an introductory talk about MCMC (link). I'll have more to say about the talk another time (hint: I think it's terrific). The speaker described a method for sampling from a probability distribution. The idea is to generate random samples from U[0,1], and then find the value of x at which the cumulative density function first exceeds the random number. The resulting samples of x are plotted as a histogram.

The Python script has three parts. In the first section we set up a weighted mixture of three normal distributions of different mean and sd, and plot that in red using x increments of 0.01. In the second section the same values are used to generate a discrete cdf for the same points. This is plotted in blue (after normalizing by the interval size).

Then we implement the algorithm for sampling. Binning for the histogram used the collections.Counter class that will be in Python 2.7. I don't have that, so I grabbed the class from here.

You might notice that the labels on the y-axis also cover the histogram. It would be better to have two separate plots, but I'm not quite sure how to handle that yet.

[UPDATE: more here]

import math, sys import numpy as np import matplotlib.pyplot as plt import Counter def normal(mu,sigma): def f(x): z = 1.0*(x-mu)/sigma e = math.e**(-0.5*z**2) C = math.sqrt(2*math.pi)*sigma return 1.0*e/C return f p1 = normal(0,2) p2 = normal(10,1) p3 = normal(18,0.5) # sum of weighted normal distributions @np.vectorize def p(x): return 0.5*p1(x) + 0.25*p2(x) + 0.25*p3(x) dx = 0.01 xmax = 25 R = np.arange(-10,xmax+dx,dx) # dashed lines plt.plot((R[0],R[-1]),(0,0),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(1,1),color='k',ls=':',lw=2) plt.plot((R[0],R[-1]),(-0.5,-0.5),color='k',ls=':',lw=2) L = p(R) plt.plot(R,L,color='r',lw=3) #=========================================== cdf = [L[0]] for e in L[1:]: cdf.append(cdf[-1] + e) cdf = np.array(cdf) cdf *= dx plt.plot(R,cdf,color='b',lw=3) ax = plt.axes() ax.set_xlim(-6,xmax) ax.set_ylim(-0.55,1.05) #=========================================== def find_first(n,L): for i,e in enumerate(L): if n < e: return i return len(L) samples = list() width = 0.4 f = 1/width for i in range(1000): n = np.random.random() # must adjust to actual range value = find_first(n,cdf)*dx - 10.0 # trick to truncate at fractional values samples.append(np.floor(f*value)/f) #samples = np.array(samples) c = Counter.Counter(samples) maxn = c.most_common(1)[0][1] for k in c: n = c[k] n = n * 0.45 / maxn r = plt.Rectangle((k,-0.5), width=width,height=n, facecolor='green') ax.add_patch(r) plt.savefig('example.png')

Peaks of density in numpy

Here is a Python script that generates peaks of user-defined density. It has two main parts. Given an existing matrix and a point x,y (in the interval 0-1), it constructs a matrix of the same shape containing the distance of each element from the designated point.

The second part uses the distance information to construct peaks of density to add to the matrix. The score function is vectorized. For each input element of a distance matrix, it computes a score based on (i) a threshold beyond which the score is zero, (ii) a height for the maximum score, and (iii) a function describing how the density changes with distance.

Nothing fancy, and a bit slow, but it does make it easy to play with the parameters and generate a density like that in the sample. At the end, we need to make sure each value is <= 1.0 and >= 0, because make_heatmap converts these directly to "gray" values to input to the Rectangle function.

Here is output for a 10 x 10 example:

0.04 0.09 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.09 0.81 0.09 0.04 0.04 0.05 0.16 0.22 0.16 0.05 0.04 0.09 0.04 0.04 0.04 0.16 0.41 0.56 0.41 0.16 0.04 0.04 0.04 0.04 0.04 0.22 0.56 1.0 0.56 0.22 0.04 0.04 0.04 0.04 0.04 0.16 0.41 0.56 0.41 0.16 0.04 0.04 0.04 0.04 0.04 0.05 0.16 0.22 0.16 0.05 0.04 0.07 0.09 0.07 0.04 0.04 0.04 0.04 0.04 0.04 0.07 0.17 0.28 0.17 0.07 0.04 0.04 0.04 0.04 0.04 0.09 0.28 0.81 0.28 0.09 0.04 0.04 0.04 0.04 0.04 0.07 0.17 0.28 0.17 0.07 0.04 0.04 0.04 0.04 0.04

import matplotlib.pyplot as plt import numpy as np import HM_Utils_frac def init_data(N=30): A = np.zeros(N**2) A.shape = (N,N) return A def eucl_distance(x0,y0,x,y): return np.sqrt((x-x0)**2+(y-y0)**2) def distance_matrix(A,x,y): # matrix is 1.0 x 1.0 in "size" # given a matrix A with shape R,C # each row,col in M is a fractional pos R,C = A.shape L = [1.0*c/C for c in range(C)] M = list() for r in range(R): for c in range(C): M.append((L[c],L[r])) # now compute each dist from given point x,y M = [eucl_distance(x0,y0,x,y) for x0,y0 in M] M = np.array(M) M.shape = A.shape return M @np.vectorize def score(d,ht,thr,func): # given a distance from x,y compute score # first, check threshold if d > thr: return 0.0 d = thr - d # scale d to func(d) and height return ht* (func(d) / func(thr)) def make_peak(A,x,y,ht,thr,func=None): M = distance_matrix(A,x,y) if not func: func = lambda x: x**1.5 return score(M,ht,thr,func) def scale(A): maxN = np.max(A) if maxN > 1.0: A = A/(1.0*maxN) return A def f(x): return x**4 N = 40 ax,fig = HM_Utils_frac.init_fig() A = init_data(N) A += make_peak(A,x=0.7,y=0.3,ht=1.0,thr=0.3) A += make_peak(A,x=0.2,y=0.8,ht=0.8,thr=0.4,func=f) A += make_peak(A,x=0.1,y=0.1,ht=0.8,thr=0.2,func=f) A = A + 0.04 A = scale(A) if N <= 10: HM_Utils_frac.show_array(A) HM_Utils_frac.make_heatmap(A,ax) plt.savefig('example.png')

Yet another heatmapper

I'm planning to try a small Python MCMC simulation, utilizing a simple 2D grid of densities. I wrote a new heatmap script as a first step. I'll post the code today and then use it in later updates.

One small thing, the original version asked matplotlib for the xlim and ylim and then calculated fractional square dimensions to use. Since it seemed a bit slow, I changed the script so that we tell the Axes object what dimensions we want (the array's shape, naturally). Then there are no fractions and the increment is +1. However, it was not significantly faster at all. Here is the timeit output:

$ python -m timeit -s 'from HM_Utils import demo' 'demo(10)' 10 loops, best of 3: 225 msec per loop $ python -m timeit -s 'from HM_Utils_frac import demo' 'demo(10)' 10 loops, best of 3: 229 msec per loop

Here's the listing:

import matplotlib.pyplot as plt import numpy as np def init_fig(): fig = plt.figure() ax = fig.add_subplot(111) ax.xaxis.set_visible(False) ax.yaxis.set_visible(False) return ax,fig def show_array(A): for row in A: L = [str(round(n,2)).ljust(4) for n in row] print ' '.join(L) # origin is lower left def make_heatmap(A,ax): R,C = A.shape ax.set_xlim(C) ax.set_ylim(R) for r in range(R): for c in range(C): n = A[r][c] assert 0.0 <= n <= 1.0 p = plt.Rectangle((r,c), width=1,height=1, # black is 0.0, w is 1.0 facecolor=str(1-n), edgecolor='w') ax.add_patch(p) def demo(N=10): ax,fig = init_fig() A = np.random.uniform(size=N**2) A.shape = (N,N) if N <= 10 and __name__ == '__main__': show_array(A) make_heatmap(A,ax) plt.savefig('example.png') if __name__ == '__main__': demo(50)

Selecting array elements in R and numpy

I was a bit confused about the behavior of vectorize in numpy (last post), so I posted a question to SO. Of course, they knew the answer immediately, as I described.

But there is one more thing that came up, and that has to do with the selection of elements from an array. I wasn't aware that you can do indexing by a test that returns a boolean vector in numpy, but you can. In R, it's common to see something like this. Here is a matrix m. We can leave a row out by "-":

> m = 1:9 > dim(m) = c(3,3) > m = t(m) > m [,1] [,2] [,3] [1,] 1 2 3 [1,] 4 5 6 [2,] 7 8 9 > m[-1,] [,1] [,2] [,3] [1,] 4 5 6 [2,] 7 8 9

The test m[,1] > 2 asks for all rows in which the first column is greater than 2:

> m[,1] [1] 1 4 7 > sel = m[,1] > 2 > sel [1] FALSE TRUE TRUE

Note: I screwed this up in the original example. We use the selector to get those rows:

> m[sel,] [,1] [,2] [,3] [1,] 4 5 6 [2,] 7 8 9

How would you do this in Python with numpy?

import numpy as np A = np.arange(1,10) A.shape = (3,3) print A

[[1 2 3] [4 5 6] [7 8 9]]

B = A[A[:,0]>2] print B

[[4 5 6] [7 8 9]]

sel = A[:,0]>2 print sel

[False True True]

print A[sel,:]

[[4 5 6] [7 8 9]]

It works!

Vectorize in numpy

I'm working on a little project using numpy and matplotlib where I found a need to write my own function that operates on numpy arrays. So, for example, we could have, as a first attempt:

import numpy as np A = np.arange(0,1,0.1) def f(n): if n > 0.3: return n return 0 print A print f(A)

$ python vec.py [ 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9] Traceback (most recent call last): File "vec.py", line 7, in print f(A) File "vec.py", line 4, in f if n > 0.3: return n ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Numpy is trying to tell us that we can't just feed an array to a function that operates on single values. The function has to be "vectorized." Say:

f = np.vectorize(f) print f(A)

[0 0 0 0 0 0 0 0 0 0]

We didn't get a ValueError, but the result is not correct. What happened? According to the help(np.vectorize):

Data-type of output of vectorized is determined by calling the function with the first element of the input.

The problem is that although the array contains floats, we specified 0 as the value to return for n < 0.3. The first n evaluated was 0.0 which is less than 0.3 and so returns the integer 0. That leads to trouble, since the vectorized function now returns all the rest of the values as ints. Or so the docs imply. (Also see here). There are a couple more issues. First:

A[0] = 2.0 print f(A) B = np.arange(1,0,-0.1) print f(B)

[2 0 0 0 0 0 0 0 0 0] [1 0 0 0 0 0 0 0 0 0]

The output above is obtained if it is run in the same script as the first part. That is, the data type to be returned is only evaluated once, when the function is first called. If the first call is commented out, then we get this:

[ 2. 0. 0. 0.3 0.4 0.5 0.6 0.7 0.8 0.9] [ 1. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0. 0. ]

One way to fix the whole problem is to return a float. The following example also uses vectorize as a decorator

@np.vectorize def f2(n): if n > 0.3: return n return 0.0 print f2(A)

[ 0. 0. 0. 0.3 0.4 0.5 0.6 0.7 0.8 0.9]

Continuing in the docs:

This can be avoided by specifying the otypes argument as either a string of typecode characters or a list of data-types specifiers.

We're supposed to do something like:

def f3(n,otypes=[np.float]): def f3(n,otypes='d'):

But, in my tests, neither of these can compensate for returning the integer 0. The real reason for the otypes argument seems to be that there is a large number of numpy dtypes.

Unifrac (8): understanding PCoA---calculate

I'm exploring Unifrac (first and second posts).

In the previous post, I showed code and some output related to our topic for this post. You'll need to generate those numbers to follow along. Of course, you could just paste 'em in:

[[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

The whole trick to PCoA is that it can take a "measure of association" between two samples, like a distance. So, let's use the values from before, calculate the distances between them, and then see how well PCoA recovers the original data. The second half of the code listing is at the bottom of the post. In the first segment of this second half, we simply define two functions to calculate the distances between all the points, using either the "manhattan" distance or the "euclidean" distance. We will only use the euclidean version for this analysis. If you want to see the results of the distance calculation, you'll have to run the code (manhattan) or modify it to print the euclidean result.

Next, we want to do the PCoA analysis. I did this two ways (to check myself). First, I looked under the hood at the PyCogent module cogent.cluster.metric_scaling, and I call the same functions used in metric_scaling.principal_coordinates_analysis. These include the following two functions as well as a numpy.linalg routine:

E = metric_scaling.make_E_matrix(eucl) F = metric_scaling.make_F_matrix(E) eigvals, eigvecs = eigh(F)

followed by a little massaging of the vectors to generate the output. As you will see if you run the code, the eigenvalues and then the last two rows of my result are:

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.3 36.34 ... [ 1.9 -1. -2. -1.9 1.2 -0.8 1.3 3.4 2.4 -4.5] [-18. -14. -10. -6. -2. 2. 6. 9.9 14. 18.1]

The other approach is to call PyCogent's function directly. The eigenvalues and the last two results of the call to

pca = metric_scaling.principal_coordinates_analysis point_matrix, eigvals = pca(eucl)

are

-0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 53.28 1320.32 ... [ 1.9 -1. -2. -1.9 1.2 -0.8 1.3 3.4 2.4 -4.5] [-18. -14. -10. -6. -2. 2. 6. 9.9 14. 18.1]

This looks pretty good, though I notice that my eigenvalues aren't matching. I'm missing a square root somewhere. I will postpone worrying about this for another day.

>>> 1320.32/53.28 24.780780780780781 >>> (36.34/7.3)**2 24.781302308125362

The original matrix is:

[[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

I think this is pretty remarkable. We've taken a matrix of distances between our sample points, done a little PCoA mumbo-jumbo, and regenerated something pretty close to the original.

However, I am a bit worried at what look to be bugs in the implementation:

• the original x-values are systematically transformed by + 2
• the original y-values are systematically transformed by * -1

This gives a shift and a reflection of the points, but the relationships revealed by the PCoA are identical to the original.

[UPDATE: <whacks self on head> Duh...I bet that in view of the last line above, if we calculate the distance we'll get the same result as we did for the original matrix. So, of course, there's no reason that PCoA could possibly tell the difference. Presumably, it's just centered its result a little better around (0,0). ]

#============================================ def delta(n1,n2): return abs(n2-n1) def manhattan_dist(x1,y1,x2,y2): return delta(x1,x2) + delta(y1,y2) def euclidean_dist(x1,y1,x2,y2): sumsq = delta(x1,x2)**2 + delta(y1,y2)**2 return np.sqrt(sumsq) eucl = list() manh = list() for i in range(len(x)): for j in range(len(x)): x1,y1 = A[0][i], A[1][i] x2,y2 = A[0][j], A[1][j] manh.append(manhattan_dist(x1,y1,x2,y2)) eucl.append(euclidean_dist(x1,y1,x2,y2)) manh = np.array(manh) manh.shape = (N,N) eucl = np.array(eucl) eucl.shape = (N,N) print 'manhattan\n', manh, '\n' #============================================ # PyCogent stuff starts here: # I call their functions E = metric_scaling.make_E_matrix(eucl) F = metric_scaling.make_F_matrix(E) eigvals, eigvecs = eigh(F) t_eigvecs = eigvecs.transpose() print 'eigvals' for e in eigvals: print round(e,2), print print 't_eigvecs' for row in t_eigvecs: print np.round(row,decimals=1) eigval_root = np.sqrt(abs(eigvals)) print 'eigval_root\n' for e in eigval_root: print round(e,2), print final = t_eigvecs * eigval_root[:,np.newaxis] print 'my calculation' for row in final: print np.round(row,decimals=1) print # They call their functions print 'PyCogent' pca = metric_scaling.principal_coordinates_analysis point_matrix, eigvals = pca(eucl) print 'eigvals' for e in eigvals: print round(e,2), print print 'point_matrix' for row in point_matrix: print np.round(row,decimals=1) print print 'original matrix (A)\n', A

Unifrac (7): understanding PCoA---setup

I'm exploring Unifrac (first and second posts).

One of the nice aspects of UniFrac is the "multivariate analysis" it can do including PCoA (Principal Coordinates Analysis). I'm trying to understand how this works, and I have two posts to go through it. This first post is similar to stuff we've done before (e.g. here), but I want to try to present a coherent picture. Sorry for the repetition. The code is presented in two parts overall, and then today's section has two parts also. We begin by getting a bunch of integers in x-coordinates over the range from -20 to +20, and in y-coordinates randomly between -5 and +5. We turn those into a 2D array and plot them using matplotlib as red circles.

A [[-20 -16 -12 -8 -4 0 4 8 12 16] [ -2 1 2 2 -1 1 -1 -3 -2 5]]

Next, the array is rotated through 45 degrees and replotted (blue circles).

In the second part (of the first half), we use standard Principal Component Analysis (PCA) on these data points to recover their long dimension (now rotated) back to the x-axis. We calculate the covariance matrix for these two sets of points, and then do np.linalg.eigh(m) to get the eigenvalues and eigenvectors.

cov [[ 74.08888889 70.35555556] [ 70.35555556 78.53333333]] eigenvalues [ 5.9204692 146.70175303] eigenvectors [[-0.71818168 0.69585564] [ 0.69585564 0.71818168]]

The eigenvalues come out with the largest one last. We grab the first eigenvector (the last in the array we received), and do

np.dot(ev1,B)

which gives back the x-values we started with.

[-20. -16. -12. -8. -4. 0. 4. 8. 12. 16.]

We can do the same with the second eigenvector to get the y-values. We can combine these two operations by re-ordering the eigenvectors:

print 'C' T2 = np.array([ev1,ev2]) C = np.dot(T2,B) for row in C: print np.round(row) print

C [-20. -16. -12. -8. -4. 0. 4. 8. 12. 16.] [-2. 1. 2. 2. -1. 1. -1. -3. -2. 5.]

In the figure, the C array has been plotted as black cross-bars, which match the original red circles pretty well. Next time, we'll take the same data points and look at them from the perspective of PCoA. We'll also see how PyCogent fares with our data, and how reasonable a picture it can estimate of what we started with.

Here's the code:

import numpy as np from pylab import * from cogent.cluster import metric_scaling plotting = len(sys.argv) > 1 np.random.seed(1357) x = range(-20,20,4) y = np.random.random_integers(-5,5,size=len(x)) A = np.array([x,y]) N = len(x) A.shape=(2,N) print 'A\n', A n = 1.0/np.sqrt(2) T = np.array([[n,-n],[n,n]]) B = np.dot(T,A) if plotting: scatter(A[0],A[1],s=250,color='r') scatter(B[0],B[1],s=250,color='b') #============================================ m = np.cov(B) print 'cov\n', m e = np.linalg.eigh(m) print 'eigenvalues\n', e[0] print 'eigenvectors\n', e[1] ev1 = e[1][1] print np.round(np.dot(ev1,B)) ev2 = e[1][0] print np.round(np.dot(ev2,B)) print print 'C' T2 = np.array([ev1,ev2]) C = np.dot(T2,B) for row in C: print np.round(row) print if plotting: scatter(C[0],C[1],s=100,c='k',marker='+') R = (-25,25) xlim(R) ylim(R) savefig('plot.png')

Confidence intervals in matplotlib

This plot comes from the same place as the previous one (Cumming 2007 PMID 17420288 ).

It shows the influence of sample size on the standard error and confidence interval, and compares these values to the standard deviation. Standard error (se) is calculated as sd / √n, and the confidence interval is calculated as the t statistic * se, where t is hard coded in the script. The value approaches 1.96 (for p-value of 0.05) as n gets large. I have to investigate more how to determine this with numpy or PyCogent.

import numpy as np from pylab import * data1 = [53,38,29] data2 = [16,29,32,34,40,44,49,51,52,53] data3 = np.random.normal(loc=39,scale=10,size=30) data = [data1,np.abs(data2),data3] def do_error_bar(x,e,lw=2,w=2): o = plot([x,x],[m+e,m-e],color='k',lw=lw) o = plot([x-w,x+w],[m+e,m+e],color='k',lw=lw) o = plot([x-w,x+w],[m-e,m-e],color='k',lw=lw) #layout from 1 to 100 N = len(data) total = 100 margin = 5 space = 15 fig_width = total - margin*2 all_spaces = (N-1)*space total_group_width = fig_width - all_spaces group_width = total_group_width * 1.0 / N dx = group_width/3.0 # 4 elements per group x = margin # start tD = {3:3.182, 10:2.228, 30:2.042} S = 50 # large circle size s = 10 # small circle size y = 17 # y pos for text for group in data: for g in group: scatter(x,g,s=s,color='k') x += dx n = len(group) m = np.mean(group) sd = np.std(group) scatter(x,m,s=S,color='k') do_error_bar(x,sd) text(x-1.5,y,'std') text(x,y-5,'n = ' + str(n)) x += dx se = sd/np.sqrt(n) t = tD[n] ci = t * se scatter(x,m,s=S,color='k') do_error_bar(x,ci) text(x-1.5,y,'CI') x += dx scatter(x,m,s=S,color='k') do_error_bar(x,se) text(x-1.5,y,'se') x += space ax = axes() ax.xaxis.set_visible(False) show()

Bar plot with data points

I saw this here (Cumming 2007 PMID 17420288 ).

There are still a few issues to work out. Most important, I don't know how to specify the limits for the axes. This one looks OK, but it doesn't always come out so, and I haven't been able to find the right call to fix them.

[UPDATE: I forgot all about it! I need to surpress the x-axis labels!]

Another issue is that the error bar thickness doesn't seem to be settable. There's no API (linewidth refers to the bar rather than the error bars), and the objects that come back don't seem to have an appropriate reference. One last detail is that the code for the distance between a set of data points and the bar is not rational. But I think it's not too bad for a novice.

Here is the listing:

import numpy as np from pylab import * data1 = np.linspace(12,22,7) data2 = np.random.normal(loc=15,scale=5,size=10) data3 = [11,12,18,19,20,26,27] data = [data1,np.abs(data2),data3] # n = number of groups def layout(n,r=7): s = r**2 # r = radius of each data point #layout from 1 to 100 margin = 5 spacer = 10 group_width = (100 - 2*margin - (n-1)*spacer)*1.0/n dot_width = r bar_width = group_width - dot_width current = margin rL = list() for i in range(n): rL.append(current) # x for point rL.append(current + 3) # x for bar current += group_width + spacer return s, bar_width, rL s, w, xlocs = layout(len(data)) for group in data: x = xlocs.pop(0) for e in group: scatter(x,e,s=s,color='k') m = np.mean(group) e = np.std(group) x = xlocs.pop(0) o = bar(x,m,width=w,color='0.6', yerr=e, ecolor='k') show()

[UPDATE: The answer on StackOverflow was to draw the error bars manually. Call this from the main loop:

def do_error_bar(x,e,lw=3): w = 3 o = plot([x,x],[m+e,m-e],color='k',lw=lw) o = plot([x-w,x+w],[m+e,m+e],color='k',lw=lw) o = plot([x-w,x+w],[m-e,m-e],color='k',lw=lw)

And do this before calling show():
axes().xaxis.set_visible(False)

The new figure:

Comparing R and Python sequences

This is a post about elementary sequence operations in R and Python. It's as much for me as for you.

The most obvious difference between sequences in R and Python is Python's use of 0-based indexing:

R:

> 1:5 [1] 1 2 3 4 5 > A = seq(1,10,by=2) > A [1] 1 3 5 7 9 > A[2] [1] 3

Python:

>>> range(1,6) [1, 2, 3, 4, 5] >>> A = range(1,11,2) >>> A [1, 3, 5, 7, 9] >>> A[1] 3

Another difference is that in R, but not in Python, one can assign to an index outside the initial range:

> m = 1:2 > m[6] = 35 > m [1] 1 2 NA NA NA 35

>>> m = range(1,3) >>> m[6] = 35 Traceback (most recent call last): File "", line 1, in IndexError: list assignment index out of range

In R, but not in regular Python, we can make the increments non-integral values:

> A = seq(0,20,by=0.1) > A[1] [1] 0 > length(A) [1] 201 > A[length(A)] [1] 20

We can use numpy to get around this restriction:

>>> import numpy as np >>> A = np.arange(0,20.1,0.1) >>> A[0] 0.0 >>> len(A) 201 >>> A[-1] 20.0

It's sometimes more convenient to specify how many numbers we want to obtain (evenly spaced in some interval):

> A = seq(0,2,length=6) > A [1] 0.0 0.4 0.8 1.2 1.6 2.0

>>> A = np.linspace(0,2,6) >>> A array([ 0. , 0.4, 0.8, 1.2, 1.6, 2. ])

Vectorized operations:

> m = 1:9 > dim(m) = c(3,3) > m [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > m = t(m) > m [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 > apply(m,1,mean) [1] 2 5 8 > apply(m,2,mean) [1] 4 5 6 > mean(m) [1] 5

>>> m = np.arange(1,10) >>> m.shape = (3,3) >>> m array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> np.mean(m, axis=0) array([ 4., 5., 6.]) >>> np.mean(m, axis=1) array([ 2., 5., 8.]) >>> np.mean(m) 5.0

Here are some examples of fancy indexing where we rearrange rows and columns both at the same time:

> m [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 > m[c(2,3,1),c(3,2,1)] [,1] [,2] [,3] [1,] 6 5 4 [2,] 9 8 7 [3,] 3 2 1

The naive implementation in Python gives something different than in R (though useful), but the sequential approach works:

>>> m array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> m[[1,2,0],[2,1,0]] array([6, 8, 1]) >>> m[[1,2,0],:][:,[2,1,0]] array([[6, 5, 4], [9, 8, 7], [3, 2, 1]])

R has a few more indexing tricks for which I don't know if there is a Python equivalent:

> m = 1:9 > dim(m) = c(3,3) > m = t(m) > m[-1,] [,1] [,2] [,3] [1,] 4 5 6 [2,] 7 8 9 > sel = m[1,] > 2 > sel [1] FALSE FALSE TRUE > m[sel] [1] 7 8 9 > y = -5:5 > y [1] -5 -4 -3 -2 -1 0 1 2 3 4 5 > y[y < 0] <- -y[y < 0] > y [1] 5 4 3 2 1 0 1 2 3 4 5 > y <- abs(y) > y [1] 5 4 3 2 1 0 1 2 3 4 5

But here, finally is an example we can do in both:

> m [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 > sel = array(c(1:3,3:1), dim=c(3,2)) > sel [,1] [,2] [1,] 1 3 [2,] 2 2 [3,] 3 1 > m[sel] = 0 > m [,1] [,2] [,3] [1,] 1 2 0 [2,] 4 0 6 [3,] 0 8 9

>>> m array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> m[[0,1,2],[2,1,0]] = 0 >>> m array([[1, 2, 0], [4, 0, 6], [0, 8, 9]])

Multiple bars in a bar-graph in matplotlib

I answered a question on Stack Overflow. It was posted 3 hr ago, and the guy's never coming back, so I'll put it up here and you can enjoy it. If you want the axes all purty and everything, check out the manual.

import numpy as np from pylab import * data = [[6.7, 1.5, 4.5], [2.0, 3.25, 5.7]] w = 0.5 xlocations = np.array(range(len(data)))+w colors = ['0.2','0.5','0.8'] oL = list() for x,d in zip(xlocations, data): for c,value in zip(colors, sorted(d,reverse=True)): b = bar(x, value, width=w, color=c) oL.extend(b) show()

Two sample t-test in Python

[UPDATE: I'm currently (November. 2010) in the middle of a series of posts about Student's t test. The first six posts are here, here, here, here, here, and here.]

I'd like to do my teaching in Python as much as possible. In keeping with that, I plan to explore methods for analyzing microarray data using Python rather than R and Bioconductor. As a preliminary matter, I'm going to use Student's t-test (two sample) as implemented in PyCogent. A simple-minded approach would be to a situation where the samples have already been classified, would be to filter them by the t-test statistic.

Here is a simple test. We draw a bunch of numbers from the normal distribution and keep them in a numpy array. We slide a window along the array, pulling out two small sets of n numbers at a time, and do a two sample t-test on them. We record the p-value. We expect that since all the numbers come from the same normal distribution the distribution of the p-value will be such that, on the average, p < 0.05 about 5% of the time.

# two sample t-tests import numpy as np import cogent.maths.stats.test as stats def oneTrial(n,f=np.random.normal): N = 500 # num of individual tests SZ = n*2*N # need this many nums draw = f(loc=50,scale=3,size=SZ) counter = 0 for i in range(0,SZ,n*2): nums1 = draw[i:i+n] nums2 = draw[i+n:i+2*n] t, prob = stats.t_two_sample(nums1,nums2) if prob < 0.05: counter += 1 return 1.0*counter / N

Now, we run the trial a bunch of times. For this run, I've chosen the sample size to be 3.

n = 3 # sample size L = list() for i in range(100): #if i and not i % 10: print 'i =', i L.append(oneTrial(n)) print 'avg %3.4f, std %3.4f' % (np.mean(L), np.std(L)), print 'min %3.4f, max %3.4f' % (min(L), max(L))

Output:

avg 0.0511, std 0.0095 min 0.0280, max 0.0780

Checkerboard in Numpy

Here's a question on SO about how to make a checkerboard in Numpy (like you see on pdfs when as "the classic representation for "no pixels", or transparent." I have no doubt that I had the correct answer, but I could not buy a vote. Even from a guy with 22.9 K rep. Sheesh.

We use hstack and vstack to stack "unit" squares of solid color:

>>> import numpy as np >>> u = 15 >>> b = np.zeros(u**2) >>> b.shape = (u,u) >>> w = b + 0x99 >>> >>> width = 20 # squares across of a single type >>> row1 = np.hstack([w,b]*width) >>> row2 = np.hstack([b,w]*width) >>> board = np.vstack([row1,row2]*width) >>> board.shape (600, 600)

Now just slice it down to what's needed

>>> width_actual = 357 >>> height_actual = 239 >>> board = board[:height_actual,:width_actual] >>> board.shape (239, 357)

Broadcasting in Numpy

The issue of array shapes matters when trying to combine them by arithmetic operations. Here is a simple array of 4 rows and 3 columns:

>>> from numpy import ones >>> A = ones((4,3)) # 4 rows x 3 cols >>> A.shape (4, 3) >>> A array([[ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.], [ 1., 1., 1.]])

Here is a second array of 4 rows x 1 col

>>> ones((4,1)) # 4 rows x 1 col array([[ 1.], [ 1.], [ 1.], [ 1.]])

We can add them, and in the process the second array is "broadcast" to the same shape as the first one:

>>> A + ones((4,1)) array([[ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.]])

An array of a single row and 3 cols works the same:

>>> ones((1,3)) # 1 row x 3 cols array([[ 1., 1., 1.]]) >>> A + ones((1,3)) array([[ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.]])

Now, it is possible to make a 1D array, e.g.:

>>> B = ones((3,)) # a 1D array >>> B array([ 1., 1., 1.]) >>> B.shape (3,)

And again, we can add A + B:

>>> A + B array([[ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.]])

But this doesn't work if the number of columns doesn't match:

>>> C = ones((4,)) # a 1D array >>> C.shape (4,) >>> C array([ 1., 1., 1., 1.]) >>> A + C Traceback (most recent call last): File "", line 1, in ValueError: shape mismatch: objects cannot be broadcast to a single shape

We need to call newaxis:

from numpy import newaxis >>> D = C[:,newaxis] >>> D.shape (4, 1) >>> A + D array([[ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.]])

Matplotlib in OS X (4)

While hunting around for colors for heatmaps, I came across the matplotlib class matplotlib.colors.LinearSegmentedColormap. The way the colors are specified seemed a bit strange at first, but then it became clear.

The first thing is that, naturally, the red, green and blue components are handled separately. They may vary independently and don't even have to be split up into the same number of "intervals." Color values are specified corresponding to the numerical endpoints of each interval, and when it is desired to encode a number as a color, this is done by linear interpolation. Intervals and colors are specified by values between 0 and 1.

An example I found here splits the interval 0 to 1 into two parts for all three colors. The red values come from this array:

r = ((0.0, 0.0, 0.0), (0.5, 1.0, 0.7), (1.0, 1.0, 1.0))

What is going on is that we have split the interval into two parts: 0 to 0.5 and 0.5 to 1.0. These are the entries [0][0], [1][0] and [2][0]. The manual calls these x. From there:

row i: x y0 y1 / / row i+1: x y0 y1 / / row i+2: x y0 y1

Suppose we encounter a number that lies between the x in row i and the x in row i+1. We grab the indicated values (y1 from row i, and y0 from row i+1), and just interpolate.

So we have 7 values of interest for this set of two intervals:

• the beginning and end of each interval
(where the beginning of the second equals the end of the first)

• red at the beginning and end of interval 1

• red at the beginning and end of interval 2

Two values are never used---y0 from the first row and y1 from the last. Check out the second link for lots more examples.

import matplotlib.pyplot as plt import matplotlib import numpy as np r = ((0.0, 0.0, 0.0), (0.5, 1.0, 0.7), (1.0, 1.0, 1.0)) g = ((0.0, 0.0, 0.0), (0.5, 1.0, 0.0), (1.0, 1.0, 1.0)) b = ((0.0, 0.0, 0.0), (0.5, 1.0, 0.0), (1.0, 0.5, 1.0)) colors = {'red':r, 'green':g, 'blue':b} f = matplotlib.colors.LinearSegmentedColormap m = f('my_color_map', colors, 256) plt.pcolor(np.random.rand(10,10),cmap=m) plt.colorbar() plt.savefig('example.png')

Matplotlib in OS X (3)

I've been playing with matplotlib, and I have to say that the figures it produces are beautiful. Although the documentation is extensive, I do find it hard to navigate to the right pages, or find helpful examples by direct search. But of course, GIYF.

Here is a proof of principle for a heat map in a few short lines of code. I still have to work on placing text, getting the right colors, etc. The fact that it's all Python is terrific. We use a scatterplot and adjust the size of the plot symbol (a square) to be as large as we need.

Probably it's obvious, but the call plt.axis('equal') is crucial because by default the axes are not typically equal and the "squares" are therefore not square. The call yL = sum(yL,[]) is a clever trick I found on Stack Overflow for flattening a list.

[UPDATE: As is not uncommon, I found a simpler API after making this post. The result is shown below. The code is the second of the two samples.]

import matplotlib.pyplot as plt import random xL = range(1,11) * 10 yL = [[j]*10 for j in range(1,11)] yL = sum(yL,[]) R = xrange(100000) cL = [random.choice(R) for i in range(len(xL))] fig = plt.figure() ax = fig.add_subplot(111,axisbg='0.9') # light gray plt.axis('equal') ax.scatter(xL, yL,c=cL, marker='s',edgecolor='w',s=1250) plt.savefig('example.pdf')

import matplotlib.pyplot as plt import numpy as np import random R = xrange(10000) data = [random.choice(R) for i in range(100)] data = np.array(data) data.shape = (10,10) fig = plt.figure() plt.hot() plt.pcolormesh(data) plt.colorbar() plt.savefig('example.png')

Matplotlib in OS X (2)

Now that I have a working intall of matplotlib, I'm going to explore its capabilities to replace R, at least for plotting stuff. The more heavy duty components of R will still be irreplaceable.

The matplotlib example for histograms is here. It's not very intuitive, at least at first sight. I took nearly everything in the example code and wrapped it in a function with a silly name and placed it in a module called H. If you want to follow this example it should be simple to figure out what I did.

[UPDATE: Although this works well, and is said to be more efficient, there is a simpler API. See example here]

The script that connects my data with the plot code is just a couple of lines:

import random from random import gauss import matplotlib.pyplot as plt from H import doWeirdStuffToMakeAHistogram random.seed(133) mu = 10 sigma = 3 L = [gauss(mu,sigma) for i in range(1000)] doWeirdStuffToMakeAHistogram( L,bins=40,x0=0,x1=20,col='blue') plt.savefig('example.pdf')

Matplotlib in OS X (1)

Although matplotlib is a sophisticated approach to plotting and it uses Python, I have used R instead for two reasons: Bioconductor uses R, and I couldn't get matplotlib installed properly. Since PyCogent uses matplotlib and I've got it working now, I thought I'd take a look.

The immediate goal is to get some numbers to use to work up a simple example of Principal Coordinates Analysis with PyCogent (and try to remember the difference with Principal Component Analysis). According to the abstract of this book chapter:

The problem is that PCA is based on the correlation or covariance coefficient, and this may not always be the most appropriate measure of association. Principal coordinate analysis (PCoA) is a method that, just like PCA, is based on an eigenvalue equation, but it can use any measure of association (Chapter 10). Just like PCA, the axes are plotted against each other in a Euclidean space, but the PCoA does not produce a biplot (a joint plot of the variables and observations).

This should all be similar to what we did previously in R. Here is a Python script, in two parts. The first part generates generates the points plotted in red as separate lists of x and y values, and makes a scatterplot. 's' is for the size of the plot character (weirdly, 'size' doesn't work).

The second part uses the array manipulation from numpy to rotate these points. The transformation matrix is t, and the multiplication is achieved by dot(t,a). The matplotlib OS X GUI doesn't work properly, so I just save the output as a PDF. (In contrast, the R OS X GUI is beautiful and always works perfectly for me. I think the matplotlib developers must be Windows or Unix guys. That would also explain why the Makefile is out of date).

The figure is a little off because I haven't figured out how to set the plot limits yet.

Here's the code:

import matplotlib.pyplot as plt from numpy import * from math import sqrt import random random.seed(1357) R = range(-100,100) def fx(): return random.choice(R) r = range(-10,10) def fy(): return random.choice(r) xL = [fx() for i in range(30)] yL = [fy() for i in range(30)] fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(xL,yL,s=250,color='r',marker='o') #---------------------------------------------- a = xL+yL a = array(xL + yL) a.shape = (2,len(xL)) z = 1.0/sqrt(2) t = [[z,-z],[z,z]] t = array(t) xL,yL = dot(t,a) ax2 = fig.add_subplot(111) ax2.scatter(xL,yL,s=250,color='b',marker='o') plt.grid(True) plt.savefig('example.pdf')