Testing the confidence regions and intervals 2

We consider a simple and quick-to-solve Random ODE to test the confidence regions and intervals. With a simple model, we can easily run a million simulations to test the statistics.

The Random ODE is a simple homogeneous linear equation in which the coefficient is a Wiener process and for which we know the distribution of the exact solution.

Now we consider an arbitrary number of mesh resolutions.

The equation

We consider the RODE

\[ \begin{cases} \displaystyle \frac{\mathrm{d}X_t}{\mathrm{d} t} = W_t X_t, \qquad 0 \leq t \leq T, \\ \left. X_t \right|_{t = 0} = X_0, \end{cases}\]

where $\{W_t\}_{t\geq 0}$ is a standard Wiener process. The explicit solution is

\[ X_t = e^{\int_0^t W_s \;\mathrm{d}s} X_0.\]

As seen in the first example of this documentation, once an Euler approximation is computed, along with realizations $\{W_{t_i}\}_{i=0}^n$ of a sample path of the noise, we consider an exact sample solution given by

\[ X_{t_j} = X_0 e^{\sum_{i = 0}^{j-1}\left(\frac{1}{2}\left(W_{t_i} + W_{t_{i+1}}\right)(t_{i+1} - t_i) + Z_i\right)},\]

for realizations $Z_i$ drawn from a normal distribution and scaled by the standard deviation $\sqrt{(t_{i+1} - t_i)^3/12}$. This is implemented by computing the integral recursively, via

\[ \begin{cases} I_j = I_{j-1} + \frac{1}{2}\left(W_{t_{j-1}} + W_{t_j}\right)(t_{j} - t_{j-1}) + Z_j, \\ Z_j = \sqrt{\frac{(t_{j} - t_{j-1})^3}{12}} R_j, \\ R_j \sim \mathcal{N}(0, 1), \\ \end{cases}\]

with $I_0 = 0$, and setting

\[ X_{t_j} = X_0 e^{I_j}.\]

Setting up the problem

First we load the necessary packages

using Plots
using Random
using Distributions
using RODEConvergence

Then we set up some variables, starting by choosing the Xoshiro256++ pseudo-random number generator, and setting its seed for the sake of reproducibility:

rng = Xoshiro(123)

We set the right hand side of the equation:

f(t, x, y, p) = y * x

Next we set up the time interval and the initial distribution law for the initial value problem, which we take it to be a standard Distributions.Normal random variable:

t0, tf = 0.0, 1.0
x0law = Normal()

Distributions.Normal{Float64}(μ=0.0, σ=1.0)

The noise is a WienerProcess starting at $y_0 = 0$:

y0 = 0.0
noise = WienerProcess(t0, tf, y0)

WienerProcess{Float64}(0.0, 1.0, 0.0)

There is no parameter in the equation, so we just set params to nothing.

params = nothing

The number of mesh points for the target solution and the approximations

ntgt = 2^12
ns = 2 .^ (4:2:10)

4-element Vector{Int64}:
   16
   64
  256
 1024

Notice we just chose two mesh sizes, so we can easily visualize the distributions.

The target solution as described above is implemented as

target_solver! = function (xt::Vector{T}, t0::T, tf::T, x0::T, f::F, yt::Vector{T}, params::Q, rng::AbstractRNG) where {T, F, Q}
    axes(xt) == axes(yt) || throw(
        DimensionMismatch("The vectors `xt` and `yt` must match indices")
    )

    n = size(xt, 1)
    dt = (tf - t0) / (n - 1)
    i1 = firstindex(xt)
    xt[i1] = x0
    integral = zero(T)
    zscale = sqrt(dt^3 / 12)
    for i in Iterators.drop(eachindex(xt, yt), 1)
        integral += (yt[i] + yt[i1]) * dt / 2 + zscale * randn(rng)
        xt[i] = x0 * exp(integral)
        i1 = i
    end
end

and with that we construct the CustomMethod that solves the problem with this target_solver!:

target = CustomUnivariateMethod(target_solver!, rng)

The method for which we want to estimate the rate of convergence is, naturally, the Euler method, denoted RandomEuler:

method = RandomEuler()

RandomEuler{Float64, Distributions.Univariate}(Float64[])

Defining helper functions

We first write some helper functions to grab the statistics, print some information, and build some plots.

