Newton's approximation for pi

Borrowing shamelessly from Journey Through Genius, we'll finish the example of Newton's version of the binomial expansion from here.

Start with a semicircle as shown in the diagram below. The equation of the circle is

(x - 1/2)2 + (y - 0)2 = (1/2)2 x2 - x + 1/4 + y2 = 1/4

Solve for y:

y = √(x - x2) = √x √(1-x)

The positive square root is the equation of the semicircle, where √(1-x) can be replaced by its binomial expansion:

1 - (1/2)x - (1/8)x2 - (1/16)x3 - (5/128)x4 ..

and multiplied by x^1/2 to give:

x1/2 - (1/2)x3/2 - (1/8)x5/2 - (1/16)x7/2 - (5/128)x9/2 ..

Now let B be the point (1/4,0) as shown.

The area of the arc of the circle ACD in the diagram comprises two areas: the triangle BCD plus the area ABD. Newton used fluxions to calculate the area of ABD. In modern terminology, we integrate the equation of the semicircle given by the binomial expansion:

(2/3)x3/2 - (1/2)(2/5)x5/2 - (1/8)(2/7)x7/2 - (1/16)(2/9)x9/2 - (5/128)(2/11)x11/2 ..

evaluated for between x = 0 and x = 1/4. Now the reason for picking these particular coordinates is revealed:

(1/4)3/2 = (√(1/4))3 = 1/8 (1/4)5/2 = (√(1/4))5 = 1/32 ..

So the series above reduces to:

1/12 - 1/160 - 1/3584 - 1/36846 .. = 0.07677

