Just Enough Math - example code + data + results

2.01.py

class Monoid (object):
    def __init__ (self, null, lift, op):
        self.null = null
        self.lift = lift
        self.op   = op
 
    def fold (self, xs):
        if hasattr(xs, "__fold__"):
            return xs.__fold__(self)
        else:
            return reduce(self.op, (self.lift(x) for x in xs), self.null)
 
    def __call__ (self, *args):
        return self.fold(args)
 
    def star (self):
        return Monoid(self.null, self.fold, self.op)


summ   = Monoid(0,  lambda x: x,      lambda a,b: a+b)
joinm  = Monoid('', lambda x: str(x), lambda a,b: a+b)
listm  = Monoid([], lambda x: [x],    lambda a,b: a+b)
tuplem = Monoid((), lambda x: (x,),   lambda a,b: a+b)
lenm   = Monoid(0,  lambda x: 1,      lambda a,b: a+b)
prodm  = Monoid(1,  lambda x: x,      lambda a,b: a*b)


## extended to define a monoid for folding Python `dict`

def dict_op (a, b):
    for key, val in b.items():
        if not key in a:
            a[key] = val
        else:
            a[key] += val

    return a


dictm  = Monoid({}, lambda x: x,      lambda a,b: dict_op(a, b))


if __name__=='__main__':
    x1 = { "a": 2, "b": 3 }
    x2 = { "b": 2, "c": 7 }

    print x1, x2
    print dictm.fold([x1, x2])

2.02.py

summ.star().fold([[10, 20, 30], [40, 50]])
    == summ.fold([10, 20, 30]) + summ.fold([40, 50])
    == (10 + 20 + 30) + (40 + 50)
    == 10 + 20 + 30 + 40 + 50

2.03.py

>>> listm.star().fold([[1, 2, 3], [4, 5, 6]])
[1, 2, 3, 4, 5, 6]
>>> summ.star().fold([[10, 20, 30], [40, 50]])
150

>>> print dictm.fold([x1, x2])
{'a': 4, 'c': 14, 'b': 10}

3.01.py

## based on "Cookie" problem                                                              
## http://www.greenteapress.com/thinkbayes/                                               

from thinkbayes import Pmf

pmf = Pmf()
pmf.Set("H1", 1 / 2.0)
pmf.Set("H2", 1 / 2.0)

pmf.Mult("H1", 30 / 40.0)
pmf.Mult("H2", 20 / 40.0)

pmf.Normalize()

for bowl in ["H1", "H2"]:
    print bowl, pmf.Prob(bowl)

3.02.py

>>> for bowl in ["H1", "H2"]:
...     print bowl, pmf.Prob(bowl)
... 
H1 0.6
H2 0.4

4.01.py

## use NumPy to calculate vector operations

import numpy as np

a = 3
b = 2
x = np.array([ 2, 3, 1, 0 ])
y = np.array([ 5, 1, 5, 0 ])

print np.add(x, a)

print np.dot(b, x)

print np.add(x, y)

print np.dot(x, y)

4.02.py

>>> import numpy as np
>>> 
>>> a = 3
>>> b = 2
>>> x = np.array([ 2, 3, 1, 0 ])
>>> y = np.array([ 5, 1, 5, 0 ])
>>> 
>>> print np.add(x, a)
[5 6 4 3]
>>> 
>>> print np.dot(b, x)
[4 6 2 0]
>>> 
>>> print np.add(x, y)
[7 4 6 0]
>>> 
>>> print np.dot(x, y)
18

4.03.py

import numpy as np

A = np.array([1, 2, 3, 4]).reshape(2, 2)
B = np.array([5, 6, 7, 8]).reshape(2, 2)
x = [2, 3]

print A
print B

print np.add(A, B)
print np.dot(A, x)

4.04.py

>>> import numpy as np
>>> 
>>> A = np.array([1, 2, 3, 4]).reshape(2, 2)
>>> B = np.array([5, 6, 7, 8]).reshape(2, 2)
>>> x = [2, 3]
>>> 
>>> print A
[[1 2]
 [3 4]]
>>> print B
[[5 6]
 [7 8]]
>>> 
>>> print np.add(A, B)
[[ 6  8]
 [10 12]]
>>> print np.dot(A, x)
[ 8 18]

4.05.py

import numpy as np
A = np.array([1, 2, 3, 4]).reshape(2, 2)
I = np.array([1, 0, 0, 1]).reshape(2, 2)

