P4DA NumPy Basics3

Linear Algebra

There is a function dot, both an array method, and a function in the numpy namespace, for matrix multiplication:

x = np.array([[1, 2, 3], [4, 5, 6]])  
y = np.array([[6, 23], [-1, 7], [8, 9]])  
x  
#Output:
#array([[1, 2, 3],
#       [4, 5, 6]])
y  
#Output:
#array([[ 6, 23],
#       [-1,  7],
#       [ 8,  9]])
x.dot(y)  
#Output:
#array([[ 28,  64],
#       [ 67, 181]])

A matrix product between a 2D array and a suitably sized 1D array results in a 1D array:

np.dot(x, np.ones(3))  
#Output:
#array([  6.,  15.])

numpy.linalg has a standard set of matrix decompositions and things like inverse and determinant.

from numpy.linalg import inv, qr  
X = randn(5, 5) #A 5*5 matrix  
mat = X.T.dot(X) # matrix X's transposition multiple X  
inv(mat) #inverse of mat  
#Output:
#array([[   3.99634927,   -8.60568268,   10.9985889 ,    4.30635948,        -19.89346638],
#       [  -8.60568268,   19.44536324,  -24.1954652 ,   -9.50208014,         44.30059249],
#       [  10.9985889 ,  -24.1954652 ,   31.12382956,   11.90611367,        -56.04639441],
#       [   4.30635948,   -9.50208014,   11.90611367,    5.11939852,        -21.66116917],
#       [ -19.89346638,   44.30059249,  -56.04639441,  -21.66116917,       102.24533701]])

mat.dot(inv(mat)) #mat multiple inverse of mat(nearly Zero)  
Out[27]:  
#array([[  1.00000000e+00,-7.19381762e-15,1.21321962e-14,7.49195900e-16,-3.08386607e-14],
#       [  5.43036456e-15,1.00000000e+00,-2.07803635e-14,-2.60814354e-15,3.19665692e-14],
#       [ -5.29575083e-15,-1.71795908e-15,1.00000000e+00,-7.93173245e-15,-1.91802987e-15],
#       [ -1.75635977e-15,2.25691753e-15,-6.99940965e-15,1.00000000e+00,-6.98859926e-15],
#       [  9.65822515e-16,-1.97044265e-14,1.59646563e-14,2.23927654e-16,1.00000000e+00]])

More about numpy.linalg

Random Number Generation

The numpy.random module supplements the built-in Python random with functions for efficiently generating whole arrays of sample values from many kinds of probability distributions.

samples = np.random.normal(size=(4, 4))  
samples  
#Output:
#array([[-0.02778426, -2.30306873,  0.44936783,  0.02236305],
#       [ 0.09586369,  0.49977571,  0.17674787,  0.76304052],
#       [ 1.78797042, -0.32190291, -0.6569864 ,  0.63058886],
#       [-0.0686503 ,  0.30850549,  0.11803601, -1.13768578]])

numpy.random is well over an order of magnitude faster for generating very large samples:

from random import normalvariate  
N = 10000  
%timeit samples = [normalvariate(0, 1) for _ in xrange(N)]
#Output:100 loops, best of 3: 7.44 ms per loop
%timeit np.random.normal(size=N)
#Output:1000 loops, best of 3: 345 ┬Ás per loop

More about numpy.random

Example: Random Walks

Let's first consider a simple random walk starting at 0 with steps of 1 and -1 occurring with equal probability.

import random

position = 0  
walk = [position]  
steps = 1000  
for i in xrange(steps):  
    step = 1 if random.randint(0, 1) else -1
    position += step
    walk.append(position)

Thus, using the np.random module to draw 1000 coin flips at one, set these to 1 and -1, and compute the cumulative sum:

nsteps = 1000  
draws = np.random.randint(0, 2, size = nsteps)  
steps = np.where(draws > 0 , 1, -1)  
walk = steps.cumsum()  
walk.min()  
#Output: -49
walk.max()  
#Output: 2