We introduce numpy and show some of its properties. Then we solve the birthday problem and simulate Chuck-A-Luck.

In [1]:

%pylab inline

Populating the interactive namespace from numpy and matplotlib

Beware when doing calculations involving division, especially using integers.

In [251]:

print 1/6, 1//6, 1./6, 1/6., 5/4, 5//4, 5./4, 5/4., 5.//4, 5.//4.

0 0 0.166666666667 0.166666666667 1 1 1.25 1.25 1.0 1.0

Example Numpy Calculations¶

The import below isn't necessary as the pylab command above imports everything from both numpy and pylab. We include just to show the command.

In [177]:

from numpy import *

Vectors and matrices are numpy arrays.

In [186]:

l = [1,2,3] #python list
a = array(l)
a

Out[186]:

array([1, 2, 3])

In [216]:

print a.shape, a.size, a.dtype, a.ndim

(3,) 3 int64 1

In [191]:

z = zeros(3)
x = ones(3)
z,x

Out[191]:

(array([ 0.,  0.,  0.]), array([ 1.,  1.,  1.]))

Arithmetic is elementwise.

In [222]:

a+1, a*2, a+x, a*a, exp(a), a>2, a*(a>2)

Out[222]:

(array([2, 3, 4]),
 array([2, 4, 6]),
 array([ 2.,  3.,  4.]),
 array([1, 4, 9]),
 array([  2.71828183,   7.3890561 ,  20.08553692]),
 array([False, False,  True], dtype=bool),
 array([0, 0, 3]))

In [223]:

A = arange(12).reshape(4,3)
print A
print A.shape, A.size, A.dtype, A.ndim

[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
(4, 3) 12 int64 2

In [195]:

2*A

Out[195]:

array([[ 0,  2,  4],
       [ 6,  8, 10],
       [12, 14, 16],
       [18, 20, 22]])

Broadcasting: if possible, smaller arrays extended to bigger array.

In [228]:

A+a, A*a

Out[228]:

(array([[ 1,  3,  5],
       [ 4,  6,  8],
       [ 7,  9, 11],
       [10, 12, 14]]),
 array([[ 0,  2,  6],
       [ 3,  8, 15],
       [ 6, 14, 24],
       [ 9, 20, 33]]))

Multiplication is element-wise. To get the linear algebra multiplication, use the dot function.

In [196]:

A*A

Out[196]:

array([[  0,   1,   4],
       [  9,  16,  25],
       [ 36,  49,  64],
       [ 81, 100, 121]])

In [197]:

dot(A,x)

Out[197]:

array([  3.,  12.,  21.,  30.])

In [257]:

dot(A.T,A)

Out[257]:

array([[126, 144, 162],
       [144, 166, 188],
       [162, 188, 214]])

In [260]:

X = eye(3)+dot(ones(3),ones(3).T)
X

Out[260]:

array([[ 4.,  3.,  3.],
       [ 3.,  4.,  3.],
       [ 3.,  3.,  4.]])

In [261]:

inv(X), dot(X,inv(X))

Out[261]:

(array([[ 0.7, -0.3, -0.3],
       [-0.3,  0.7, -0.3],
       [-0.3, -0.3,  0.7]]),
 array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [ -2.22044605e-16,   1.00000000e+00,   0.00000000e+00],
       [ -2.22044605e-16,   0.00000000e+00,   1.00000000e+00]]))

In [282]:

n = 365
days = arange(1,n+1)
p = 1.0*(n+1-days)/n
probnopair = cumprod(p)
sum(probnopair>0.5)

Out[282]:

In [285]:

xlim = 30
cutoff = 23
plot(days[:xlim],probnopair[:xlim],days[:xlim],probnopair[:xlim],'.')
axvline(x=cutoff,color='r')
axhline(y=0.5)
xlabel('Number of People in Group')
ylabel('Probability of No Match')

Out[285]:

<matplotlib.text.Text at 0x111cb2750>

In [284]:

probnopair[21],probnopair[22] #python arrays start from 0, not 1.

Out[284]:

(0.5243046923374497, 0.4927027656760144)

Chuck-A-Luck is a dice game that is sometimes played at carnivals and charity events. The player pays 1 USD, selects a number between 1 and 6, and rolls 3 dice. If the player's number comes up once, the player wins 1 USD; if it comes up twice, the player wins 2 USD; and if thrice, wins 3 USD (in one variation, 10 USD).

(Sorry for the European notation, but there doesn't seem to be a way of getting dollar signs properly displayed. Ipython uses dollar signs to delimit inline math, e.g., \(x^2+y^2=z^2\).)

In [272]:

dice = 3
throws = 1000
p = 1./6

In [274]:

brv = stats.bernoulli.rvs
rolls = brv(p,size=dice*throws).reshape(throws,dice)
counts = sum(rolls,1)

We use what python calls "fancy indexing" to transform the number of correct dice in each throw to the amount won or lost on each throw. It uses an implied for loop over the entries in counts. If the i'th entry in counts is k, then payouts[k] is selected. (You can also use the numpy take function to get the same effect.)

In [275]:

payouts = array([-1,1,2,10])
#wins = take(payouts,counts) #alternative function
wins = payouts[counts]

In [276]:

totals = cumsum(wins)
plot(totals)
axhline(y=0.0)

Out[276]:

<matplotlib.lines.Line2D at 0x110ed15d0>

In [277]:

averages = 1.0 * totals/arange(1,throws+1)
plot(averages)
axhline(y=0.0)

Out[277]:

<matplotlib.lines.Line2D at 0x111b67410>

In [278]:

averages[-1]

Out[278]:

-0.047

Birthday and Chuck-A-Luck¶

Python Division¶

Example Numpy Calculations¶

Birthday Problem¶

Chuck-A-Luck¶