print I
print np.dot(A, I)

4.06.py

>>> import numpy as np
>>> A = np.array([1, 2, 3, 4]).reshape(2, 2)
>>> I = np.array([1, 0, 0, 1]).reshape(2, 2)
>>> 
>>> print I
[[1 0]
 [0 1]]
>>> print np.dot(A, I)
[[1 2]
 [3 4]]
>>>

4.07.py

>>> import numpy as np
>>> A = np.array([1, 2, 3, 4]).reshape(2, 2)
>>> 
>>> print np.linalg.det(A)
-2.0

4.08.py

>>> import numpy as np
>>> 
>>> A = np.array([3, 9, 4, 8]).reshape(2, 2)
>>> y = [ 5, 12 ]
>>> invA = np.linalg.inv(A)
>>> 
>>> print np.linalg.det(A)
-12.0
>>> print invA
[[-0.66666667  0.75      ]
 [ 0.33333333 -0.25      ]]
>>> print np.dot(invA, y)
[ 5.66666667 -1.33333333]

4.09.py

>>> import numpy as np
>>> 
>>> x = np.array([0, 1, 2, 3])
>>> y = np.array([-1.0, 0.2, 0.9, 2.1])
>>> A = np.vstack([x, np.ones(len(x))]).T
>>> 
>>> (m, b), resid, rank, singular = np.linalg.lstsq(A, y)
>>> 
>>> print A
[[ 0.  1.]
 [ 1.  1.]
 [ 2.  1.]
 [ 3.  1.]]
>>> print (m, b), resid
(1.0, -0.95) [ 0.05]

4.10.py

>>> import numpy as np
>>> 
>>> A = np.array([0, 1, -2, -3]).reshape(2, 2)
>>> 
>>> print np.linalg.eig(A)
(array([-1., -2.]), array([[ 0.70710678, -0.4472136 ],
       [-0.70710678,  0.89442719]]))

4.11.py

## @JustGlowing 2011-07-22
## http://glowingpython.blogspot.com/2011/07/principal-component-analysis-with-numpy.html

import numpy as np

def princomp (A):
    """
    principal components analysis on the N x P data matrix A
    rows correspond to observations, columns to variables
    """

    # subtract the mean (along columns)
    M = (A - np.mean(A.T, axis=1)).T

    # NB: results are not always sorted
    [latent, coeff] = np.linalg.eig(np.cov(M))

    # projection of the data in the new space
    score = np.dot(coeff.T, M)

    return coeff, score, latent


## analyze each judge's scores for our bartendrs
A = np.array([
        [ 5.1, 4.9, 4.7, 4.6, 5.7, 5.7, 6.2, 5.1, 5.7, 6.5, 7.7, 7.7, 6.0, 6.9, 5.6 ],
        [ 3.5, 3.0, 3.2, 3.1, 3.0, 2.9, 2.9, 2.5, 2.8, 3.0, 3.8, 2.6, 2.2, 3.2, 2.8 ],
        [ 1.4, 1.4, 1.3, 1.5, 4.2, 4.2, 4.3, 3.0, 4.1, 5.5, 6.7, 6.9, 5.0, 5.7, 4.9 ],
        [ 0.2, 0.2, 0.2, 0.2, 1.2, 1.3, 1.3, 1.1, 1.3, 1.8, 2.2, 2.3, 1.5, 2.3, 2.0 ],
        ])

coeff, score, latent = princomp(A.T)

print coeff
print score
print latent

portion_variance = map(lambda x: x / sum(latent), latent)
print portion_variance

4.12.py

>>> A = np.array([
...         [ 5.1, 4.9, 4.7, 4.6, 5.7, 5.7, 6.2, 5.1, 5.7, 6.5, 7.7, 7.7, 6.0, 6.9, 5.6 ],
...         [ 3.5, 3.0, 3.2, 3.1, 3.0, 2.9, 2.9, 2.5, 2.8, 3.0, 3.8, 2.6, 2.2, 3.2, 2.8 ],
...         [ 1.4, 1.4, 1.3, 1.5, 4.2, 4.2, 4.3, 3.0, 4.1, 5.5, 6.7, 6.9, 5.0, 5.7, 4.9 ],
...         [ 0.2, 0.2, 0.2, 0.2, 1.2, 1.3, 1.3, 1.1, 1.3, 1.8, 2.2, 2.3, 1.5, 2.3, 2.0 ],
...         ])
>>> 
>>> coeff, score, latent = princomp(A.T)
>>> 
>>> print coeff
[[ 0.41517419  0.57141489  0.62117063  0.33950326]
 [-0.01541612  0.77708865 -0.59860493 -0.19382391]
 [ 0.84516352 -0.19642892 -0.12061365 -0.48225165]
 [ 0.33629059 -0.1761645  -0.49119399  0.78396631]]
