Confidence ellipsoids for SDEs with Bridge and Makie

ellipses.png

ellipsoid.jl

using Makie
using Bridge
using StaticArrays
using Statistics
using LinearAlgebra
using Colors

T = 10
K = 10000
t = 0:0.01:T
stept = 100
stepk = 100
R3 = SVector{3,Float64}
W = [sample(t, Wiener{R3}()) for k in 1:K]

b(t, x) = 5*R3(0.0, 0.0, 1.0)
σ(t, x) = SDiagonal(1.0, 0.4, 2.0)
#for k in 1:K
#    W[k].yy .+= [5* t * R3(0.0, 0.0, 1.0) for t in t]
#end
X = [solve(Euler(), R3(0.0, 0.0, 0.0), W[k], (b, σ)) for k in 1:K]

kith(X, k, i) = X[k].yy[i]
m = mean(Xᵢ.yy[1:stept:end] for Xᵢ in X)
v = [cov([kith(X, k, i) for k in 1:K]) for i in 1:stept:length(t)]
sv = svd.(v)
scal = map(sv -> sqrt.(sv.S), sv)
svu = [sv.U for sv in sv]
rot = Bridge.quaternion.(svu)
rot = map(x-> AbstractPlotting.Quaternionf0(x...), rot)

sphere = Sphere(Point3f0(0,0,0), 1.0f0)

scene = lines(X[1].yy)
for k in 1:stepk:K
    k != 1 && lines!(scene, X[k].yy, linewidth=0.5)
end
display(scene)
meshscatter!(scene, m, markersize = scal, rotations = rot, marker = sphere, color=RGBA(1., 0.0, 0.0, 0.2))
display(scene)