"This work was done by Adrian Cowan"
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{This work was done by Adrian Cowan}
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{This work was done by Adrian Cowan}
|
"Part A"
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part A}
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part A}
|
forget()
y = var('y');
o_X = var('sigma_X');
o_Y = var('sigma_Y');
p = var('rho');
u_X = var('mu_X');
u_Y = var('mu_Y');
|
|
f(x, y, o_X, o_Y, p, u_X, u_Y)=(1/(2*pi*o_X*o_Y*sqrt(1-p^2)))*e^((-1/(2*(1-p^2)))*(((x-u_X)/o_X)^2-2*p*((x-u_X)/o_X)*((y-u_Y)/o_Y)+((y-u_Y)/o_Y)^2));f
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{e^{\left(-\frac{\frac{2 \, {\left(u_{Y} - y\right)} {\left(u_{X} - x\right)} p}{o_{X} o_{Y}} - \frac{{\left(u_{Y} - y\right)}^{2}}{o_{Y}^{2}} - \frac{{\left(u_{X} - x\right)}^{2}}{o_{X}^{2}}}{2 \, {\left(p^{2} - 1\right)}}\right)}}{2 \, \sqrt{-p^{2} + 1} \pi o_{X} o_{Y}}
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{e^{\left(-\frac{\frac{2 \, {\left(u_{Y} - y\right)} {\left(u_{X} - x\right)} p}{o_{X} o_{Y}} - \frac{{\left(u_{Y} - y\right)}^{2}}{o_{Y}^{2}} - \frac{{\left(u_{X} - x\right)}^{2}}{o_{X}^{2}}}{2 \, {\left(p^{2} - 1\right)}}\right)}}{2 \, \sqrt{-p^{2} + 1} \pi o_{X} o_{Y}}
|
plot3d(f(x,y,1,1,.5,0,0),(x,-5,5),(y,-5,5))
|
|
"Part B"
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part B}
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part B}
|
assume(-1<=p<=1)
assume(o_X>0)
assume(o_Y>0)
part = integral((x-u_X)*(y-u_Y)*f,(x, -oo, oo), algorithm='maxima');part
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{\sqrt{-2 \, p^{2} + 2} {\left(u_{Y} - y\right)} {\left(o_{X} p u_{Y} - o_{X} p y\right)} e^{\left(-\frac{u_{Y}^{2} - 2 \, u_{Y} y + y^{2}}{2 \, o_{Y}^{2}}\right)}}{2 \, \sqrt{-p^{2} + 1} \sqrt{\pi} o_{Y}^{2}}
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{\sqrt{-2 \, p^{2} + 2} {\left(u_{Y} - y\right)} {\left(o_{X} p u_{Y} - o_{X} p y\right)} e^{\left(-\frac{u_{Y}^{2} - 2 \, u_{Y} y + y^{2}}{2 \, o_{Y}^{2}}\right)}}{2 \, \sqrt{-p^{2} + 1} \sqrt{\pi} o_{Y}^{2}}
|
cov = integral(part, (y, -oo, oo)).simplify_full();cov
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ o_{X} o_{Y} p
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ o_{X} o_{Y} p
|
cor = cov/sqrt(o_X^2*o_Y^2);cor
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{o_{X} o_{Y} p}{\sqrt{o_{X}^{2} o_{Y}^{2}}}
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ \frac{o_{X} o_{Y} p}{\sqrt{o_{X}^{2} o_{Y}^{2}}}
|
cor.simplify_full()
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ p
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, y, o_{X}, o_{Y}, p, u_{X}, u_{Y} \right) \ {\mapsto} \ p
|
"Part C"
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part C}
\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{Part C}
|
print "Rayleigh distribution"
ray_pdf(x, sigma) = (x/(sigma^2))*e^(-x^2/(2*sigma^2));ray_pdf
Rayleigh distribution
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, \sigma \right) \ {\mapsto} \ \frac{x e^{\left(-\frac{x^{2}}{2 \, \sigma^{2}}\right)}}{\sigma^{2}}
Rayleigh distribution
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, \sigma \right) \ {\mapsto} \ \frac{x e^{\left(-\frac{x^{2}}{2 \, \sigma^{2}}\right)}}{\sigma^{2}}
|
print "Rayleigh cumulative distribution"
o = var('sigma')
ray_cdf(x, o) = 1-e^(-x^2/(2*o^2));ray_cdf
Rayleigh cumulative distribution
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, o \right) \ {\mapsto} \ -e^{\left(-\frac{x^{2}}{2 \, o^{2}}\right)} + 1
Rayleigh cumulative distribution
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, o \right) \ {\mapsto} \ -e^{\left(-\frac{x^{2}}{2 \, o^{2}}\right)} + 1
|
ray(x,o) = o*sqrt(-2*ln(x));ray
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, o \right) \ {\mapsto} \ o \sqrt{-\log\left(x\right)} \sqrt{2}
\newcommand{\Bold}[1]{\mathbf{#1}}\left( x, o \right) \ {\mapsto} \ o \sqrt{-\log\left(x\right)} \sqrt{2}
|
plot(ray(x,1),(x,0,1))
|
|
ray(0,1)
\newcommand{\Bold}[1]{\mathbf{#1}}+\infty
\newcommand{\Bold}[1]{\mathbf{#1}}+\infty
|
ray(1,1)
\newcommand{\Bold}[1]{\mathbf{#1}}0
\newcommand{\Bold}[1]{\mathbf{#1}}0
|
a = []
for i in (0..10000):
a.append([i,ray(random(),1)])
scatter_plot(a, markersize=1)
|
|
a = []
for i in (0..1000000):
a.append([ray(random(),1),random()*2*pi])
len(a)
Traceback (click to the left of this block for traceback) ... NameError: name 'ray' is not defined Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_4.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("YSA9IFtdCmZvciBpIGluICgwLi4xMDAwMDAwKToKICAgIGEuYXBwZW5kKFtyYXkocmFuZG9tKCksMSkscmFuZG9tKCkqMipwaV0pCmxlbihhKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
File "", line 1, in <module>
File "/tmp/tmpnEkF7T/___code___.py", line 5, in <module>
a.append([ray(random(),_sage_const_1 ),random()*_sage_const_2 *pi])
NameError: name 'ray' is not defined
|
b = []
for item in a:
b.append([item[0]*cos(item[1]),item[0]*sin(item[1])])
Traceback (click to the left of this block for traceback) ... NameError: name 'a' is not defined Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_3.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("YiA9IFtdCmZvciBpdGVtIGluIGE6CiAgICBiLmFwcGVuZChbaXRlbVswXSpjb3MoaXRlbVsxXSksaXRlbVswXSpzaW4oaXRlbVsxXSldKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
File "", line 1, in <module>
File "/tmp/tmptcBNWE/___code___.py", line 4, in <module>
exec compile(u'for item in a:\n b.append([item[_sage_const_0 ]*cos(item[_sage_const_1 ]),item[_sage_const_0 ]*sin(item[_sage_const_1 ])])
File "", line 1, in <module>
NameError: name 'a' is not defined
|
scatter_plot(b, markersize=2)
|
|
bin=.3
s=1/bin;index = []
for i in (-4*s..4*s):
index.append(i/s)
index
def bin(x):
return int((round(x*s,0))+4*s)
plot(bin, (x,-4,4))
|
|
h=[]
for i in (0 .. 8*s):
h.append([])
for j in (0 .. 8*s):
h[int(i)].append(0)
for item in b:
i = bin(item[0])
j = bin(item[1])
h[int(i)][int(j)] += 1
|
|
list_plot3d(h,interpolation_type='nn')
|
|
ray(x,o) = o*sqrt(-2*ln(x));ray
a = []
for i in (0..1000000):
a.append([ray(random(),1),random()*2*pi])
len(a)
b = []
for item in a:
b.append([item[0]*cos(item[1]),item[0]*sin(item[1])])
bin=.2
s=1/bin;index = []
for i in (-4*s..4*s):
index.append(i/s)
index
def bin(x):
return int((round(x*s,0))+4*s)
plot(bin, (x,-4,4))
h=[]
for i in (0 .. 8*s):
h.append([])
for j in (0 .. 8*s):
h[int(i)].append(0)
for item in b:
i = bin(item[0])
j = bin(item[1])
h[int(i)][int(j)] += 1
list_plot3d(h,interpolation_type='nn')
Traceback (click to the left of this block for traceback) ... __SAGE__ ^CTraceback (most recent call last): for item in a:
File "", line 1, in <module>
File "/tmp/tmpD8uQFS/___code___.py", line 6, in <module>
a.append([ray(random(),_sage_const_1 ),random()*_sage_const_2 *pi])
File "expression.pyx", line 3473, in sage.symbolic.expression.Expression.__call__ (sage/symbolic/expression.cpp:15647)
File "/home/sage/sage_install/sage-alpha/local/lib/python2.6/site-packages/sage/symbolic/callable.py", line 451, in _call_element_
return SR(_the_element.substitute(**d))
File "expression.pyx", line 3324, in sage.symbolic.expression.Expression.substitute (sage/symbolic/expression.cpp:15026)
File "pynac.pyx", line 606, in sage.symbolic.pynac.py_rational_power_parts (sage/symbolic/pynac.cpp:5818)
File "rational.pyx", line 270, in sage.rings.rational.rational_power_parts (sage/rings/rational.c:5586)
File "rational.pyx", line 280, in sage.rings.rational.rational_power_parts (sage/rings/rational.c:5837)
File "integer.pyx", line 3137, in sage.rings.integer.Integer.factor (sage/rings/integer.c:18989)
File "integer_ring.pyx", line 1071, in sage.rings.integer_ring.factor (sage/rings/integer_ring.c:9546)
File "/home/sage/sage_install/sage-alpha/local/lib/python2.6/site-packages/sage/rings/arith.py", line 2324, in factor
return factorization.Factorization(__factor_using_trial_division(n), unit)
File "/home/sage/sage_install/sage-alpha/local/lib/python2.6/site-packages/sage/rings/arith.py", line 2108, in __factor_using_trial_division
while n%p == 0:
KeyboardInterrupt
__SAGE__
|
list_plot3d(h,interpolation_type='nn')
Traceback (click to the left of this block for traceback) ... NameError: name 'h' is not defined Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_3.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("bGlzdF9wbG90M2QoaCxpbnRlcnBvbGF0aW9uX3R5cGU9J25uJyk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
File "", line 1, in <module>
File "/tmp/tmpSWUUfk/___code___.py", line 2, in <module>
exec compile(u"list_plot3d(h,interpolation_type='nn')" + '\n', '', 'single')
File "", line 1, in <module>
NameError: name 'h' is not defined
|
|
|