>>> print score
[[ -2.89330144e+00  -2.96862822e+00  -3.13926264e+00  -3.01020574e+00
    6.62595725e-02   1.01430244e-01   3.93533691e-01  -1.22296216e+00
    1.84555042e-02   1.69888585e+00   3.33347445e+00   3.55463555e+00
    9.80162718e-01   2.19905031e+00   8.88472315e-01]
 [  6.73661053e-01   1.70833750e-01   2.31611394e-01   5.74752546e-02
   -9.81998222e-02  -1.93525137e-01   7.25394171e-02  -5.76261923e-01
   -2.51591109e-01  -2.12420839e-03   9.99064073e-01   9.65546015e-03
   -7.58438762e-01   2.54491445e-01  -5.89190886e-01]
 [  4.19865554e-02   2.17054897e-01  -1.48388497e-02  -4.12181496e-02
   -1.14920820e-01  -1.04179726e-01   1.94344222e-01   5.53505357e-03
   -3.22578670e-02  -6.94984616e-02  -1.44191636e-01   5.00892154e-01
    3.06465195e-01  -2.10470923e-01  -5.34701644e-01]
 [  4.96901878e-02   7.87014923e-02   2.02612221e-02  -9.07570427e-02
   -2.16034211e-01  -1.18255189e-01   3.27127732e-03   1.77481138e-01
   -5.06476323e-02  -1.00978962e-01  -1.13749635e-01   1.00785362e-01
   -1.09735529e-01   2.91590384e-01   7.83771368e-02]]
>>> print latent
[ 5.07038591  0.21236083  0.05816678  0.01775315]
>>> portion_variance = map(lambda x: x / sum(latent), latent)
>>> print portion_variance
[0.94620289391157297, 0.039629415391115773, 0.010854711199668675, 0.0033129794976426439]

4.13.py

>>> import numpy as np
>>> 
>>> A = np.array([0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1]).reshape(4, 4)
>>> 
>>> print A
[[0 1 0 1]
 [1 0 1 1]
 [0 1 0 1]
 [1 1 0 1]]
>>> print np.linalg.eig(A)
(array([ 2.65109341, -1.37720285, -0.27389055,  0.        ]),
 array([[ -4.25648128e-01,  -4.56381963e-01,  -2.61003301e-01,
          4.53246652e-17],
       [ -5.42229129e-01,   7.53529250e-01,  -6.20457793e-01,
         -5.77350269e-01],
       [ -4.25648128e-01,  -4.56381963e-01,  -2.61003301e-01,
         -5.77350269e-01],
       [ -5.86203816e-01,  -1.24998709e-01,   6.91944132e-01,
          5.77350269e-01]]))

5.01.py

from glpk import LPX

# set up to maximize the objective function
lp = LPX()
lp.name = 'example'
lp.obj.maximize = True

# append 3 rows, named p, q, r
row_names = ["p", "q", "r"]
lp.rows.add(len(row_names))

for r in lp.rows:
    r.name = row_names[r.index]

lp.rows[0].bounds = None, 100.0
lp.rows[1].bounds = None, 600.0
lp.rows[2].bounds = None, 150.0

# append 3 cols, named x0, x1, x2
lp.cols.add(3)

for c in lp.cols:
    c.name = 'x%d' % c.index
    c.bounds = 0.0, None

# set the objective coefficients and
# non-zero entries of the constraint matrix
lp.obj[:] = [ 5.0, 3.0, 2.0 ]
lp.matrix = [ 1.0, 1.0, 3.0, 10.0, 4.0, 5.0, 1.0, 1.0, 3.0 ]

# report the objective function value and structural variables
lp.simplex(msg_lev=LPX.MSG_ALL)
print 'Z = %g;' % lp.obj.value,
print '; '.join('%s = %g' % (c.name, c.primal) for c in lp.cols)

5.02.py

$ python ch5_lp.py
*     0: obj =   0.000000000e+00  infeas =  0.000e+00 (0)
*     2: obj =   3.666666667e+02  infeas =  0.000e+00 (0)
OPTIMAL SOLUTION FOUND
Z = 366.667; x0 = 33.3333; x1 = 66.6667; x2 = 0

5.03.py

## based on an analytic solution, we expect that the local minimum of
## the polynomial function f(x) occurs at x = 9/4 = 2.25
 
def f (x):
    return x**4 - 3 * x**3 + 2

def f_prime (x):
    """first derivative of f(x)"""
    return 4 * x**3 - 9 * x**2

def g (x):
    return x**2

def g_prime (x):
    """first dervative of g(x)"""
    return 2 * x


## inputs: f, starting value x1, termination tolerances

def gradient_descent (x1, toler, step_size, max_iter, func, func_deriv):
    x0 = 0
    print "\t".join([ "i", "x1", "error", "func(x1)" ])

    for i in xrange(max_iter):
        error = abs(x1 - x0)
        print "\t".join([ "%d %0.4e %0.4e %3.4e" % (i, x1, error, func(x1)) ])

        if abs(x1 - x0) <= toler:
            return "local minimum:", x1

        x0 = x1
        x1 = x0 - step_size * func_deriv(x0)

    return "halted"


## NB: step size is set here, and is not adaptive

x1 = 6
toler = 0.00001
step_size = 0.01
max_iter = 93

print gradient_descent(x1, toler, step_size, max_iter, f, f_prime)

5.04.py

>>> x1 = 6
>>> toler = 0.00001
>>> step_size = 0.01
>>> max_iter = 93
>>> print gradient_descent(x1, toler, step_size, max_iter, f, f_prime)
i	x1	error	func(x1)
0 6.0000e+00 6.0000e+00 6.5000e+02
1 6.0000e-01 5.4000e+00 1.4816e+00
2 6.2376e-01 2.3760e-02 1.4233e+00
3 6.4907e-01 2.5309e-02 1.3571e+00
4 6.7605e-01 2.6978e-02 1.2819e+00
5 7.0482e-01 2.8774e-02 1.1964e+00
6 7.3553e-01 3.0704e-02 1.0989e+00
7 7.6830e-01 3.2773e-02 9.8789e-01
8 8.0328e-01 3.4985e-02 8.6137e-01
9 8.4062e-01 3.7341e-02 7.1727e-01
10 8.8046e-01 3.9837e-02 5.5332e-01
11 9.2293e-01 4.2467e-02 3.6711e-01
12 9.6815e-01 4.5216e-02 1.5620e-01
13 1.0162e+00 4.8060e-02 -8.1809e-02
14 1.0672e+00 5.0964e-02 -3.4906e-01
15 1.1211e+00 5.3883e-02 -6.4723e-01
16 1.1778e+00 5.6753e-02 -9.7725e-01
17 1.2373e+00 5.9495e-02 -1.3389e+00
18 1.2993e+00 6.2014e-02 -1.7305e+00
19 1.3635e+00 6.4199e-02 -2.1485e+00
20 1.4294e+00 6.5925e-02 -2.5872e+00
21 1.4965e+00 6.7066e-02 -3.0389e+00
22 1.5640e+00 6.7499e-02 -3.4937e+00
23 1.6311e+00 6.7121e-02 -3.9405e+00
24 1.6970e+00 6.5862e-02 -4.3677e+00
25 1.7607e+00 6.3702e-02 -4.7644e+00
26 1.8214e+00 6.0675e-02 -5.1215e+00
27 1.8782e+00 5.6878e-02 -5.4328e+00
28 1.9307e+00 5.2460e-02 -5.6956e+00
29 1.9783e+00 4.7609e-02 -5.9105e+00
30 2.0208e+00 4.2533e-02 -6.0807e+00
31 2.0583e+00 3.7433e-02 -6.2117e+00
32 2.0908e+00 3.2490e-02 -6.3098e+00
33 2.1186e+00 2.7843e-02 -6.3815e+00
34 2.1422e+00 2.3590e-02 -6.4327e+00
35 2.1620e+00 1.9788e-02 -6.4686e+00
36 2.1784e+00 1.6456e-02 -6.4933e+00
37 2.1920e+00 1.3584e-02 -6.5101e+00
38 2.2032e+00 1.1143e-02 -6.5214e+00
39 2.2123e+00 9.0928e-03 -6.5289e+00
40 2.2196e+00 7.3880e-03 -6.5338e+00
41 2.2256e+00 5.9814e-03 -6.5370e+00
42 2.2305e+00 4.8286e-03 -6.5391e+00
43 2.2343e+00 3.8887e-03 -6.5405e+00
44 2.2375e+00 3.1257e-03 -6.5414e+00
45 2.2400e+00 2.5085e-03 -6.5420e+00
46 2.2420e+00 2.0107e-03 -6.5423e+00
47 2.2436e+00 1.6100e-03 -6.5426e+00
48 2.2449e+00 1.2882e-03 -6.5427e+00
49 2.2459e+00 1.0300e-03 -6.5428e+00
50 2.2467e+00 8.2311e-04 -6.5429e+00
51 2.2474e+00 6.5751e-04 -6.5429e+00
52 2.2479e+00 5.2506e-04 -6.5429e+00
53 2.2483e+00 4.1918e-04 -6.5429e+00
54 2.2487e+00 3.3457e-04 -6.5430e+00
55 2.2489e+00 2.6700e-04 -6.5430e+00
56 2.2492e+00 2.1305e-04 -6.5430e+00
57 2.2493e+00 1.6998e-04 -6.5430e+00
58 2.2495e+00 1.3560e-04 -6.5430e+00
59 2.2496e+00 1.0817e-04 -6.5430e+00
60 2.2497e+00 8.6287e-05 -6.5430e+00
61 2.2497e+00 6.8826e-05 -6.5430e+00
62 2.2498e+00 5.4896e-05 -6.5430e+00
63 2.2498e+00 4.3785e-05 -6.5430e+00
64 2.2499e+00 3.4921e-05 -6.5430e+00
65 2.2499e+00 2.7852e-05 -6.5430e+00
66 2.2499e+00 2.2213e-05 -6.5430e+00
67 2.2499e+00 1.7716e-05 -6.5430e+00
68 2.2499e+00 1.4129e-05 -6.5430e+00
69 2.2500e+00 1.1268e-05 -6.5430e+00
70 2.2500e+00 8.9864e-06 -6.5430e+00
('local minimum:', 2.2499646074278457)
>>> 

