In [1]:
from numpy import *

In [2]:
from mth309 import *


## Is det(AB) = det(A)det(B) a plausible conjecture?¶

In [3]:
n = 9
A = Matrix(random.randint(-9,10,(9,9)))
A

Out[3]:
 -1  2 -4 -9  5  0 -3  4  7
-9  5  3  1  3 -2  4  6  0
6 -2  1  0 -8  4 -1 -5  4
0 -2  4 -5 -5  0 -4 -6  9
9 -5  8  1 -9  7  3  9 -1
-8 -9  9 -4  1 -4 -5  5  2
1  1 -5 -7 -5  4 -9 -8  7
4  7  0 -9  9 -2  7 -1  7
-6 -7 -7 -3  8  5  8  1 -8

In [4]:
B = Matrix(random.randint(-9,10,(9,9)))
B

Out[4]:
 -1 -8 -6 -5  7 -4 -5 -2  5
-3  9  2  7 -4 -4  9 -5 -4
4  4  0 -7 -4 -5 -7 -3  3
2  3  3  1 -8 -8  8  4  5
-5  4 -4 -9 -6 -5  9  7 -5
6 -4 -5  3 -7 -2 -5 -4 -7
-3 -4 -9 -9 -5 -7 -1 -5 -4
-8 -2 -5  4  4  6  1 -6 -4
-9  9  7  1 -3 -2  6  0  0

In [5]:
AB = Matrix(dot(A,B))
AB

Out[5]:
 -150   70   19   43   53  94   73   -6  -99
-79  124   -1   15 -103 -10  116  -39  -92
75  -60   34   26   39 -20 -127  -42   77
16   72  130   19   31  27  -79   -1   60
55 -185 -139  -25   72   2 -254 -173   42
-9   55   68  -44   37 119  -85   36   26
39   39  145  143   73  91  -39   22   16
-176   85  -72 -127  -28 -82   96  -37 -108
23 -146 -177 -106  -57  27  -16   46 -149

In [6]:
det(A)

Out[6]:
30350008
In [7]:
det(B)

Out[7]:
-682607336
In [8]:
det(AB)

Out[8]:
-20717138108458688
In [9]:
det(A)*det(B)

Out[9]:
-20717138108458688

On the basis of this random test, yes! It is plausible (though not proved by one example!).

## Misled by uniform sampling¶

In [11]:
from mystery import f

In [12]:
f(.03)

Out[12]:
2.2295279470212641

Estimate the average of f over [0,1]. Sample at 17 uniformly spaced points in [0,1] and average.

In [16]:
x = linspace(0,1,17)
y = f(x)
print(y)
y.mean()

[ 1.91012145  1.877126    1.84426586  1.81157721  1.77909604  1.74685812
1.71489896  1.68325373  1.6519573   1.62104412  1.59054823  1.56050321
1.53094215  1.5018976   1.47340155  1.44548536  1.41817979]

Out[16]:
1.6565386296040427

True value is 2. Why are we so far off?

In [17]:
import matplotlib.pyplot as plt
%matplotlib inline

In [20]:
plt.plot(x,y,'ro');


Looks linearish decreasing. Look at a finer sample mesh ...

In [22]:
xx = linspace(0,1,1000)
yy = f(xx)
plt.plot(xx,yy)
plt.plot(x,y,'ro')

Out[22]:
[<matplotlib.lines.Line2D at 0x7f6479c8f748>]

Random sampling as an alternative ...

In [23]:
x = random.rand(17)
y = f(x)
#print(y)
y.mean()

Out[23]:
2.1063333401531952
In [24]:
ntrials = 20
for i in range(ntrials):
x = random.rand(17)
y = f(x)
#print(y)
print(y.mean())

2.06073520275
2.01139293142
2.10673784302
2.01237498507
1.83447307077
2.07979189399
1.89866626868
2.3379136183
2.23098457839
1.98148662999
1.89484059003
1.98688403936
1.97742268344
1.70987173803
1.88713209721
2.24816574599
2.10924688215
2.08831564206
1.85461137609
2.28262063841


## Source of "random" numbers?¶

In [26]:
# recall mayfly model
b = 4
def g(r): return b*(1-r)*r

In [27]:
r = .1
for i in range(20):
print(r)
r = g(r)

0.1
0.36000000000000004
0.9216
0.28901376000000006
0.8219392261226498
0.5854205387341974
0.970813326249438
0.11333924730376121
0.4019738492975123
0.9615634951138128
0.1478365599132853
0.5039236458651636
0.9999384200124991
0.0002463047816241282
0.0009849764623147091
0.00393602513473358
0.015682131363489303
0.061744808477550275
0.23172954841448365
0.7121238592244125


This is irregular (technically "chaotic"), non-periodic - it is sort-of random.

But .. no good for our purposes here. (i) very non-uniform