(more precision shown in Dunham's book, achieved easily by Newton). Using the Pythagorean theorem instead we get that

BD = √((1/2)2 - (1/4)2) = √(3/16) = √3/4

And the area of the triangle BCD is

(1/2) √3/4 (1/4) = √3/32

Since the base of BCD is 1/4 and the radius is 1/2, the cosine of the angle ACD (BCD) is 0.5 and ACD equals 60 degrees, meaning that the total area of ACD is 1/6 of the circle:

π/24

Using the value obtained from the series above:

π = 24(0.07677 + √3/32) = 3.1415..

A contemporary recorded this account of Newton's Lucasian lectures: "so few went to hear him, & fewer yet understood him, that oftimes he did in a manner, for want of Hearers, read to ye Walls." The commentator added that Newton's lectures would last for half an hour except when there was no one at all in the audience, in that case he would stay only 15 minutes.

EMBOSS: first explorations

This is a quick update on my exploration of EMBOSS. It's not a set of detailed notes to follow or anything like that. I'm just trying to get a feeling for how to use it, and recording my impressions.

As mentioned last time, I can launch a window with the JEMBOSS GUI like so:

$ /usr/local/share/EMBOSS/jemboss/runJemboss.sh

On Mac OS X, it installed into /usr/local/share... I found a tutorial here:

/usr/local/share/EMBOSS/doc/tutorials/emboss_tut.tar.gz

Unpack and grab emboss_tutorial.pdf. Nice short tutorial. We get info about a program like say wossname by one of these:

wossname -opt wossname -help wossname -outfile ~/Desktop/file.txt wossname > help.txt

Even I can figure out that what wossname does is to search for info on a specific topic, say "restriction":

$ wossname Finds programs by keywords in their short description Text to search for, or blank to list all programs: restriction SEARCH FOR 'RESTRICTION' rebaseextract Process the REBASE database for use by restriction enzyme applications recoder Find restriction sites to remove (mutate) with no translation change redata Retrieve information from REBASE restriction enzyme database remap Display restriction enzyme binding sites in a nucleotide sequence restover Find restriction enzymes producing a specific overhang restrict Report restriction enzyme cleavage sites in a nucleotide sequence showseq Displays sequences with features in pretty format silent Find restriction sites to insert (mutate) with no translation change

I put a FASTA-formatted sequence file on my desktop and do:

$ remap example.fasta Display restriction enzyme binding sites in a nucleotide sequence Comma separated enzyme list [all]: Minimum recognition site length [4]: Output file [da247.remap]: EMBOSS An error in remap.c at line 244: Cannot locate enzyme file. Run REBASEEXTRACT

$ REBASEEXTRACT Process the REBASE database for use by restriction enzyme applications REBASE database withrefm file: Error: Input file is required

Google led me to REBASE. I was confused because the example input files shown on a help page here don't match the three files I downloaded. But then I noticed that (some of) the example output files match. REBASE/embossre.enz looks like link_emboss_e.txt.

Actually reading the docs: the input file must be the "withrefm" file of a REBASE distribution. For example, the withrefm file for REBASE version 005 is at: ftp://ftp.neb.com/pub/rebase/withrefm.005. But the link is dead.

Google: site:rebase.neb.com withrefm

No matches. So, I can't find the required input files on the REBASE site.

[UPDATE: the files are there, see the end of the post]

I found what look like equivalents for two of the four output files including REBASE/embossre.enz. The reference file format looks like it's changed. And the fourth file REBASE/embossre.equ isn't there.

Check out the files in

$ ls /usr/local/share/EMBOSS/data/REBASE/ dummyfile embossre.enz embossre.equ embossre.ref embossre.sup

Some look good, like embossre.equ, but some are dummy files, e.g. embossre.enz.

Try this:

sudo mv embossre.enz embossre.enz.old sudo cp ~/Desktop/link_emboss_e.txt /usr/local/share/EMBOSS/data/REBASE/embossre.enz $ remap -sequence example.fasta Display restriction enzyme binding sites in a nucleotide sequence Comma separated enzyme list [all]: Minimum recognition site length [4]: 6 Output file [da247.remap]: results.txt

The result looks like things are working:

DA247 NheI | Cac8I Cac8I | | BmtI | StuI NspI XmnI \ \ \ \ \ \ \ GATGAACGCTAGCGGCAGGCCTAACACATGCAAGTCGAGGGAGAAGCCCTTCGGGGCGGA 10 20 30 40 50 60 ----:----|----:----|----:----|----:----|----:----|----:----| CTACTTGCGATCGCCGTCCGGATTGTGTACGTTCAGCTCCCTCTTCGGGAAGCCCCGCCT / / / / / / / | | NheI | StuI NspI XmnI | Cac8I Cac8I BmtI D E R * R Q A * H M Q V E G E A L R G G M N A S G R P N T C K S R E K P F G A E * T L A A G L T H A S R G R S P S G R K ----:----|----:----|----:----|----:----|----:----|----:----| S S R * R C A * C M C T S P S A R R P P X H V S A A P R V C A L R P L L G E P R I F A L P L G L V H L D L S F G K P A S

So, what I might do is look harder for withrefm, but it seems to work, at least. Not crazy about the output format, they should do something like:

. . . . . 60 GATGAACGCTAGCGGCAGGCCTAACACATGCAAGTCGAGGGAGAAGCCCTTCGGGGCGGA CTACTTGCGATCGCCGTCCGGATTGTGTACGTTCAGCTCCCTCTTCGGGAAGCCCCGCCT

[UPDATE: Should've looked harder. The files I need are right there on the REBASE downloads page (#5 and #31).

$ sudo rebaseextract Password: Process the REBASE database for use by restriction enzyme applications REBASE database withrefm file: ~/Desktop/link_withrefm.txt REBASE database proto file: ~/Desktop/link_proto.txt c-98-236-78-154:Desktop telliott_admin$ cat /usr/local/share/EMBOSS/data/REBASE/embossre.enz ..

Repeat remap and it looks good. It would be nice if the output of enzymes with the number of sites listed the position of each cut.

# Enzymes that cut Frequency Isoschizomers Acc65I 1 Asp718I AflIII 2 ApoI 1 AcsI,XapI BbvII* 1 BpiI,BpuAI,BstV2I,BbsI BmtI 1 BspOI BsaAI 1 BstBAI,Ppu21I ..

Newton's expansion of the binomial (and more)

According to Journey Through Genius, in a letter to Leibniz Newton wrote the expression given in the first red box above. It's an expansion, not just for the binomial (a + b)², but for any power m/n.

Each of the coefficients A, B, C, D and so on consists of the entire previous term, as indicated by the definitions in the two blue boxes in the second group.

P^m/n can be factored out of each term on both sides, so we do that, yielding the expression in the second red box. The coefficients have been expanded to show the increasing powers of Q. Also each of the fractions beyond m/n can be rearranged as indicated in the bottom group. The final form isn't shown explicitly, but if we let e = m/n, we see that the coefficient for Q⁴ has the form:

e * (e-1) * (e-2) * (e-3) / 4 * 3 * 2

matching the standard form. It is trivial to verify the results for (1 + x)³, so let's think about something more interesting namely,

(1 + x)-3

substituting for m/n = -3 we obtain:

1 - 3x + (-3)(-4)/2 x2 + .. 1 - 3x + 6x2 -10x3 + ..

Expanding:

(1 + x)-3 1/(1 + 3x + 3x2 + x3)

and moving it to the right side we get:

1 = (1 + 3x + 3x2 + x3)(1 - 3x + 6x2 -10x3 + ..)

which can be checked by multiplying out to see that the rhs does equal 1.

In the same way, we can work with

√(1 - x) = (1 - x)1/2

We have Q = -x and m/n = 1/2 so that the expansion is:

1 - x/2 + x2/8 - x3/16 ..

Newton used this method to calculate roots of various kinds, e.g.

7 = 9 (7/9) = 9 (1 - 2/9) √7 = 3 (1 - 2/9)1/2

and expand as before with Q = -2/9.

As Dunham says:

In one sense, there is nothing terribly surprising about the fact that √7 can be approximated by a sum of six fractions. The truly amazing thing about this whole procedure is that Newton's binomial theorem shows us precisely which fractions to use, and generates them in an utterly mechanical fashion, devoid of the need for any particular insight or ingenuity on our part. It was a remarkably efficient and clever way to get roots of any order.

Next: an approximation for pi.

Pascal's Triangle

I'm sure you already know about Pascal's triangle. The screenshot above is from an animation in the wikipedia article.

Among other things, Pascal's triangle provides the coefficients for the binomial expansion (p + q)n. And of course, today we would generate any particular coefficient using the formula for combinations:

C(n,k) = n! / [(n-k)! k!]

So, for example, the 4th element in the expansion for (p + q)⁹ is:

9! / (5! * 4!) =

>>> 9 * 8 * 7 * 6 / (4 * 3 * 2) 126

(0-based indexing).

I'm working my way (slowly) through a wonderful book called Journey Through Genius, and I want to explore one of the chapters which involves Newton, his form of the expansion, and the connection to infinite series and finally an approximation for π.

But, let's start by simply generating the numbers as one would fill out the triangle, but using Python to do the arithmetic.

output:

1 1 1 1 2 1 1 3 3 1 1 4 6 4 1 1 5 10 10 5 1 1 6 15 20 15 6 1 1 7 21 35 35 21 7 1 1 8 28 56 70 56 28 8 1 1 9 36 84 126 126 84 36 9 1

code listing:

N = 10 L = [[1]] for i in range(1,N): sL = L[-1] R = range(len(sL)-1) pL = [sum(sL[i:i+2]) for i in R] L.append([1] + pL + [1]) sp = ' ' j = len(str(max(L[-1]))) for i,sL in enumerate(L): print sp*(N-i-1)*j, sL = [str(e).center(j) for e in sL] print (sp*j).join(sL) print

La Taqueria

Came across this pic taken in the summer on my last trip to San Francisco. They are not kidding and they're probably not wrong either. Wish I could go every day. The other place I like in the City is El Gordo. It hasn't changed since I first went there in 1982.
For the rest of us, there is only cochinita pibil (video).
Anyone know how to get the accented i in a title on Blogger?

EMBOSS

I made a stab at building EMBOSS but failed and succeeded, I think.

According to the adminguide, it requires zlib and libpng as well as gd. Since EMBOSS is set up for Linux, I will probably have to edit the config file to show where these libraries are on my machine.

I installed zlib and libpng to get matplotlib up and running (according to Gavin Huttley's wiki instructions).

Let's get gd. Downloaded gd-2.0.35.tar.gz from this page. According to the readme, it's the usual:

./configure make install

(using sudo). And it has a problem:

gdft.c:1403:35: error: fontconfig/fontconfig.h: No such file or directory gdft.c:1466: error: expected '=', ',', ';', 'asm' or '__attribute__' before '*' token gdft.c:1479: error: expected ')' before '*' token gdft.c: In function 'font_pattern': gdft.c:1516: error: 'FcPattern' undeclared (first use in this function) gdft.c:1516: error: (Each undeclared identifier is reported only once gdft.c:1516: error: for each function it appears in.) gdft.c:1516: error: 'font' undeclared (first use in this function) gdft.c:1517: error: 'FcChar8' undeclared (first use in this function) gdft.c:1517: error: 'file' undeclared (first use in this function) gdft.c:1518: error: 'pattern' undeclared (first use in this function) gdft.c:1525: warning: assignment makes pointer from integer without a cast gdft.c:1531: error: expected ')' before 'FcChar8' gdft.c:1538: error: 'FC_FILE' undeclared (first use in this function) gdft.c:1538: error: 'FcResultMatch' undeclared (first use in this function) make[1]: *** [gdft.lo] Error 1 make: *** [install-recursive] Error 1

It seems to be related to a library called Fontconfig.

Get fontconfig-2.8.0.tar.bz2 from here. According to INSTALL do:

./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man

but it fails because there is no configure, why not? From a web search it seems that

autoconf processes configure.in to produce a configure script

so I try this:

$ autoconf configure.in:36: error: possibly undefined macro: AM_INIT_AUTOMAKE If this token and others are legitimate, please use m4_pattern_allow. See the Autoconf documentation. configure.in:38: error: possibly undefined macro: AM_MAINTAINER_MODE configure.in:59: error: possibly undefined macro: AM_CONFIG_HEADER configure.in:64: error: possibly undefined macro: AC_LIBTOOL_WIN32_DLL configure.in:65: error: possibly undefined macro: AM_PROG_LIBTOOL configure.in:78: error: possibly undefined macro: AM_CONDITIONAL

and then do as before:

$ sudo ./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man ./configure: line 1728: syntax error near unexpected token `fontconfig,' ./configure: line 1728: `AM_INIT_AUTOMAKE(fontconfig, 2.8.0)'

It looks like I'm supposed to use autoconf, but first I have to fix the above error.

After another web search, try adding this to configure.in

m4_pattern_allow([AM_INIT_AUTOMAKE]) m4_pattern_allow([AM_MAINTAINER_MODE]) m4_pattern_allow([AM_CONFIG_HEADER]) m4_pattern_allow([AC_LIBTOOL_WIN32_DLL]) m4_pattern_allow([AM_PROG_LIBTOOL]) m4_pattern_allow([AM_CONDITIONAL])

autoconf goes fine. Now:

$ sudo ./configure --sysconfdir=/etc --prefix=/usr --mandir=/usr/share/man ./configure: line 1736: syntax error near unexpected token `fontconfig,' ./configure: line 1736: `AM_INIT_AUTOMAKE(fontconfig, 2.8.0)'

And I'm stuck. You can't feed fontconfig to AM_INIT_AUTOMAKE.No EMBOSS for me, it seems.

[UPDATE: I found instructions for installing GD here. It turns out I already had /usr/X11R6/include/fontconfig. The test of GD looks fine.

Followed simple instructions for EMBOSS and did:

./configure make make install

And after searching around for 10 minutes I finally ran:

/usr/local/share/EMBOSS/jemboss/runJemboss.sh

Now I have a window to play with (graphic at top of post)! ]

Normal distribution construction in Python

Here is an example with the normal distribution that will seem trivial after the t-distribution (here).

The basic form of the normal is exp {-x2/2}. We define that as a Python function f(x), vectorize it, and construct an array X of discrete points from -10 to +10 with interval dx = 0.001. We apply the vectorized function to the array to get the relative densities. So that we obtain the correct area under the curve, we multiply the height (value of f(x)) of each piece by its width, dx. When we sum up all the pieces, the total is equal to √2π as seen in the printout. We divide by this value to normalize the distribution so that its total area is equal to 1 and it becomes a pdf.

The form of the normal that includes a term for the standard deviation is normal(x). Everything is as before, except we substitute z for x, and at the end we find that the normalizing constant is 1/σ√2π. We plot it to have something pretty to look at.

[UPDATE:
In my discussions of probability distributions of late I've played a little fast and loose with terminology. The pdf (probability density function) has a value for any x. At x = μ the exponential term is equal to 1, and the value is 1/σ√2π. For the standard normal, this is about 0.4.

>>> 1.0/sqrt(2*pi) 0.3989422804014327

Of course, this is somewhat misleading, since the probability that x = any particular point approaches zero, because x is a real number with infinitely many values. The way we actually use the pdf is to ask what is the probability for a particular window, a range of values for a < x < b, and at least conceptually, this is done by integrating the function between these limits. That's where the discrete version that I developed comes in handy. We evaluate p(x) for a small enough interval and then multiply by the interval size to get a probability for that slice. To integrate between limits, just add the included slices. But for this to work, we have to "normalize" the pdf so that the total of all the little rectangles of width dx is equal to 1. It also helps to have the discrete version in developing the cdf (cumulative distribution function), since we can just accumulate the pdf as we move along the values of X generating the cdf as we go.

Technically, the pdf(x) is the slope of the cdf(x), which gets around the issue mentioned above for a continuous function. ]

output:

2.507 2.507 5.013 5.013

code listing:

from __future__ import division import numpy as np import math import matplotlib.pyplot as plt @np.vectorize def f(x): return math.e**(-0.5*(x**2)) @np.vectorize def normal(x,mu=0,sigma=1): z = (x-mu)/sigma return math.e**(-0.5*(z**2)) #================================== dx = 0.001 X = np.arange(-10,10+dx,dx) pdf = f(X) pdf *= dx print round(sum(pdf),3), S = math.sqrt(2*math.pi) print round(S,3), pdf /= S sigma = 2 pdfn = normal(X,0,sigma) pdfn *= dx print round(sum(pdfn),3), Sn = math.sqrt(2*math.pi) * sigma print round(Sn,3), pdfn /= Sn #================================== plt.plot(X,pdf,color='r',lw=4) plt.plot(X,pdfn,color='k',lw=4) ax = plt.axes() ax.set_xlim(-6,6) m = max(pdf) ymin,ymax = -m/100,m*1.1 ax.set_ylim(ymin,ymax) plt.text(3,0.75*ymax,s='$\sigma = 1$', color='r',fontsize=24) plt.text(3,0.65*ymax,s='$\sigma = 2$', color='k',fontsize=24) plt.savefig('example.png')

One more thing

One thing I forgot to do is test the two_sample_t function. Here we draw samples of different sizes (m=2,n=4) from distributions A and B (normal with the same mean and standard deviation). Since A and B have the same mean, the expected difference is 0. We record how frequently p < 0.05 and it is, as expected, about 5%. The corresponding t-statistic is -2.133, which is also expected, for one-tailed application of the t-test with 4 degrees of freedom.

output:

0.0491 -2.133 0.0516 -2.133 0.0524 -2.133 0.0506 -2.133 0.0502 -2.133

code listing:

from __future__ import division import numpy as np from t_test import two_sample_t mu,sigma = 10,3 m,n = 2,4 alpha = 0.05 N = 10000 tmax = -N for i in range(5): counter = 0 for j in range(N): A = np.random.normal(mu,sigma,m) B = np.random.normal(mu,sigma,n) t,p = two_sample_t(A,B) if p < alpha: counter += 1 if t > tmax: tmax = t print '%3.4f' % (counter/N), print round(tmax,3)

Student's t-test again 8

wikimedia

Top level: tails and errors

This is the eighth and last in a series Student's t test. The first five links are here here, plus here and here.

If you're just looking for some code, it is in the sixth post.

What's left? Error checking is important. You should look at the PyCogent module test.py for examples of that. I want to focus today on "tails."

A central idea of the frequentist approach to statistics is that the population mean and standard deviation are fixed but of unknown value. Some consequences:

• we must pick the p-value we will consider significant in advance of seeing the data (in practice, nearly all biologists use 0.05, so that p < 0.05 is signficant. The reason is that this is the largest p-value they think they can get away with!)

• we can specify in which direction we expect an expected difference to lie

Let's expand the second point. The one-sample t-test requires input of a sample and an expected mean. If we specify in advance of seeing the data that the expected mean is either lower or higher than the sample mean, we can use a "one-tailed test" where the t-statistic is significant if a value of 0.05 or less is obtained from the t-distribution. If we are uncertain or agnostic about the direction of the change, then the t-statistic is tested against 0.025, so that the difference between the sample mean and the expected mean must be larger to be significant. In my examples, I used the first choice (because it was simpler to program), as can be seen from the R versions which specified "alternative='less'."

So, to actually use my versions of the t-test, you will need top-level code to handle the tails.

Another issue that looks simple but could be tricky to code is whether the t-statistic is positive or negative.

In the one-sample t-test, if the expected mean is larger than the sample mean, then t < 0, since we do np.mean(A) - mu. If we had specified that it would be smaller than the sample mean, then we're OK ("alternative='less'"). If we specified it to be larger, then we need to switch the sign of t. Similarly with the paired t-test.

The last point is that we've made the (unspoken) assumption for the two-sample t-test that the variances are the same for the two populations from which samples were drawn. We should explicitly state whether this is the case ("var.equal=TRUE" in R), and if it is not, then the test becomes more complicated and the values from the function as given here will not be correct.

Student's t-test again 7

Interlude: calculating the sum of squares

This is the seventh in a series of posts about Student's t test. The last post has a list of all the previous posts here.

The PyCogent two-sample t-test uses an interesting method to calculate the sum of squares (X - mean(X))**2. Here is my transcription and a comparison to a more usual approach:

from __future__ import division import numpy as np np.random.seed(153) def pycogent_sumsq(X): return sum(X**2) - sum(X)**2 / len(X) def my_sumsq(X): X_bar = np.mean(X) squares = (X - X_bar)**2 return sum(squares) for i in range(5): X = np.random.normal(0,2,5) print pycogent_sumsq(X), my_sumsq(X)

We do five trials and compare the results.

15.0864873521 15.0864873521 14.2528488226 14.2528488226 8.86248689793 8.86248689793 22.3695874949 22.3695874949 37.0783435439 37.0783435439

Clearly, it works. The question is, how?
Using the code listed at the bottom of the post, examine the calculation.

pycogent_sumsq X 3.003 0.99 -2.319 5.699 -0.03 X**2 9.019 0.98 5.379 32.483 0.001 sum(X**2) 47.86 sum(X) 7.34 sum(X)**2 53.92 sum(X)**2/len(X) 10.78 sum(X**2) - sum(X)**2/len(X) 37.08 pycogent_sumsq(X) 37.08 my_sumsq X - X_bar 1.535 -0.478 -3.788 4.231 -1.499 (X - X_bar)**2 2.355 0.229 14.348 17.9 2.247 sum((X - X_bar)**2) 37.08 my_sumsq(X) 37.08

But I am stuck in the process of figuring out why they are equivalent.

(1) sum(X**2) = x12 + x22 + ... (2) sum(X)**2/n = (x1 + x2 + ... )2 / n m = mean(X) = sum(X) / n = (x1 + x2 + ... ) / n (1) - (2) = (x12 + x22 + ... ) - m (x1 + x2 + ... ) = x12 - mx1 + x22 - mx2 + ... = x1(x1 - m) + x2(x2 - m) + ...

I don't seem to be able to reach:

= (x1 - m)2 + (x2 - m)2 + ...

Any ideas?

[UPDATE: I got an answer from Srikant Vadali here, and will discuss it sometime soon.]

Remainder of the code listing:

def show(s,f): print s.ljust(30), try: print str(round(f,2)).rjust(6) except TypeError: for n in f: print str(round(n,3)).rjust(6), print print print 'pycogent_sumsq' show('X', X) show('X**2', X**2) show('sum(X**2)', sum(X**2)) show('sum(X)', sum(X)) show('sum(X)**2', sum(X)**2) show('sum(X)**2/len(X)',sum(X)**2/len(X)) show('sum(X**2) - sum(X)**2/len(X)', sum(X**2) - sum(X)**2 / len(X)) show('pycogent_sumsq(X)',pycogent_sumsq(X)) print print 'my_sumsq' X_bar = np.mean(X) show('X - X_bar', X - X_bar) show('(X - X_bar)**2', (X - X_bar)**2) show('sum((X - X_bar)**2)', sum((X - X_bar)**2)) show('my_sumsq(X)', my_sumsq(X))

Student's t-test again 6

Code

This is the sixth in a series of posts about Student's t test. The previous posts are here, here, here, here, and here.

In this post I'm going to talk about how I would test my home-grown functions to do the t-test in Python. I'm only going to do this for the one-sample test, because I'm running out of energy. But you could easily extend this to the other functions.

We could check using SciPy if it is installed, and that would be easy. But I'm going to use R, partly because it's the gold standard for statistics, and partly because it's a useful technique that I've used in my real work. While we could use RPy as I mentioned the other day, here I am going to use R to generate the sample data, run the test of interest and then write the results to disk.

What we're going to do is run R in batch mode. If it gives an error, which frequently happens when building such a test, we want to show the error. We also clean up extra files when we're done. The code to do this has been added to utils.py and is given below. It should be self-explanatory.

When we run R this way, we feed it a text file with the commands to run. This file is constructed by the function write_R_code. It's not pretty, but I hope you get the idea. Perhaps a better method, which I've also used, is to write the code to a text file, then load it into Python, substitute the file names and specific functions etc., then write that to disk as the file to give to R.

In the last step, we load the results that R has given us and evaluate how well our functions do on the same data. As you can see, we match very well as a rule, but have rounding errors at the extremes of the t-distribution. I believe this could be fixed by extending the range over which we built the distribution (see here). I think it looks pretty good.

[UPDATE: Just to make it clear what we did: we do 1000 runs and look at the details for the first five, then test the rest for deviation from the R results by more than ε = 0.0001. ]

Output:

run (1626, 0) clean .RCode removed .RCode.Rout removed .RData removed .RHistory not found test results: N = 1000 R gave: t = -1.858651 I get: t = -1.85865011386 R gave: t = -4.250644 I get: t = -4.25063601226 R gave: t = -1.743078 I get: t = -1.74307933578 R gave: t = 0.1286087 I get: t = 0.128612751636 R gave: t = -9.36424 I get: t = -9.36429217504 error with: [ 8.946167 8.941678 8.940817] -1843.026 -1843.23012185 error with: [ 8.895525 8.686258 8.761813] -52.60937 -52.6095381366 error with: [ 10.54477 10.42989 10.66126] -21.77912 -21.7796812829 error with: [ 11.05575 10.89208 10.68675] -10.50883 -10.5090054298 error with: [ 10.02618 10.07471 10.16116] -48.45138 -48.4527009573

R_test.py

import numpy as np from t_test import one_sample_t import utils def write_R_code(D): L = [ 'fn <- file("' + D['results_fn'] + '","w")', 'n = ' + str(D['n']), 'mu = ' + str(D['mu']), 'sigma = ' + str(D['sigma']), 'for (i in 1:' + str(D['N']) + ') {', 'A = rnorm(n,mu,sigma)', 'result = ' + D['func'] + str(D['mu_est']) + ')', 'cat(A,file = fn,sep = "\n")', 'cat(result$statistic,file = fn,sep="\n")', 'cat("\n",file = fn)', '}', 'close(fn)' ] FH = open(D['code_fn'], 'w') FH.write('\n'.join(L)) FH.close() def compare_to_my_result(L): print 'test results:' print 'N = ', len(L) for i,e in enumerate(L): A = np.array(e[:-1]) result = one_sample_t(A,mu=12) if i < 5: print 'R gave: t =', e[-1] print 'I get: t =', result[0] print else: if e[-1] - result[0] > 0.0001: print 'error with:' print A, e[-1], result[0] print if __name__ == '__main__': directory = '/Users/telliott_admin/Desktop/' D = { 'code_fn':directory + '.RCode', 'results_fn':directory + 'results.txt', 'func':'t.test(A,alternative="less",mu=', 'n':3,'mu':10,'sigma':2,'mu_est':12,'N':1000 } write_R_code(D) utils.runR(D['code_fn'],v=True) utils.cleanup_R_files(D['code_fn'],directory,v=True) L = utils.load_data(D['results_fn']) compare_to_my_result(L)

utils.py

from __future__ import division import os,subprocess import numpy as np def unbiased_var(X): n = len(X) sample_SS = sum(X**2) - sum(X)**2 / n return sample_SS/ (n-1) def unbiased_std(X): return np.sqrt(unbiased_var(X)) def runR(code_fn,v=False): if v: print 'run' cmd = 'R CMD BATCH ' + code_fn p = subprocess.Popen(cmd, shell=True) sts = os.waitpid(p.pid, 0) if v: print sts if sts[1] != 0: show_R_error() def show_R_error(): fn = directory + '.Rcode.Rout' FH = open(fn, 'r') msg = FH.read() FH.close() print '\n'.join(msg.split('\n')[-3:]).strip() def cleanup_R_files(code_fn,directory,v=False): if v: print 'clean' L = ['.RCode.Rout','.RData', '.RHistory'] L = [code_fn] + [directory + e for e in L] for fn in L: try: os.remove(fn) if v: print fn.split('/')[-1], 'removed' except OSError: if v: print fn.split('/')[-1], 'not found' print def load_data(results_fn): def convert(s): L = s.split('\n') return [float(n) for n in L] FH = open(results_fn, 'r') data = FH.read().strip() FH.close() L = data.split('\n\n') rL = [convert(s) for s in L] return rL

Student's t-test again 5

Code

This is the fifth in a series of posts about Student's t test. The previous posts are here , here, here, and here.

Last time I showed a version of the one-sample t-test in Python. As we saw, the t-statistic is constructed similar to a standard error (like the sem). It is something like a z-score using the expected value for the population mean, and the observed (unbiased) sample standard deviation, then adjusted by multiplying this "z-score" by √n. (The same function, slightly edited, is also given below).

The paired t-test is simply a one-sample t-test done on the differences between paired values, comparing them to the expected difference rather than the expected mean. The paired values arise, for example, in repeated tests of the same individuals. The idea is that variance due to differences between individuals is the same for each value in a single pair.

The two-sample t-test takes account of the expected difference in means:

diff = (np.mean(A) - np.mean(B) - expected_diff)

The code I copied (and massaged) takes the observed sample variance and turns it back into the sum of squares:

sum_sq = (var(A)*(na-1) + var(B)*(nb-1))

There's a super-duper factor which is constructed from the number of values in each sample (and the dependent degrees of freedom):

na = len(A) nb = len(B) df = na + nb - 2 f = (1/na + 1/nb)/df t = diff/np.sqrt(sum_sq*f)

The t-statistic is as shown. I offer no justification for these procedures. (If you want more theory and don't like your book, start here: pdf). What we will do is test the output from these functions. Today we simply check that R gets the same t-statistic for one (or two) examples. Next time we will construct a test harness that does more. And after that we may do some more testing that shows these functions do give the correct distribution with simulated data.

I feel like I do understand the tests better after this project.

The functions for unbiased variance and standard deviation in utils.py were given last time. There's something a bit unusual about the variance calculation, which comes from PyCogent. I'll try to remember to deal with that when I talk about error checking, and "tails", in a future post. The functions here are just the basic t-statistic calculations themselves.

The source module for the t distribution:

from transcendental import stdtr

was explained here.

Here is the output if the module is run as "__main__":

t-statistic p-value -2.378 0.049 -6.455 0.004 -2.611 0.040 -2.574 0.021

R output using the code as given at the end of the listing (and as shown here):

> A = c(3.1,2.3,2.1,1.7) > result = t.test(A, alternative='less',mu=3) > result$statistic t -2.377782 > B = c(2.1,1.8,2.7,2.4) > result=t.test(B, alternative='less',mu=3.5) > result$statistic t -6.454972 > C = c(3.1,4.3,4.1,2.7) > result=t.test(A,C,alternative='less', + paired=TRUE,var.equal=FALSE) > result$statistic t -2.611165 > result=t.test(A,C,alternative='less',var.equal=TRUE) > result$statistic t -2.573993 >

code listing:

from __future__ import division import numpy as np from transcendental import stdtr from utils import unbiased_var as var from utils import unbiased_std as std def one_sample_t(A,mu): n = len(A) df = n-1 z = (np.mean(A) - mu) / std(A) t = z * np.sqrt(n) return t, stdtr(df,t) def paired_t(A,B,expected_diff=0): return one_sample_t(A - B,expected_diff) def two_sample_t(A,B,expected_diff=0): diff = (np.mean(A) - np.mean(B) - expected_diff) na = len(A) nb = len(B) df = na + nb - 2 sum_sq = (var(A)*(na-1) + var(B)*(nb-1)) f = (1/na + 1/nb)/df t = diff/np.sqrt(sum_sq*f) return (t, stdtr(df,t)) #========================================== def test(func,D): if func == one_sample_t: result = one_sample_t(D['X'],D['mu']) elif func == paired_t: result = paired_t(D['X'],D['Y']) elif func == two_sample_t: result = two_sample_t(D['X'],D['Y']) if D['verbose']: print '%5.3f %5.3f' % result return result if __name__ == '__main__': print 't-statistic p-value' A = np.array([3.1,2.3,2.1,1.7]) B = np.array([2.1,1.8,2.7,2.4]) test(one_sample_t,{'X':A,'mu':3,'verbose':True}) test(one_sample_t,{'X':B,'mu':3.5,'verbose':True}) C = np.array([3.1,4.3,4.1,2.7]) test(paired_t,{'X':A,'Y':C,'verbose':True}) test(two_sample_t,{'X':A,'Y':C,'verbose':True}) ''' R code: A = c(3.1,2.3,2.1,1.7) result = t.test(A, alternative='less',mu=3) result$statistic B = c(2.1,1.8,2.7,2.4) result=t.test(B, alternative='less',mu=3.5) result$statistic C = c(3.1,4.3,4.1,2.7) result=t.test(A,C,alternative='less', paired=TRUE,var.equal=FALSE) result$statistic result=t.test(A,C,alternative='less',var.equal=TRUE) result$statistic '''

Student's t-test again 4

One Sample t-test

This is the fourth in a series of posts about Student's t test. The previous posts are here , here and here.

There are several different versions of the t-test, but the simplest is the one sample test. We will have a single sample with a (relatively) small number of values. We calculate the mean and the variance and then the standard deviation of the sample values. Importantly, the variance and standard deviation are the unbiased versions, in which the sum of squares is divided by n-1 rather than n.

It's easy to see that we should not use numpy's var for this test:

>>> import numpy as np >>> A = [1,2,3] >>> np.var(A) 0.66666666666666663

As help(np.var) indicates:

Notes ----- The variance is the average of the squared deviations from the mean, i.e., var = mean(abs(x - x.mean())**2). The computed variance is biased, i.e., the mean is computed by dividing by the number of elements, N, rather than by N-1.

In the code below we do:

z = np.mean(A) - mu z /= unbiased_std(A)

So, at least conceptually, at this point we have a z-score for the sample mean. With a standard z-score, we would compare it to the normal distribution. Here we do two things differently: we multiply by √n (the square root of the number of samples), and we find the value of the resulting t statistic in the t distribution.

t-statistic p-value -2.378 0.049

The result indicates that, given a prior choice of a limit of 0.05 (for this one-sided test), the null hypothesis that the mean of the sample values is equal to 3 is not supported (just barely).

You can read much more about the background of the test here.

The second set of values in the code below gives:

t-statistic p-value -6.455 0.004

It's reassuring that a One Sample t-test in R gives a similar result. I'll have more to say about "tails" and sidedness in another post.

> A = c(3.1,2.3,2.1,1.7) > t.test(A, alternative='less',mu=3) One Sample t-test data: A t = -2.3778, df = 3, p-value = 0.04891 alternative hypothesis: true mean is less than 3 95 percent confidence interval: -Inf 2.992811 sample estimates: mean of x 2.3 > B = c(2.1,1.8,2.7,2.4) > t.test(B, alternative='less',mu=3.5) One Sample t-test data: B t = -6.455, df = 3, p-value = 0.003771 alternative hypothesis: true mean is less than 3.5 95 percent confidence interval: -Inf 2.705727 sample estimates: mean of x 2.25 >

Obviously, we could use a lot more testing.

code listing:

from __future__ import division import numpy as np from transcendental import stdtr from utils import unbiased_std def one_sample_t(A,mu): n = len(A) df = n - 1 z = np.mean(A) - mu z /= unbiased_std(A) t = z * np.sqrt(n) return t, stdtr(df,t) def test(A,mu): result = one_sample_t(A,mu) print 't-statistic p-value' print '%5.3f %5.3f' % result if __name__ == '__main__': A = np.array([3.1,2.3,2.1,1.7]) B = np.array([2.1,1.8,2.7,2.4]) test(A,mu=3) test(B,mu=3.5) ''' R code: A = c(3.1,2.3,2.1,1.7) t.test(A, alternative='less',mu=3) B = c(2.1,1.8,2.7,2.4) t.test(B, alternative='less',mu=3.5) '''

utils.py

from __future__ import division import numpy as np def unbiased_var(X): n = len(X) sample_SS = sum(X**2) - sum(X)**2 / n return sample_SS/ (n-1) def unbiased_std(X): return np.sqrt(unbiased_var(X))

Student's t-test again 3

This is the third in a series of posts about Student's t test. The previous posts are here and here.

While looking at the code in PyCogent I noticed a reference to "Cephes." The Cephes Mathematical Library (here)

It contains lots and lots of functions (like the t distribution and the incomplete beta and so on). The analogous functions in PyCogent appear to be just straight translations of the C code.

I found a Python module that makes the Cephes library available directly. It was put together by Michiel de Hoon (here, old home page here), when he was at the University of Tokyo. I've also looked at his clustering software before (here).

According to the last page, he's currently at Columbia University. Just two seconds ago, I saw this, some kind of manual he's written for doing statistics in Python. Whoaah! Looks like it should be very interesting!

Anyhow, I downloaded the source and did

$python setup.py config ERROR: unknown if big-endian or little-endian edit the config file before compiling

The config file looks like this

/* config.h file created by setup.py script Fri Nov 5 15:55:52 2010
* run on darwin */
#define UNK 1
#define BIGENDIAN UNKNOWN /* replace with 0 or 1 */

I modified it appropriately and then did as the instructions say:

python setup.py build python setup.py install

test by:

>>> from transcendental import ndtr >>> ndtr(2.) - ndtr(-2.) 0.95449973610364158

I figured the C version should be a lot faster, and it is:

$ python -m timeit -s 'from transcendental import stdtr; import random' 't = random.random(); stdtr(3,t)' 1000000 loops, best of 3: 1.26 usec per loop $ python -m timeit -s 'import distribution; from distribution import stdtr; import random' 't = random.random(); stdtr(3,t)' 10000 loops, best of 3: 24.2 usec per loop

But ... it doesn't seem to be the slow step in the t test code.

---

Calculating the cdf and testing the result

I looked up the definitions for the t distribution's pdf and cdf in wikipedia (here).

I translated the x-dependent part into Python and "vectorized" and used it to constuct a pdf. It's important to set really wide limits for the vector in this step, otherwise numerous small probability values are discarded, which then leads to small errors in the eventual product.

Since the normalization constants (gamma functions) depend on df but not on x, I was able to normalize the pdf simply by doing:

pdf /= sum(pdf)

I also used a very small step size in the Numpy array. The second section of the code listing ends by printing values from the distributions for whole numbers between -5 and 5:

pdf -5.000 0.00422 -4.000 0.00916 -3.000 0.02297 -2.000 0.06751 -1.000 0.20675 0.000 0.36755 1.000 0.20675 2.000 0.06751 3.000 0.02297 4.000 0.00916 5.000 0.00422 calculated cdf -5.000 0.00770 -4.000 0.01401 -3.000 0.02884 -2.000 0.06970 -1.000 0.19560 0.000 0.50018 1.000 0.80460 2.000 0.93037 3.000 0.97118 4.000 0.98600 5.000 0.99231

As a test, I show that I get very similar values using stdtr from either PyCogent or Cephes:

stdtr from transcendental -5.000 0.00770 -4.000 0.01400 -3.000 0.02883 -2.000 0.06966 -1.000 0.19550 0.000 0.50000 1.000 0.80450 2.000 0.93034 3.000 0.97117 4.000 0.98600 5.000 0.99230 stdtr from PyCogent -5.000 0.00770 -4.000 0.01400 -3.000 0.02883 -2.000 0.06966 -1.000 0.19550 0.000 0.50000 1.000 0.80450 2.000 0.93034 3.000 0.97117 4.000 0.98600 5.000 0.99230

And that all matches what I get from R. So I'm pretty confident that I have the distributions correct and the functions are working as they should.

> A = seq(-5,5,by=1) > pt(A,3) [1] 0.007696219 0.014004228 0.028834443 0.069662984 [5] 0.195501109 0.500000000 0.804498891 0.930337016 [9] 0.971165557 0.985995772 0.992303781

code listing:

import numpy as np import distribution # also imports special import transcendental dx = 0.001 X = np.arange(-100,100+dx,dx) @np.vectorize def p(df,x): base = 1 + x**2/df exponent = (df+1)/2 return base**(-exponent) df = 3 pdf = p(df,X) # normalization depends on df but not x # so this works pdf /= sum(pdf) #============================================= def show(F,multiplier=1): N = 95000 for n in range(N,len(X)-N,1000): t = X[n], multiplier*F[n] print '%5.3f %6.5f' % t print # we include a multiplier for the pdf only print 'pdf' show(pdf,multiplier=1000) cdf = [pdf[0]] for n in pdf[1:]: cdf.append(n + cdf[-1]) print 'calculated cdf' show(cdf) #============================================= @np.vectorize def f(x): return transcendental.stdtr(df,x) cdf2 = f(X) print 'stdtr from transcendental' show(cdf2) @np.vectorize def f(x): return distribution.stdtr(df,x) cdf3 = f(X) print 'stdtr from PyCogent' show(cdf3)

Student's t-test again 2

This is the second in a series of posts about Student's t test. The previous post is here. If you want to do a t-test for real data, I would urge you to use either R or SciPy.

---

R

R is great software. It was written by statisticians for statistical work and is thoroughly tested. My only problem is that I find it very difficult to write R code. R was not designed for text or character manipulation, where Python excels. I've been very impressed with how well the GUI works on Mac OS X.

Note: I posted about the t-test in R before (here, here and here).

---

SciPy

SciPy is a library that provides all kinds of goodies useful for scientific applications. It has statistical functions too, and lots of 'em, e.g. the two sample t test . If your goal is to use Python for statistical programming, that is probably where you should go. I posted about my difficulties installing SciPy: here, here and here. But I finally got it to work, even on my office computer, described briefly here. Resist the urge to "roll your own."

You could also try Sage. It should be much easier to install than SciPy (a binary with everything but the kitchen sink). However, I could only find elementary statistics in a quick search.

---

PyCogent

Our goal here is to understand how the tests work. It's not complicated. In order to do this I started with the statistical functions included as part of PyCogent. I have a bunch of posts about PyCogent here, including the one on the Two sample t-test in Python here.

If you don't want to install the whole PyCogent package (or you just want to take the modules apart like I did) you can download (download source) and modify three modules from the source: these are test.py, distribution.py and special.py. You can find them in /PyCogent-1.4.1/cogent/maths/stats. I copied them out onto my desktop, and then I stripped out everything I could, while still allowing my test of t_two_sample to run.

Let's look at the names defined in the stripped-down modules:

>>> def show(module): ... for n in dir(module): ... if not '_' == n[0]: ... print n ... >>> import test >>> show(test) array asarray division isinf isnan mean sqrt sum t_high t_low t_one_observation t_one_sample t_paired t_tailed_prob t_two_sample tprob var >>> import distribution >>> show(distribution) MACHEP PI atan betai division exp sqrt stdtr stdtri t_high t_low tprob >>> import special >>> show(special) GP GQ Gamma Gamma_small MACHEP MAXGAM MAXLOG MINLOG betai betai_result division exp floor log polevl pseries sin sqrt

The modules are hierarchical. At the top level (test.py) are the functions we would call from our scripts, including t_two_sample. These call down into distribution.py to stdtr, which is the standard t distribution. This distribution is in turn computed using the integral of the incomplete beta function (betai) and the Gamma function in special.py.

stdtr is a cumulative distribution function. I used a trick to compute the pdf from it, and then plot the pdf for various values of the degrees of freedom (df). This figure is at the top of the post. As you can see, the smaller the df, the "fatter" the tails on the distribution. I made a second plot for the cdf.



script.py

import math import numpy as np import matplotlib.pyplot as plt from distribution import stdtr PI = math.pi # needed by stdtr dx = 0.01 X = np.arange(-4,4+dx,dx) colors = 'rkkkm' plot_pdf = True for i,df in enumerate([1,2,5,10,30]): cdf = [stdtr(df,x) for x in X] pdf = [0] for j in range(1,len(cdf)): pdf.append(cdf[j] - cdf[j-1]) Y2 = np.array(pdf) Y2 /= sum(Y2) if plot_pdf: plt.plot(X,Y2,color=colors[i],lw=2) else: plt.plot(X,cdf,color=colors[i],lw=2) plt.savefig('example.png')

ttest.py

# two sample t-tests import numpy as np import test as stats def oneTrial(n,f=np.random.normal): N = 50 # 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 n = 3 # sample size L = list() for i in range(10): #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 from the second script:

avg 0.0515, std 0.0094 min 0.0300, max 0.0720

The output shows that on the average, the t statistic is < 0.05 about 5% of the time, as expected. Compare to the results from here.

Still to come: computing the t distribution, and setting up the tests.

Student's t-test again 1

Periodically, I've wondered about whether anyone reads these posts. In a meta sort of way, I've even posted about it.

As ABOUT ME says, a major objective here is to provide a record so that "google can organize my head." (bbum now says "index"---well, OK)

Still, one wonders. And I have to admit that I'm NOT repeat NOT an expert on most of these topics. I'm just a student, and a true believer in Python. Although I'm sure Ruby would satisfy me as well (if I hadn't met Python first), but I think the possible advantage would be marginal.

And, please, if you find an error here (or even something interesting), speak up.

So I was excited to discover this AM that blogger has stats. You can't see my stats, but I can. And the stats say that 166 people loaded pages from my blog today. 166 != 0. That's great. A perennial winner is a post about t-tests.

That smells like homework. Nevertheless, it put a thought in my tiny brain, and it won't go away. So I started thinking about the t-test: the what, the why, the how, and the Python.

Note: I posted about the t-test in R before (here, here and here).

Anyway, I'm going to look into the subject a bit more seriously. But for starters, let's model the basic problem. If we sample repeatedly from a population of known mean and variance, then the accuracy of the observed means (the standard error of the mean) depends markedly on sample size. The normal distribution of the set of means is guaranteed by the miraculous Central Limit Theorem (here, and wikipedia).

The t-test deals with samples of small size, which are ubiquitous in classical biology. And even in microarray studies, if you consider the data gene-by-gene. Of course there is more to it, for example, whether the size of the samples is the same, and their individual variances.

Here is a simple Python simulation that shows the problem. The first set of samples is from a normal distribution, the second from a uniform distribution. You can clearly see the dependence of the standard deviation of the sample means on the sample size. The column labeled stdev is the computed standard deviation for the observed means, and the column labeled sem is the predicted standard deviation based on the population standard deviation and n:
sem = σ / √n.

num mean stdev sem 2 10.09 3.56 3.54 3 9.99 2.86 2.89 4 10.03 2.53 2.50 5 9.97 2.23 2.24 10 10.01 1.58 1.58 15 10.02 1.29 1.29 20 10.02 1.11 1.12 25 10.01 1.00 1.00 30 10.00 0.91 0.91 2 10.01 4.08 3 10.02 3.33 4 10.00 2.90 5 10.01 2.57 10 9.99 1.80 15 10.00 1.48 20 9.99 1.28 25 10.00 1.16 30 10.01 1.06

Code listing:

import numpy as np import random, math N = int(1e6) mu = 10 sigma = 5 A_normal = np.random.normal(loc=mu,scale=sigma,size=N) A_uniform = np.random.uniform(low=0,high=mu*2,size=N) def sample(A,SZ): rL = list() for j in range(10000): L = [random.choice(A) for k in range(SZ)] rL.append(np.mean(L)) return (np.mean(rL),np.std(rL)) def show(A,extra=True): R = range(2,6) + range(10,31,5) for SZ in R: t = sample(A,SZ) pL = ['%2d' % SZ, '%5.2f' % t[0], '%5.2f' % t[1]] if extra: pL.append('%5.2f' % (sigma/math.sqrt(SZ))) print (' '*5).join(pL) print print 'num mean stdev sem' show(A_normal) show(A_uniform,extra=False)

R from Python, baby steps

I've been playing around with RPy a bit this morning. As the main page says:

rpy2 is a redesign and rewrite of rpy. It is providing a low-level interface to R, a proposed high-level interface, including wrappers to graphical libraries, as well as R-like structures and functions.

I just used easy_install

$ easy_install rpy2 Searching for rpy2 Reading http://pypi.python.org/simple/rpy2/ Reading http://rpy.sourceforge.net Best match: rpy2 2.1.7 Downloading http://pypi.python.org/packages/source/r/rpy2/rpy2-2.1.7.tar.gz#md5=e8e8db05f13644ce04784888156af471 Processing rpy2-2.1.7.tar.gz ... error: /Library/Python/2.6/site-packages/easy-install.pth: Permission denied

For some reason, root was the owner of the .pth file. So I changed it, and then got:

Using /Library/Python/2.6/site-packages/rpy2-2.1.7_20101104-py2.6-macosx-10.6-universal.egg Processing dependencies for rpy2 Finished processing dependencies for rpy2

The example I chose to run was described in more detail in this post. If we run it from R, it looks like this:

> library(Bolstad) Warning message: package 'Bolstad' was built under R version 2.10.1 > result = binobp(68,200,1,1) Posterior Mean : 0.3415842 Posterior Variance : 0.0011079 Posterior Std. Deviation : 0.0332852 Prob. Quantile ------ --------- 0.005 0.2591665 0.01 0.2666906 0.025 0.2779134 0.05 0.287724 0.5 0.3410604 0.95 0.3972323 0.975 0.4082264 0.99 0.4210788 0.995 0.4298666 > result$mean [1] 0.3415842 > class(result$mean) [1] "numeric" >

In R, the variable result is a list of numeric vectors with names:
$posterior, $likelihood, $prior, $pi (990 elements each), $mean, $var, $sd, $quantiles.

In the Python interpreter:

>>> import rpy2.robjects as robjects >>> robjects.r['pi'][0] 3.1415926535897931 >>> >>> from rpy2.robjects.packages import importr >>> importr('Bolstad') Warning message: package 'Bolstad' was built under R version 2.10.1 >>> >>> binobp = robjects.r['binobp'] >>> result = binobp(68,200,1,1) Posterior Mean : 0.3415842 Posterior Variance : 0.0011079 Posterior Std. Deviation : 0.0332852 Prob. Quantile ------ --------- 0.005 0.2591665 0.01 0.2666906 0.025 0.2779134 0.05 0.287724 0.5 0.3410604 0.95 0.3972323 0.975 0.4082264 0.99 0.4210788 0.995 0.4298666

and it opens X11 (rather than Quartz, not sure why) to do the plot. Getting the individual values from the result is a slight pain, but not too bad:

>>> result <Vector - Python:0x100544290 / R:0x100c6c610> >>> L = str(result).split('\n') >>> L[0] '$posterior' >>> L[1][:35] ' [1] 1.516820e-80 8.663662e-78 ' >>> result[0][0] 1.5168201745820013e-80

To get the names of the vectors, we need to parse str(result) after splitting on double newlines. If you already know the index of the value you want you can just grab it directly as shown. And of course, that's better, since the value is a float rather than a string.

[UPDATE: As the first comment says, using names is the way to do this. Docs here. And see examples in later posts.]

>>> L = str(result).split('\n\n') >>> str(L[0]).split('\n')[0] '$posterior' >>> str(L[0]).split('\n')[1][:35] ' [1] 1.516820e-80 8.663662e-78 ' >>> str(L[4]).split('\n')[:2] ['$mean', '[1] 0.3415842'] >>> result[4][0] 0.34158415841584161 >>> type(result[4][0]) <type 'float'>

On the value of silicon

Computers are no more a threat to mathematicians than food
processors are a threat to cooks.

Andrei Okounkov as quoted here

Poisson simulation

This is a very simple simulation of the Poisson distribution using Python.

First a function event() is defined which returns either 1 or 0, with 1 occuring at frequency f (here f = 0.1). A list is constructed by calling event() N times.

Then the list is chopped into groups using a method described and discussed previously (here and at Stack Overflow here).

In the example, there are 10 events in each group, so the expectation for the average number of successes is 1. This should generate a Poisson distribution with mean λ = 1.

The Counter class is used to evaluate the results. It is included in Python 2.7 but not in my 2.6 installation. Then we output the results for various k:

0 3533 1 3801 2 1959 3 578 4 108 5 19 6 2

The problem I'm having is that the results are not quite correct. We obtain the expected values when n = 9 rather than n = 10. If f is made smaller (say f = 0.01) it looks better, but still the expected distribution is more closely approximated by n = 99 than n = 100.

I'd be grateful for any ideas about what's wrong!

[UPDATE: I think the problem here is the same as I had the other day. The Poisson is an approximation. If we do int(1e7) events in groups of n=1000, it looks as expected. ]

from itertools import izip_longest import numpy as np import Counter def event(): f = 0.1 r = np.random.random() if r <= f: return 1 return 0 N = int(1e5) L = [event() for i in range(N)] def groups(iterable, n=3, padvalue=0): "groups('abcde', 3, 'x') --> ('a','b','c'), ('d','e','x')" return izip_longest(*[iter(iterable)]*n, fillvalue=padvalue) rL = [sum(g) for g in groups(L,n=10)] C = Counter.Counter(rL) for i in range(max(C.keys())+1): print str(i).rjust(2), C[i]

Don't forget D.C.

I spend much more time on politics than I really should. More precisely, on trying to "suss out" the remarkably large differences between what politicians profess to believe and what they actively support. Also, to try to understand which policies are good ones and what the prospects are for those to be enacted.

One of my favorite sources is Nate Silver's blog FiveThirtyEight, even though he's been absorbed into the currently undistinguished national newspaper of record.

The question arises: why 538?

Because that is the number of electors in the Electoral College. OK. And why is that not equal to the number of Congress critters (435) + Senators (100)? It isn't because of the non-voting delegates (American Samoa, the District of Columbia, the Virgin Islands, Guam, and the Northern Mariana Islands) and one Commissioner (Puerto Rico).

It's because the 23rd Amendment to the U.S. Constitution says:



wikimedia