5.05.py

/* foobartendr.io  pricing.mod */

var x0 >= 0;
var x1 >= 0;
var x2 >= 0;

maximize
value: 5.0 * x0 + 3.0 * x1 + 2.0 * x2;

subject to

ingrednt:  1.0 * x0 + 1.0 * x1 + 3.0 * x2 <= 100;
bartendr: 10.0 * x0 + 4.0 * x1 + 5.0 * x2 <= 600;
delivery:  1.0 * x0 + 1.0 * x1 + 3.0 * x2 <= 150;

end;

5.06.py

$ glpsol --math pricing.mod 
Reading model section from pricing.mod...
16 lines were read
Generating value...
Generating ingrednt...
Generating bartendr...
Generating delivery...
Model has been successfully generated
glp_simplex: original LP has 4 rows, 3 columns, 12 non-zeros
glp_simplex: presolved LP has 3 rows, 3 columns, 9 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  1.000e+01  ratio =  1.000e+01
Problem data seem to be well scaled
Crashing...
Size of triangular part = 3
*     0: obj =   0.000000000e+00  infeas =  0.000e+00 (0)
*     2: obj =   3.666666667e+02  infeas =  0.000e+00 (0)
OPTIMAL SOLUTION FOUND
Time used:   0.0 secs

5.07.py

## http://pythonhosted.org//neurolab/lib.html
import numpy as np
import neurolab as nl

# create a neural network with 3 inputs, 4 neurons
# in the hidden layer, and 1 in the output layer
net = nl.net.newff([[0, 1],[0, 1],[0, 1]], [4, 1])
input = [[0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0]]
target = [[0], [0], [1], [1]]

# train
net.trainf = nl.train.train_gdx
error = net.train(input, target, show=25)

# test
print "expect 0", net.sim([[0, 1, 0]])
print "expect 1", net.sim([[1, 1, 0]])

# visualize
import pylab as pl
pl.plot(error)
pl.xlabel('Epoch number')
pl.ylabel('Train error')
pl.grid()
pl.show()

5.08.py

$ python ch5_nn.py
Epoch: 25; Error: 0.714363135207;
Epoch: 50; Error: 0.6290077045;
Epoch: 75; Error: 0.481297387105;
Epoch: 100; Error: 0.382371205263;
Epoch: 125; Error: 0.32463469683;
Epoch: 150; Error: 0.0345094306294;
The goal of learning is reached
expect 0 [[-0.01642426]]
expect 1 [[ 0.90354536]]