import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.special import erf
from scipy.optimize import fsolve
from scipy.integrate import quad
plt.rc('font', size=14)
def make_plot():
μ = 2
σ = 0.5
x = np.linspace(0,4,100)
g = 1/np.sqrt(2*np.pi*σ**2)*np.exp(-1/2/σ**2*(x-μ)**2)
G = 0.5*(1+erf(1/σ/np.sqrt(2)*(x-μ)))
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,5))
ax1.plot(x,g)
ax1.set_xlabel('x'); ax1.set_ylabel('g(x)')
ax1.set_xlim([0,4]); ax1.set_ylim([-0.1,1])
ax2.plot(x,G)
ax2.set_xlabel('x'); ax2.set_ylabel('G(x)')
ax2.set_xlim([0,4]); ax2.set_ylim([-0.1,1])
def Ginv(xFind, Gwant):
return 0.5*(1+erf(1/σ/np.sqrt(2)*(xFind-μ))) - Gwant
for yy in np.linspace(0.01,0.99,15):
xx = fsolve(Ginv, μ, args=yy)[0]
ax2.arrow(4, yy, xx-4+0.1, 0, head_width=0.02, head_length=0.1, fc='gray', ec='gray', linestyle='dashed')
ax2.arrow(xx, yy, 0, -yy+0.05-0.1, head_width=0.05, head_length=0.05, fc='gray', ec='gray')
plt.tight_layout()
make_plot()
n = 2000 # number of points
#---------- Define g(x)
μ = 2
σ = 0.5
maxG = 1.0
def g(x):
return 1/np.sqrt(2*np.pi*σ**2)*np.exp(-1/2/σ**2*(x-μ)**2)
#---------- get g(x) distributed random x points
rG = np.random.rand(n)*maxG # uniform samples on G
rx = np.empty(n) # corresponding x locations found next
def Ginv(xFind, Ggiven): # inverse function: G --> x
return 0.5*(1+erf(1/σ/np.sqrt(2)*(xFind-μ))) - Ggiven
for i, rGi in enumerate(rG):
rx[i] = fsolve(Ginv, μ, args=rGi)
#---------- get y points below g(x)
ry = np.empty(n)
for i in range(n):
ry[i] = np.random.rand()*g(rx[i])
#---------- plot results
x = np.linspace(0,2*μ,100)
fig,ax = plt.subplots(figsize=(10,5))
ax.plot(x,g(x), 'k-');
ax.plot(rx, ry, 'b.', markersize=2);
ax.plot(rx, np.random.rand(len(rx))*0.04-0.02-0.04, '.', markersize=1.0, color='red')
ax.set_xlabel('x');
ax.set_ylabel('g(x)');
rx
points above gives:plt.hist(rx, 40, rwidth=0.8);
def make_plot2():
n = 8000 # random numbers to try (not all kept)
μ = 2
σ = 0.5
def g(x):
return 1/np.sqrt(2*np.pi*σ**2)*np.exp(-1/2/σ**2*(x-μ)**2)
def F(x): # bounding function, here a simple uniform profile
return 1
#-----------------
nkeep = 0 # total number kept: incremented below
rx_all = np.empty(0) # storage for all points: expanded below
for i in range(n):
rx = np.random.rand()*2*μ
Paccept = g(rx)/F(rx)
r = np.random.rand()
if r<Paccept:
rx_all=np.append(rx_all, rx)
nkeep += 1
#-----------------
plt.hist(rx_all, 40, rwidth=0.8);
plt.xlim([0,4])
make_plot2()
| 1 2 3 4 |
|-----|-----|-----|-----|
| dt dt dt dt |
$\phantom{xxx}t_0\phantom{xxxxxxxxxxxxxxxxxx}t_0+\Delta t_e$
Hence, $$P(t_{\text{next event}}>\Delta t_e) = e^{-\lambda \Delta t_e}.$$ Or, $$P(t_{\text{next event}}<\Delta t_e) = 1- e^{-\lambda \Delta t_e}.$$
This is an exponential distribution.
So, to generate exponentially-distributed events, randomly sample $\Delta t_e$ values from this distribution.
$$\Delta t_e = -\frac{\ln(1-P)}{\lambda}.$$
def make_plot3():
λ = 1
n = 2000
rP = np.random.rand(n)
Δt_e = -np.log(1-rP)/λ
x = np.linspace(0,8,1000)
y = λ*np.exp(-λ*x)
times = np.cumsum(Δt_e[:40])
plt.rc('font', size=14);
fig,(ax1,ax2) = plt.subplots(1,2,figsize=(10,5))
ax1.hist(Δt_e, 80, density=True, rwidth=0.8);
ax1.plot(x,y, linewidth=3);
ax1.set_xlabel(r'$\Delta t_e$');
ax1.set_ylabel(r'$p(\Delta t_e)$');
ax1.set_xlim([0,5])
ax2.plot(times, np.zeros(len(times)), '.', markersize=5)
for t in times:
ax2.axvline(x=t, ymin=0.45, ymax=0.55, linewidth=1, color='gray')
ax2.set_xlabel('time')
ax2.set_ylim([-1,1])
plt.tight_layout()
make_plot3()