Parametrized Monte-Carlo¶
Nikolai Shokhirev¶
In [1]:
import numpy as np
\[ I(p)=\int_{0}^{\pi}cos(px)sin(x)dx = \frac{1+cos(\pi p)}{1 - p^2} \]
\[ I(p) \approx \frac{\pi^2 q}{4}\left(\frac{q}{2} - 1\right) , \; q = p - 1\]
error: \[ \epsilon = \frac{\pi^2 q^3}{16}\left(\frac{\pi^2}{3} - 1\right) \approx 1.412505788140299 \; q^3 \]
\[ \frac {\pi^2} {4} = 2.467401100272340 \]
In [2]:
def int_cs(p):
''' exact integral (for comparison) '''
q = p - 1
Pi = math.pi
if (abs(q) > 0.001):
return (1+cos(Pi*p))/(1-p*p)
else:
return 2.4674011*q*(0.5*q - 1)
Make int_cs work on vectors
In [3]:
vint_cs = np.vectorize(int_cs)
In [4]:
pp = np.linspace(0.0,3.0,151)
Exact integral generation
In [5]:
jj = vint_cs(pp)
plot(pp,jj)
title('Exact integral')
xlabel('$p$')
ylabel('$I(p)$')
Out[5]:
Monte Carlo integration \[ I(p)=\int_{a}^{b}f(x; p)dx \approx \frac{1}{N}\sum_{i=1}^{N}f(x_{i}; p) \]
In [6]:
def MC(p, f, a, b, n, reset = False):
s = 0.0
if reset:
np.random.seed(1)
x = np.random.uniform(a, b, n)
y = f(x, p)
return np.average(y)*(b-a)
In [7]:
def cs(x, p):
return cos(p*x)*sin(x)
Test
In [8]:
I = MC(0.0, cs, 0.0, math.pi, 1000)
print("I(0) = ", I, " - expected 2")
Generate \(I(p)\) without seed reset
In [9]:
def IMC(p, reset=False):
return MC(p, cs, 0.0, math.pi, 500, reset)
In [10]:
vIMC = np.vectorize(IMC)
In [11]:
ii = vIMC(pp)
plot(pp,jj, 'r', label=r'$ I(p) $')
plot(pp,ii, label=r'$ I_{MC}(p) $')
title('Monte Carlo without seed reset')
legend();
xlabel('$p$')
#ylabel('Imc(p), I(p)')
Out[11]:
Generate \(I(p)\) with seed reset
In [12]:
kk = vIMC(pp, True)
plot(pp,jj,'r', label=r'$ I(p) $')
plot(pp,kk, label=r'$ I_{MC}^{(r)}(p) $')
title('Monte Carlo with seed reset')
legend();
xlabel('$p$')
#ylabel('Imc(p, reset=True), I(p)')
Out[12]:
Monte Carlo error \(I_{MC}(p) - I(p)\)
In [13]:
plot(pp,ii-jj)
title('Monte Carlo error')
ylabel('$I_{MC}(p) - I(p)$')
xlabel('$p$')
Out[13]:
Error for Monte Carlo with seed reset: \(I_{MC}^{(r)}(p) - I(p)\)
In [14]:
plot(pp,kk-jj)
title('Error for Monte Carlo with seed reset')
xlabel('$p$')
ylabel('$I_{MC}^{(r)}(p) - I(p)$')
Out[14]:
In [15]:
e = ii-jj
In [16]:
eMC = e.std()
print(eMC)
In []: