added convergence plots of main poisson to main.tex

experiments/tests_period_7/local_main_poisson.jl

@@ -1,3 +1,9 @@
+using Plots
+using Statistics
+using Distributions
+using GLM, DataFrames
+
+
 function Y(t, sig, A::Function, f::Function, y0)
  (s = -log(rand()) / sig; sol = y0)
  while s < t
@@ -7,8 +13,6 @@ function Y(t, sig, A::Function, f::Function, y0)
  sol
 end
 
-using Plots
-using Statistics
 sig = 1.0
 A(s) = s
 f(s) = -s^3 + s
@@ -68,7 +72,74 @@ plot(ttss, yyss, st=:scatter, xscale=:log10, yscale=:log10, label="Y(t)")
 
 
 
-while false
- println("test")
+# plot of realizations
+using Plots
+using Statistics
+
+
+function Yplot(t, sig, A::Function, f::Function, y0)
+ (s = -log(rand()) / sig; sol = y0)
+ path = []
+ while s < t
+ push!(path, (s, sol))
+ sol += (A(s) * sol .+ f(s)) ./ sig
+ s -= log(rand()) / sig
+ end
+ if length(path) > 0
+ push!(paths, path)
+ end
+ sol
 end
-println("hah")
+
+A(s) = s
+f(s) = -s^3 + s
+q = 1.0
+sol(s) = s^2 + 1
+
+# A(s) = 1
+# f(s) = 0
+# q = 1.0
+# sol(s) = exp(s)
+
+t = 1.0
+sig = 1000.0
+nsim1 = 2 * 10^1
+paths = []
+for _ in 1:nsim1
+ Yplot(t, sig, A, f, q)
+end
+
+# Collect all points from all paths into two arrays
+x_values = Float64[]
+y_values = Float64[]
+
+for path in paths
+ append!(x_values, [p[1] for p in path])
+ append!(y_values, [p[2] for p in path])
+end
+
+# Perform OLS regression
+ols = lm(@formula(y ~ x + x^2), DataFrame(x=x_values, y=y_values))
+
+
+
+tt = range(0, t, length=100)
+pl = plot()
+
+for path in paths
+ # plot!(pl, [p[1] for p in path], [p[2]-sol(p[1]) for p in path], label="", st=:scatter, color=:blue, alpha=0.5)
+ plot!(pl, [p[1] for p in path], [p[2] - sol(p[1]) for p in path], label="", alpha=0.5)
+end
+# plot!(pl, tt, sol.(tt), label="sol(t)", linewidth=6, color=:orange)
+plot!(pl, tt, predict(ols, DataFrame(x=tt)) - sol.(tt), label="OLS y ~ x+x^2", linewidth=3, color=:red)
+title!("error of different realizations")
+display(pl)
+
+
+histogram([p[1] for p in points], bins=100, label="Y(t)")
+histogram([p[end] for p in paths], bins=100, label="Y(t)")
+
+
+endpoints = [p[end][2] for p in paths]
+n = length(endpoints)
+plot(sort(randn(n)), sort(endpoints), st=:scatter, alpha=0.3, markersize=3)

latex/main paper/main.tex

@@ -168,7 +168,7 @@ \subsection{Monte Carlo Integration}
 \end{related}
 
 For smooth $1$-dimensional integration, the linear trade-off between
-variance and average simulation time is not even close to optimal see
+variance and average simulation time is not close to optimal see
 Theorem \ref{thrm:order trap}. \\
 When it comes to comparing better trade-offs,
 Information-Based Complexity (IBC) is often employed.
@@ -194,7 +194,8 @@ \subsection{Monte Carlo Integration}
 
 \subsection{Recursive Monte Carlo}
 In this subsection, we introduce Recursive Monte Carlo (RMC)
-with an example of initial value problem (IVP) and Walk on Spheres.
+with an example of an Initial Value Problem (IVP) and
+Walk on Spheres.
 
 \begin{notation}[$U,U_{j}$]
  We will frequently use the uniform distribution, so we will abbreviate it
@@ -786,14 +787,34 @@ \subsection{Recursion}
  Sampling $\tau$ uniformly is equivalent to sampling the next event in a Poisson process.
 \end{definition}
 
-\begin{julia}[implementation of Example \ref{def:main poisson}] \label{jl:main poisson}
+\begin{julia}[implementation of \ref{def:main poisson}] \label{jl:main poisson}
  (\ref{eq:poisson main}) can be turned into a RRVE by sampling $\tau \sim U$, no Russian roulette is needed for
  termination because sampling the first integral ends recursion. No recursion is used
  for the implementation as we can sample the Poisson process in the reverse order.
 
  \juliacode{julia_code/main_poisson.jl}
 \end{julia}
 