function getstatistics(rng, suite, ns, nk, m)
    ps = zeros(nk)
    pmins = zeros(nk)
    pmaxs = zeros(nk)
    allerrors = zeros(nk, length(ns))
    allstderrs = zeros(nk, length(ns))
    @time for k in 1:nk
        resultk = solve(rng, suite)
        ps[k] = resultk.p
        allerrors[k, :] .= resultk.errors
        allstderrs[k, :] .= resultk.stderrs
        pmins[k] = resultk.pmin
        pmaxs[k] = resultk.pmax
    end
    meanerror = mean(allerrors, dims=1)
    pmean = mean(ps)

    percent_ei_in = 100 * [ count(( meanerror[i] .> allerrors[:, i] .- 1.96allstderrs[:, i] ) .& ( meanerror[i] .< allerrors[:, i] .+ 1.96allstderrs[:, i] )) for i in eachindex(axes(meanerror, 2), axes(allstderrs, 2)) ] / nk

    clevel = 0.95 # 95% confidence level
    elevel = clevel^(1/length(suite.ns)) # assuming independence
    elevel = 1.0 - ( 1.0 - clevel ) / length(suite.ns) # not assuming independence, using Bonfaroni inequality
    escore = quantile(Normal(), (1 + elevel) / 2)

    percent_e_in = 100 * count(
        reduce(
            &,
            ae - escore * ast .< meanerror' .< ae + escore * ast
        ) for (ae, ast) in zip(eachrow(allerrors), eachrow(allstderrs))
    ) / nk

    nkji = div(size(allerrors,1), size(allerrors,2))
    allerrorssplit = hcat((allerrors[nkji*(i-1)+1:nkji*i, i] for i in axes(allerrors, 2))...)
    allstderrssplit = hcat((allstderrs[nkji*(i-1)+1:nkji*i, i] for i in axes(allstderrs, 2))...)

    percent_e_split_in = 100 * count(
        reduce(
            &,
            ae - escore * ast .< meanerror' .< ae + escore * ast
        ) for (ae, ast) in zip(eachrow(allerrorssplit), eachrow(allstderrssplit))
    ) / nkji

    allerrorsdealigned = hcat([circshift(allerrors[:, i], i-1) for i in 1:4]...)
    allstderrsdealigned = hcat([circshift(allstderrs[:, i], i-1) for i in 1:4]...)

    percent_e_dealigned_in = 100 * count(
        reduce(
            &,
            ae - escore * ast .< meanerror' .< ae + escore * ast
        ) for (ae, ast) in zip(eachrow(allerrorsdealigned), eachrow(allstderrsdealigned))
    ) / nk

    deltas = (suite.tf - suite.t0) ./ suite.ns
    A = [one.(deltas) log.(deltas)]
    L = inv(A' * A) * A'

    Llnerrors = L * log.(allerrors')

    Llnerrorsdealigned = L * log.(allerrorsdealigned')

    percent_p_dealigned_in = 100 * count( ( pmean .> Llnerrorsdealigned[2, :] .- (ps .- pmins) ) .& ( pmean .< Llnerrorsdealigned[2, :] .+ (pmaxs .- ps) ) ) / nk

    percent_p_in = 100 * count(( pmean .> pmins ) .& ( pmean .< pmaxs )) / nk

    pstd = std(ps)
    percent_p_alt_in = 100 * count(( pmean .> ps .- 1.96pstd ) .& ( pmean .< ps .+ 1.96pstd )) / nk

    pdlgnstd = std(Llnerrorsdealigned[2, :])
    percent_p_alt_dealigned_in = 100 * count(( pmean .> Llnerrorsdealigned[2, :] .- 1.96pdlgnstd ) .& ( pmean .< Llnerrorsdealigned[2, :] .+ 1.96pdlgnstd )) / nk

    return ps, escore, allerrors, allstderrs, allerrorssplit, allerrorsdealigned, meanerror, pmean, Llnerrors, Llnerrorsdealigned, percent_p_in, percent_p_dealigned_in, percent_p_alt_dealigned_in, percent_ei_in, percent_e_in, percent_e_split_in, percent_e_dealigned_in, L
end

function printpercents(
    percent_p_in, percent_p_dealigned_in, percent_p_alt_dealigned_in, percent_ei_in, percent_e_in, percent_e_split_in, percent_e_dealigned_in
)
    println("percent p in: $percent_p_in%")
    println("percent p dealigned in: $percent_p_dealigned_in%")
    println("percent p alt dealigned in: $percent_p_alt_dealigned_in%")
    for i in eachindex(percent_ei_in)
        println("percent E$i in: $(percent_ei_in[i])%")
    end
    println("percent E in: $percent_e_in%")
    println("percent E in split: $percent_e_split_in%")
    println("percent E in dealigned: $percent_e_dealigned_in%")
end

function showplots(
    ps, escore, allerrors, allerrorssplit, allerrorsdealigned, Llnerrors, Llnerrorsdealigned, pmean, result, m, nk, percent_ei_in, percent_e_in, percent_e_split_in, percent_p_dealigned_in, percent_e_dealigned_in, L
)

    rect = Shape(
        [
            (result.errors[1] - escore * result.stderrs[1], result.errors[2] - escore * result.stderrs[2]),
            (result.errors[1] - escore * result.stderrs[1], result.errors[2] + escore * result.stderrs[2]),
            (result.errors[1] + escore * result.stderrs[1], result.errors[2] + escore * result.stderrs[2]),
            (result.errors[1] + escore * result.stderrs[1], result.errors[2] - escore * result.stderrs[2])
        ]
    )

    plt_errors = plot(title="Errors ϵ₁ and ϵ₂(m=$m, nk=$nk)", titlefont=10, xlabel="ϵ₁", ylabel="ϵ₂")
    begin
        scatter!(plt_errors, allerrors[:, 1], allerrors[:, 2], alpha=0.2, label="errors ($(round(percent_e_in, digits=2))% in CI)")
        scatter!(plt_errors, allerrorsdealigned[:, 1], allerrorsdealigned[:, 2], alpha=0.2, label="errors dealigned ($(round(percent_e_dealigned_in, digits=2))% in CI)")
        scatter!(plt_errors, Tuple(mean(view(allerrors, :, 1:2), dims=1)), markersize=4, label="error mean")
        plot!(plt_errors, rect, alpha=0.2, label="CI")
    end

    plt_errors_split = plot(title="Errors split (m=$m, nk=$nk) \n ($(round(percent_e_split_in, digits=2))% in CI)", titlefont=10, xlabel="ϵ₁", ylabel="ϵ₂")
    begin
        scatter!(plt_errors_split, allerrorssplit[:, 1], allerrorssplit[:, 2], alpha=0.2, label="errors")
        scatter!(plt_errors_split, Tuple(mean(view(allerrors, :, 1:2), dims=1)), markersize=4, label="error mean")
        plot!(plt_errors_split, rect, alpha=0.2, label="CI")
    end

    plt_errors_dealigned = plot(title="Errors dealigned (m=$m, nk=$nk) \n ($(round(percent_e_dealigned_in, digits=2))% in CI)", titlefont=10, xlabel="ϵ₁", ylabel="ϵ₂")
    begin
        scatter!(plt_errors_dealigned, allerrorsdealigned[:, 1], allerrorsdealigned[:, 2], alpha=0.2, label="errors")
        scatter!(plt_errors_dealigned, Tuple(mean(view(allerrors, :, 1:2), dims=1)), markersize=4, label="error mean")
        plot!(plt_errors_dealigned, rect, alpha=0.2, label="CI")
    end

    plt_errors3d = plot(title="Errors ϵ₁, ϵ₂, and ϵ₃ (m=$m, nk=$nk)", titlefont=10, xlabel="ϵ₁", ylabel="ϵ₂", zlabel="ϵ₃")
    begin
        scatter3d!(plt_errors3d, allerrors[:, 1], allerrors[:, 2], allerrors[:, 3], alpha=0.2, label="errors ($(round(percent_e_in, digits=2))% in CI)")
        scatter3d!(plt_errors3d, allerrorsdealigned[:, 1], allerrorsdealigned[:, 2], allerrorsdealigned[:, 3], alpha=0.2, label="errors dealigned ($(round(percent_e_dealigned_in, digits=2))% in CI)")
        scatter3d!(plt_errors3d, Tuple(mean(view(allerrors, :, 1:3), dims=1)), markersize=4, label="error mean")
    end

    plt_hist_e = [
        begin
            plot(title="Histogram of ϵ$i (m=$m, nk=$nk) \n ($(round(percent_ei_in[i], digits=2))% in CI)", titlefont=10, xlabel="ϵ$i")
            histogram!(allerrors[:, i], label="error ϵ$i")
            vline!([mean(allerrors[:, i])], color=:steelblue, linewidth=4, label="mean")
            vline!([result.errors[i]], label="sample")
            vline!([result.errors[i] - 2result.stderrs[i], result.errors[i] + 2result.stderrs[i]], label="CI from sample")
        end
        for i in eachindex(axes(allerrors, 2))
    ]

    #sn = 50
    #s1 = L * log.(max.(0.0, [result.errors[1] .+ escore * result.stderrs[1] * range(-1, 1, length=sn) ( result.errors[2] - escore * result.stderrs[2] ) .* ones(sn)]'))
    #s2 = L * log.(max.(0.0, [( result.errors[1] + escore * result.stderrs[1] ) .* ones(sn) result.errors[2] .+ escore * result.stderrs[2] * range(-1, 1, length=sn)]'))
    #s3 = L * log.(max.(0.0, [result.errors[1] .+ escore * result.stderrs[1] * reverse(range(-1, 1, length=sn)) ( result.errors[2] + escore * result.stderrs[2] ) .* ones(sn)]'))
    #s4 = L * log.(max.(0.0, [( result.errors[1] - escore * result.stderrs[1] ) .* ones(sn) result.errors[2] .+ escore * result.stderrs[2] * range(-1, 1, length=sn)]'))
    #sides = hcat(s1, s2, s3, s4)

    temean = L * log.(mean(allerrors, dims=1)')

    plt_Cp = plot(title="(C, p) sample from ϵ (m=$m, nk=$nk)", titlefont=10, xlabel="C", ylabel="p")
    begin
        scatter!(plt_Cp, Llnerrors[1, :], Llnerrors[2, :], alpha=0.2, label="correlated")
        scatter!(plt_Cp, Llnerrorsdealigned[1, :], Llnerrorsdealigned[2, :], alpha=0.2, label="dealigned")
        #plot!(plt_Cp, sides[1, :], sides[2, :], label="transformed errors CI")
        scatter!(plt_Cp, Tuple(temean), markersize=4, color=:orange, label="transformed error mean")
        hline!(plt_Cp, [pmean], label="p mean")
        hline!(plt_Cp, [result.pmin, result.pmax], label="sample p CI ($(round(percent_p_dealigned_in, digits=2))% in CI)")
        hline!(plt_Cp, [result.p], label="sample p")
    end

    plt_hist_p = plot(title="Histogram of p (m=$m, nk=$nk) \n ($(round(percent_p_dealigned_in, digits=2))% in CI)", titlefont=10, xlabel="ϵ₁")
    begin
        histogram!(plt_hist_p, Llnerrorsdealigned[2, :], label="p dealigned")
        histogram!(plt_hist_p, ps, label="p")
        vline!(plt_hist_p, [pmean], linewidth=4, label="p mean")
        vline!(plt_hist_p, [result.pmin, result.pmax], label="sample p CI ($(round(percent_p_dealigned_in, digits=2))% in CI)")
        vline!(plt_hist_p, [result.p], label="sample p")
    end

    plts = (
        errors = plt_errors,
        split = plt_errors_split,
        dealigned = plt_errors_dealigned,
        errors3d = plt_errors3d,
        hist = plt_hist_e,
        cp = plt_Cp,
        histp = plt_hist_p
    )

    return plts
end

showplots (generic function with 1 method)

Statistics

Now, with the helper functions, we run a loop varying the number $m$ of samples in each run and the number $nk$ of test runs, showing some relevant statistics.

ms = (400, 800, 1600)
nks = (1000, 1000, 1000)

@assert all(iseven, nks)

allplts = Any[]

for (nrun, m, nk) in zip(eachindex(ms), ms, nks)

    @info "==="
    @info "Run $nrun with m=$m and nk=$nk"
    suite = ConvergenceSuite(t0, tf, x0law, f, noise, params, target, method, ntgt, ns, m)

    ps, escore, allerrors, allstderrs, allerrorssplit, allerrorsdealigned, meanerror, pmean, Llnerrors, Llnerrorsdealigned, percent_p_in, percent_p_dealigned_in, percent_p_alt_dealigned_in, percent_ei_in, percent_e_in, percent_e_split_in, percent_e_dealigned_in, L = getstatistics(rng, suite, ns, nk, m)

    @show cor(allerrors) # strongly correlated!

    @show cor(allerrorssplit) # weakly correlated

    @show cor(allerrorsdealigned) # weakly correlated

    @show cor(Llnerrors') # somehow correlated

    @show cor(Llnerrorsdealigned') # weakly correlated

    printpercents(percent_p_in, percent_p_dealigned_in, percent_p_alt_dealigned_in, percent_ei_in, percent_e_in, percent_e_split_in, percent_e_dealigned_in)

    result = solve(rng, suite)

    plts = showplots(ps, escore, allerrors, allerrorssplit, allerrorsdealigned, Llnerrors, Llnerrorsdealigned, pmean, result, m, nk, percent_ei_in, percent_e_in, percent_e_split_in, percent_p_dealigned_in, percent_e_dealigned_in, L)

    append!(allplts, [plts])
end

[ Info: ===
[ Info: Run 1 with m=400 and nk=1000
 26.449965 seconds (172.11 k allocations: 206.868 MiB, 0.18% gc time)
cor(allerrors) = [1.0 0.9705505172416334 0.9643541989718961 0.9621934611085015; 0.9705505172416334 1.0 0.9719800031589322 0.9656392560760294; 0.9643541989718961 0.9719800031589322 1.0 0.974224307854825; 0.9621934611085015 0.9656392560760294 0.974224307854825 1.0]
cor(allerrorssplit) = [1.0 -0.0025306980685457737 0.09734647609277154 0.017989350057743023; -0.0025306980685457737 1.0 0.03326863459545171 0.002960563104343979; 0.09734647609277154 0.03326863459545171 1.0 0.02242727328070293; 0.017989350057743023 0.002960563104343979 0.02242727328070293 1.0]
cor(allerrorsdealigned) = [1.0 0.012353046792147404 0.003915869842998666 -0.04138440572727088; 0.012353046792147404 1.0 -0.002367455575510119 0.007290954805841701; 0.003915869842998666 -0.002367455575510119 1.0 -0.0016691741448486693; -0.04138440572727088 0.007290954805841701 -0.0016691741448486693 1.0]
cor(Llnerrors') = [1.0 0.12584535004221312; 0.12584535004221312 1.0]
cor(Llnerrorsdealigned') = [1.0 0.9527747119161571; 0.9527747119161571 1.0]
percent p in: 100.0%
percent p dealigned in: 99.9%
percent p alt dealigned in: 95.1%
percent E1 in: 90.3%
percent E2 in: 90.0%
percent E3 in: 90.1%
percent E4 in: 90.6%
percent E in: 92.2%
percent E in split: 80.8%
percent E in dealigned: 81.7%
[ Info: ===
[ Info: Run 2 with m=800 and nk=1000
 53.636477 seconds (172.00 k allocations: 206.863 MiB, 0.05% gc time)
cor(allerrors) = [1.0 0.9669697627885926 0.9617651711382407 0.9571961313452868; 0.9669697627885926 1.0 0.9721375828510537 0.9620252283661136; 0.9617651711382407 0.9721375828510537 1.0 0.9704562109515474; 0.9571961313452868 0.9620252283661136 0.9704562109515474 1.0]
cor(allerrorssplit) = [1.0 0.016007593342482078 -0.10278286886810155 -0.031129141200714073; 0.016007593342482078 1.0 0.016320472436240463 0.0672565677902107; -0.10278286886810155 0.016320472436240463 1.0 0.005980291576012199; -0.031129141200714073 0.0672565677902107 0.005980291576012199 1.0]
cor(allerrorsdealigned) = [1.0 0.025859150977547804 -0.02076452987680885 -0.06302663997582882; 0.025859150977547804 1.0 0.021329369937355454 -0.03702961015347029; -0.02076452987680885 0.021329369937355454 1.0 0.022819977000100022; -0.06302663997582882 -0.03702961015347029 0.022819977000100022 1.0]
cor(Llnerrors') = [1.0 0.09434209905968942; 0.09434209905968942 1.0]
cor(Llnerrorsdealigned') = [1.0 0.9557771117564445; 0.9557771117564445 1.0]
percent p in: 100.0%
percent p dealigned in: 100.0%
percent p alt dealigned in: 94.2%
percent E1 in: 92.5%
percent E2 in: 92.3%
percent E3 in: 92.1%
percent E4 in: 91.9%
percent E in: 95.4%
percent E in split: 88.8%
percent E in dealigned: 87.3%
[ Info: ===
[ Info: Run 3 with m=1600 and nk=1000
105.525813 seconds (172.00 k allocations: 206.863 MiB, 0.01% gc time)
cor(allerrors) = [1.0 0.9693307552466505 0.9602301613993494 0.9583449495799153; 0.9693307552466505 1.0 0.9670334452892287 0.9645906145837028; 0.9602301613993494 0.9670334452892287 1.0 0.9697758863055619; 0.9583449495799153 0.9645906145837028 0.9697758863055619 1.0]
cor(allerrorssplit) = [1.0 -0.023744608648916678 -0.0031170699517114376 -0.02326692334400742; -0.023744608648916678 1.0 0.16023960262516954 -0.025227488823415058; -0.0031170699517114376 0.16023960262516954 1.0 -0.021536232268612343; -0.02326692334400742 -0.025227488823415058 -0.021536232268612343 1.0]
cor(allerrorsdealigned) = [1.0 -0.03367578841166664 -0.04875731666505965 0.02609664742074762; -0.03367578841166664 1.0 -0.0327049515570865 -0.043973488506517285; -0.04875731666505965 -0.0327049515570865 1.0 -0.02691280494602071; 0.02609664742074762 -0.043973488506517285 -0.02691280494602071 1.0]
cor(Llnerrors') = [1.0 0.13993619131030968; 0.13993619131030968 1.0]
cor(Llnerrorsdealigned') = [1.0 0.9539650727071125; 0.9539650727071125 1.0]
percent p in: 100.0%
percent p dealigned in: 100.0%
percent p alt dealigned in: 94.7%
percent E1 in: 93.6%
percent E2 in: 94.0%
percent E3 in: 94.1%
percent E4 in: 93.0%
percent E in: 96.5%
percent E in split: 91.2%
percent E in dealigned: 92.0%

Visualizations

We now visualize some statistics, including histograms and sample distribution. In those plots, we report the percentage of confidence intervals and regions that include the mean.

Histograms of the marginal strong errors

We start with the histograms of each of the strong errors at each mesh resolution. These are the marginals of the joint distribution $(\epsilon_1, \ldots, \epsilon_{i_{\max}).$

plot(size=(800, 400*div(length(ns), 2)), allplts[1].hist...)

plot(size=(800, 400*div(length(ns), 2)), allplts[2].hist...)

plot(size=(800, 400*div(length(ns), 2)), allplts[3].hist...)

Density of the joint distribution of strong errors of the transformed distributions

First we visualise the joint distribution of the samples of the first three strong errors $(\epsilon_1, \epsilon_2, \epsilon_3).$

plot(allplts[1].errors3d)

plot(allplts[2].errors3d)

plot(allplts[3].errors3d)

Now we plot, on the left panel, the sample points of the joint distribution $(\epsilon_1, \epsilon_2)$ of the strong errors of the first two mesh resolutions.

On the right panel, we see the corresponding transformed samples $(C, p) = (A^{\textrm{tr}}A)^{-1}A^{\textrm{tr}}(\epsilon_1, \ldots, \epsilon_{i_{\max}})$.

The confidence interval for the order of convergence $p$ is the projection, onto the $p$ axis, of the confidence region in the $(C, p)$ plane. It includes not only the samples within the confidence region but all of those in the band $p_{\min} \leq p \leq p_{\max},$ increasing considerably the confidence level.

plot(size=(800, 400), allplts[1].errors, allplts[1].cp)

plot(size=(800, 400), allplts[2].errors, allplts[2].cp)

plot(size=(800, 400), allplts[3].errors, allplts[3].cp)

The plot for the errors only displayed the first two errors. We may also plot a 3d scatter plot of three of the strong errors, i.e. corresponding to three of the four meshes, as, for example,

Histrogram of the order of convergence

Finally, we plot the histograms for $p$, obtained both from the correlated and the decorrelated strong errors. Notice that the distributions for $p$ resembles a normal distribution even for low samples, and building a CI from the decorrelated samples works fine, in this example, despite the fact that the theory does not guarantee that. But it works only with the uncorrelad samples! Nevertheless, it requires a lot more samples, being computationally quite expensive, especially with more complicate equations. The CI from the push-forward method underestimates the confidence level, but it is more trustworthy and less demanding.

plot(allplts[1].histp)

plot(allplts[2].histp)

plot(allplts[3].histp)

This page was generated using Literate.jl.