using Plots
using SpecialFunctions
using NLsolve
using Random
using LaTeXStrings
function make_plot()
μ = 2.0
σ = 0.5
x = LinRange(0,4,100)
g = @. 1/sqrt(2*pi*σ^2)*exp(-1/2/σ^2*(x-μ)^2)
G = @. 0.5*(1+erf(1/σ/sqrt(2)*(x-μ)))
#fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,5))
plot(x,g, lw=3)
plot!(xlabel="x", ylabel="g(x)")
p1=plot!(xlim=[0,4], ylim=[-0.1,1], legend=nothing)
plot(x,G, lw=3)
plot!(xlabel="x", ylabel="G(x)")
p2=plot!(xlim=[0,4], ylim=[-0.1,1], legend=nothing)
function Ginv(xFind, Gwant)
return [0.5*(1+erf(1/σ/sqrt(2)*(xFind[1]-μ))) - Gwant]
end
GR.setarrowsize(0.6)
for yy in LinRange(0.01,0.98,15)
xx = nlsolve((F,x)->F[:]=Ginv(x,yy), [μ]).zero[1]
p2=plot!([4,xx],[yy,yy],arrow=:closed, arrowheadsize=0.1,color=:gray,lw=1,label="")
p2=plot!([xx,xx],[yy,-0.1],arrow=:closed, arrowheadsize=0.1,color=:gray,lw=1,label="")
end
p=plot(p1, p2, layout=(1,2), size=(800,400))
display(p)
end
make_plot()
n = 2000 # number of points
#---------- Define g(x)
μ = 2.
σ = 0.5
maxG = 1.0
g(x) = 1/sqrt(2*pi*σ^2)*exp(-1/2/σ^2*(x-μ)^2)
#---------- get g(x) distributed random x points
rG = rand(n)*maxG # uniform samples on G
rx = zeros(n) # corresponding x locations found next
function Ginv(xFind, Ggiven) # inverse function: G --> x
return [0.5*(1+erf(1/σ/sqrt(2)*(xFind[1]-μ))) - Ggiven]
end
rx[1] = nlsolve((F,x)->F[:]=Ginv(x,rG[1]), [μ]).zero[1]
for (i, rGi) in enumerate(rG)
rx[i] = nlsolve((F,x)->F[:]=Ginv(x,rGi), [μ]).zero[1]
end
#---------- get y points below g(x)
ry = [rand()*g(rx[i]) for i in 1:n]
#---------- plot results
x = LinRange(0,2*μ,100)
plot(x, g.(x), lw=2, color="black")
scatter!(rx, ry, ms=1, color="blue", markerstrokewidth=0)
scatter!(rx, rand(length(rx))*0.04 .- 0.06, ms=1, color="red", markerstrokewidth=0)
plot!(xlabel="x", ylabel="g(x)", legend=nothing, size=(600,300))
rx
points above gives:histogram(rx, bins=40, legend=nothing)
function make_plot2()
n = 8000 # random numbers to try (not all kept)
μ = 2.
σ = 0.5
g(x) = 1/sqrt(2*pi*σ^2)*exp(-1/2/σ^2*(x-μ)^2)
F(x) = 1 # bounding function, here a simple uniform profile
#-----------------
nkeep = 0 # total number kept: incremented below
rx_all = Float64[] # storage for all points: expanded below
for i in 1:n
rx = rand()*2*μ
Paccept = g(rx)/F(rx)
r = rand()
if r<Paccept
append!(rx_all, rx)
nkeep += 1
end
end
#-----------------
histogram(rx_all, bins=40, legend=nothing)
p=plot!(xlim=[0,4])
display(p)
end
make_plot2()
| 1 2 3 4 |
|-----|-----|-----|-----|
| dt dt dt dt |
t0 t0+Dt_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}.$$
function make_plot3()
λ = 1
n = 2000
rP = rand(n)
Δt_e = -log.(1 .-rP)/λ
x = LinRange(0,8,1000)
y = λ*exp.(-λ*x)
times = cumsum(Δt_e[1:40])
histogram(Δt_e, bins=80, normalize=:pdf)
plot!(x,y, lw=3, legend=nothing, color="orange")
plot!(xlabel=L"\Delta t_e")
plot!(ylabel=L"p(\Delta t_e)")
p1=plot!(xlim=[0,5], ylim=[0,1.1])
plot()
for t in times
plot!([t,t],[-0.15,0.15], lw=1, color="gray")
end
scatter!(times, zeros(length(times)), ms=2, markerstrokewidth=0)
p2=plot!(xlabel="time", ylim=[-1,1], legend=nothing)
plot!(p1,p2, layout=(1,2), size=(600,300))
end
make_plot3()