In [28]:
rs = [r]
for k in range(10000):
r = g(r)
rs.append(r)

In [29]:
rs[:10]

Out[29]:
[0.8200138733909665,
0.5903644833492417,
0.9673370405960985,
0.1263843619475224,
0.4416454200105602,
0.9863789719770243,
0.05374198247429209,
0.2034151271760999,
0.6481496528481241,
0.9122067214439212]
In [31]:
from myhistogram import myhistogram
myhistogram(rs,0,1,50)

Out[31]:
(array([ 0.02,  0.04,  0.06,  0.08,  0.1 ,  0.12,  0.14,  0.16,  0.18,
0.2 ,  0.22,  0.24,  0.26,  0.28,  0.3 ,  0.32,  0.34,  0.36,
0.38,  0.4 ,  0.42,  0.44,  0.46,  0.48,  0.5 ,  0.52,  0.54,
0.56,  0.58,  0.6 ,  0.62,  0.64,  0.66,  0.68,  0.7 ,  0.72,
0.74,  0.76,  0.78,  0.8 ,  0.82,  0.84,  0.86,  0.88,  0.9 ,
0.92,  0.94,  0.96,  0.98,  1.  ]),
array([ 4.48455154,  1.96980302,  1.56484352,  1.16988301,  1.04489551,
1.06489351,  1.07989201,  0.88491151,  0.81491851,  0.88991101,
0.86491351,  0.71492851,  0.64993501,  0.69993001,  0.70492951,
0.64493551,  0.70992901,  0.63493651,  0.62993701,  0.65493451,
0.71992801,  0.65493451,  0.76492351,  0.66493351,  0.70992901,
0.58994101,  0.55494451,  0.59994001,  0.62493751,  0.66993301,
0.63993601,  0.65493451,  0.71992801,  0.61493851,  0.75492451,
0.61493851,  0.68493151,  0.69993001,  0.66993301,  0.75492451,
0.74992501,  0.79992001,  0.84991501,  0.98490151,  1.01489851,
1.13488651,  1.19988001,  1.46985301,  1.97480252,  4.57954205]))
In [34]:
len(rs)
rs = rs[:-1]
len(rs)

Out[34]:
10000

and successive pairs confined to the graph of g:

In [35]:
x = rs[0::2]
y = rs[1::2]
plt.plot(x,y,'b.')

Out[35]:
[<matplotlib.lines.Line2D at 0x7f6479a6e4a8>]

## Linear Congruential Generators¶

In [47]:
# my LCG
# Linear congruential pseudo-random number generator

rng_a = 427419669081    # multiplier
rng_m = 999999999989    # mod

def mylcg():
global rng_a, rng_m, rng_i
rng_i = (rng_a*rng_i)%rng_m
return rng_i/rng_m

rng_i = 17             # seed

rs = []
n = 10000
for k in range(n): rs.append(mylcg())

#rs


## Test for uniformity¶

In [44]:
myhistogram(rs,0,1,100);


## test for "filling the dartboard"¶

In [49]:
x = rs[0::2]
y = rs[1::2]
plt.subplot(111,aspect=1)
plt.plot(x,y,'b.',alpha=0.1);


looks uniform!

## Randu - very bad, do NOT use¶

In [50]:
# my LCG
# Linear congruential pseudo-random number generator

rng_i = 17          # seed
rng_a = 2**16+3         # multiplier
rng_m = 2**31           # mod

def mylcg():
global rng_a, rng_m, rng_i
rng_i = (rng_a*rng_i)%rng_m
return rng_i/rng_m

rng_i = 17             # seed

rs = []
n = 10000
for k in range(n): rs.append(mylcg())

#rs

In [51]:
x = rs[0::2]
y = rs[1::2]
plt.subplot(111,aspect=1)
plt.plot(x,y,'b.',alpha=0.1);


Looks ok here, but suppose we use it to generate points in 3D ...

# Monte Carlo integration¶

In [70]:
t = linspace(0,2*pi,1000)
r = sqrt(2 + cos(7*t))
x = r*cos(t)
y = r*sin(t)
R = sqrt(3)
plt.subplot(111,aspect=1)
plt.fill(x,y,'#bb50aa')
rs = []
n = 20
rng_i = 194             # seed
for k in range(n): rs.append(mylcg())
rs = array(rs)
rx = 2*sqrt(3)*rs[0::2]-sqrt(3)
ry = 2*sqrt(3)*rs[1::2]-sqrt(3)
plt.plot(rx,ry,'b.',alpha=1)
plt.xlim(-R,R)
plt.ylim(-R,R);

In [74]:
rr = sqrt(rx**2 + ry**2)
rr
rt = arctan2(ry,rx)
rt
rf = sqrt(2+cos(7*rt))
rr<rf
(rr<rf).sum()/(n/2)
# n/2 because we used n random numbers to generate n/2 points

Out[74]:
0.80000000000000004