Instantly share code, notes, and snippets.

# ceteri/2.01.pySecret Last active Aug 29, 2015

Just Enough Math - example code + data + results
 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])
 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
 >>> 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}
 ## 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)
 >>> for bowl in ["H1", "H2"]: ... print bowl, pmf.Prob(bowl) ... H1 0.6 H2 0.4
 ## 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)
 >>> 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
 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)
 >>> 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]
 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)
 >>> 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]] >>>
 >>> import numpy as np >>> A = np.array([1, 2, 3, 4]).reshape(2, 2) >>> >>> print np.linalg.det(A) -2.0
 >>> 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]
 >>> 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]
 >>> 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]]))
 ## @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
 >>> 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]
 >>> 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]]))
 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)
 \$ 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
 ## 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)
 >>> 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) >>>
 /* 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;
 \$ 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
 ## 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()
 \$ 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]]