+\begin{figure}[h!]
+ \centering
+ \includegraphics[width=\textwidth]{julia_plots/main_poisson_error.pdf}
+ \caption{
+ Intermediate points of $5$ realizations of the error ($Y(t)-y(t)$) of \ref{jl:main poisson} for (\ref{couple recu ex1}) and (\ref{couple recu ex2}) with $a=1, \sigma = \{10,1000\}$.
+ Regression over the realizations and time is done with ordinary least square onto
+ a second degree polynomial. The uniform qqplot of the time realizations indicate that the time realizations
+ are uniformly distributed.
+ }
+ \label{fig:main poisson error}
+\end{figure}
+
+\begin{figure}[h!]
+ \centering
+ \includegraphics[width=\textwidth]{julia_plots/main_poisson_convergence.pdf}
+ \caption{Realizations of norm of the error of (split) \ref{jl:main poisson} for (\ref{couple recu ex1}) and (\ref{couple recu ex2}) with $a=1$,
+ $nsim =1$ and $\sigma$ variable or $nsim$ variable and $\sigma = 1$.}
+ \label{fig:main poisson convergence}
+\end{figure}
+
 \begin{related}[main Poisson]
  A very similar estimator is used to solve ODEs derived from the
  $1$-dimensional telegraph equations in \cite{acebron_monte_2016}.

latex/main paper/plots_julia_code/main_poisson_convergence.jl

@@ -0,0 +1,51 @@
+using Plots
+using Random
+using LinearAlgebra
+using Plots.PlotMeasures
+
+function Y(t, sig, A::Function, f::Function, y0)
+ (s = -log(rand()) / sig; sol = y0)
+ while s < t
+ sol += (A(s) * sol .+ f(s)) ./ sig
+ s -= log(rand()) / sig
+ end
+ sol
+end
+
+a = 1
+A(s) = [a 0; 1 a]
+q = [1, 0]
+f(s) = zero(q)
+sol(s) = [exp(a * s), s * exp(a * s)]
+sig = 1.0
+t = 1.0
+
+
+Random.seed!(2234)
+xticks = 10.0 .^ range(1, 4, step=1)
+yticks = 10.0 .^ range(-4, 0, step=1)
+sigs = exp10.(range(1, 4, length=1000))
+errors = [norm(Y(t, sig, A, f, q) - sol(t)) for sig in sigs]
+p1 = plot(sigs, errors, st=:scatter, xscale=:log10, yscale=:log10, label="error", alpha=0.5)
+plot!(p1, sigs, 5 * sigs .^ -0.5, label="\$O(sig^{-0.5})\$", linewidth=4)
+plot!(p1, sigs, 5 * sigs .^ -1, label="\$O(sig^{-1})\$", linewidth=4)
+xticks!(p1, xticks) # Set xticks
+yticks!(p1, yticks) # Set yticks
+title!(p1, "error vs sig")
+xlabel!(p1, "sig")
+ylabel!(p1, "error")
+
+nsims = exp10.(range(1, 4, length=1000))
+errors = [norm(sum(Y(t, sig, A, f, q) for _ in 1:nsim) / nsim - sol(t, sig)) for nsim in nsims]
+p2 = plot(nsims, errors, st=:scatter, xscale=:log10, yscale=:log10, label="error", alpha=0.5)
+plot!(p2, nsims, 5 * nsims .^ -0.5, label="\$O(nsim^{-0.5})\$", linewidth=4)
+plot!(p2, nsims, 5 * nsims .^ -1, label="\$O(nsim^{-1})\$", linewidth=4)
+xticks!(p2, xticks) # Set xticks
+yticks!(p2, yticks) # Set yticks
+title!(p2, "error vs nsim")
+xlabel!(p2, "nsim")
+ylabel!(p2, "error")
+
+p = plot(p1, p2, layout=(1, 2), size=(1000, 400), bottom_margin=5mm, left_margin=5mm) # Combine both plots
+display(p)
+savefig(p, "latex/main paper/julia_plots/main_poisson_convergence.pdf")

latex/main paper/plots_julia_code/main_poisson_error.jl

@@ -0,0 +1,71 @@
+using Plots
+using Random
+using GLM, DataFrames
+using Plots.PlotMeasures
+
+function Yplot(t, sig, A::Function, f::Function, y0)
+ (s = -log(rand()) / sig; sol = y0)
+ path = []
+ while s < t
+ push!(path, (s, sol))
+ sol += (A(s) * sol .+ f(s)) ./ sig
+ s -= log(rand()) / sig
+ end
+ if length(path) > 0
+ push!(paths, path)
+ end
+ sol
+end
+
+a = 1
+A(s) = [a 0; 1 a]
+q = [1, 0]
+f(s) = zero(q)
+sol(s) = [exp(a * s), s * exp(a * s)]
+
+Random.seed!(2234)
+paths = []
+t = 1.0
+sigs = [10.0, 1000.0]
+nsim1 = 5
+
+markers = [:circle, :square, :diamond, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :hexagon, :heptagon, :octagon, :star4, :star5, :star6, :star7, :star8, :vline, :hline]
+colors = [:red, :green, :blue, :orange, :purple, :brown, :black, :grey]
+pl = plot(layout=(2, 3), size=(1300, 800), xlim=(0, 1),
+ bottom_margin=5mm, left_margin=5mm, right_margin=4mm) # Change layout to 2 rows and 3 columns
+ylims = [(-1, 1), (-1 / 4, 1 / 4)]
+for i in [1, 2]
+ paths = []
+ sig = sigs[i]
+ y = sum(Yplot(t, sig, A, f, q) for _ in 1:nsim1) / nsim1
+ for c in [1, 2]
+ ylims!(pl[(i-1)*3+c], ylims[i])
+ xx = [p[1] for path in paths for p in path]
+ yy = [p[2][c] for path in paths for p in path]
+ ols = lm(@formula(y ~ x + x^2), DataFrame(x=xx, y=yy))
+ Random.seed!(6421)
+ for path in paths
+ marker = rand(markers)
+ color = rand(colors)
+ if i == 1
+ for p in path
+ scatter!(pl[(i-1)*3+c], [p[1]], [p[2][c] - sol(p[1])[c]], label="", marker=marker, color=color, alpha=0.5)
+ end
+ end
+ plot!(pl[(i-1)*3+c], [p[1] for p in path], [p[2][c] - sol(p[1])[c] for p in path], label="", alpha=0.5, color=color, linewidth=2)
+ xlabel!(pl[(i-1)*3+c], "time")
+ ylabel!(pl[(i-1)*3+c], "error")
+ end
+ tt = range(0, t, length=100)
+ hline!(pl[(i-1)*3+c], [0], label="0 error", color=:black, linewidth=3)
+ plot!(pl[(i-1)*3+c], tt, predict(ols, DataFrame(x=tt)) .- [sol(t)[c] for t in tt], label="OLS y ~ x+x^2", linewidth=3, color=:red)
+ title!(pl[(i-1)*3+c], "error, sig=$sig, component=$c", titiefontsize=10)
+ #qqplot of time
+ plot!(pl[(i-1)*3+3], sort(rand(length(xx))), sort(xx), title="uniform qqplot of time sig=$sig", label="", linewidth=3)
+ xlabel!(pl[(i-1)*3+3], "uniform")
+ ylabel!(pl[(i-1)*3+3], "Poisson process")
+ end
+end
+
+display(pl)
+savefig(pl, "latex/main paper/julia_plots/main_poisson_error.pdf")

latex/main paper/plots_julia_code/trap_example.jl

@@ -20,12 +20,12 @@ end
 Random.seed!(4181)
 f(x) = exp(x)
 nn = round.(Int, exp10.(range(0, 6.5, length=400)))
-sol = exp(1) - 1
-trapezium_error = abs.(trapezium.(f, nn) .- sol)
-mc_trapezium_error = abs.(MCtrapezium.(f, nn, 100) .- sol)
+sol1 = exp(1) - 1
+trapezium_error = abs.(trapezium.(f, nn) .- sol1)
+mc_trapezium_error = abs.(MCtrapezium.(f, nn, 100) .- sol1)
 
 yticks = 10.0 .^ range(-16, 0, step=2) # adjust the range and step as needed
-p = plot(nn, trapezium_error .+ eps(), label="OG trap", xscale=:log10, yscale=:log10, legendfontsize=12, yticks=yticks, xlabel="n", ylabel="Error", bottom_margin=2mm)
+p = plot(nn, trapezium_error .+ eps(), label="OG trap", xscale=:log10, yscale=:log10, legendfontsize=12, yticks=yticks, xlabel="n", ylabel="error", bottom_margin=2mm)
 plot!(p, nn, mc_trapezium_error .+ eps(), label="MC trap", xscale=:log10, yscale=:log10)
 plot!(p, nn, nn .^ -2 .+ eps(), label="\$O(n^{-2})\$", xscale=:log10, yscale=:log10, linestyle=:dash)
 plot!(p, nn, nn .^ -2.5 .+ eps(), label="\$O(n^{-2.5})\$", xscale=:log10, yscale=:log10, linestyle=:dash)