hunting bugs in prediction

Author: smidl

docs/src/lecture_13/ode.jl

@@ -48,6 +48,40 @@ plot!(MXe[2,1:30:end],label="y",color=:red,errorbar=SXe[2,1:30:end])
 
 
 using LinearAlgebra
+function solve_res(f,x0::AbstractVector,sqΣ0, θ,dt,N,Nr)
+   n = length(x0)
+   n2 = 2*length(x0)
+   Qp = sqrt(n)*[I(n) -I(n)]
+
+  X = hcat([zero(x0) for i=1:N]...)
+  S = hcat([zero(x0) for i=1:N]...)
+  X[:,1]=x0
+  Xp = x0 .+ sqΣ0*Qp
+  sqΣ = sqΣ0
+  Σ = sqΣ* sqΣ'
+  S[:,1]= diag(Σ)
+  for t=1:N-1
+    if rem(t,Nr)==0
+      Xp .= X[:,t] .+ sqΣ * Qp
+    end
+    for i=1:n2 # all quadrature points
+      Xp[:,i].=Xp[:,i] + [dt*f(Xp[1:2,i],[Xp[3,i];θ[2:end]]);0]
+    end
+    mXp=mean(Xp,dims=2)
+    X[:,t+1]=mXp
+    Σ=Matrix((Xp.-mXp)*(Xp.-mXp)'/n2)
+    S2=cov(Xp,dims=2)
+    @show Σ
+    @show S2
+    S[:,t+1]=sqrt.(diag(Σ))
+    # @show Σ
+
+    sqΣ = cholesky(Σ).L
+
+  end
+  X,S
+end
+
 function solve(f,x0::AbstractVector,Σ0, θ,dt,N)
    n = length(x0)
    n2 = 2*length(x0)
@@ -69,7 +103,7 @@ function solve(f,x0::AbstractVector,Σ0, θ,dt,N)
     end
     mXp=mean(Xp,dims=2)
     X[:,t+1]=mXp
-    Σ=Matrix((Xp.-mXp)*(Xp.-mXp)'/n2)
+    Σ=cov(Xp,dims=2)
     S[:,t+1]=sqrt.(diag(Σ))
     # @show Σ
   end
@@ -84,23 +118,41 @@ plot!(QX[2,1:30:end],label="y",color=:red,errorbar=QS[2,1:30:end])
 
 savefig("LV_Quadrics.svg")
 
+QXr,QSr=solve_res(f,[1.0,1.0,0.1],diagm([0.1,0.1,0.01]),θ0,0.1,1000,100)
+plot(QXr[1,1:30:end],label="x",color=:blue,errorbar=QSr[1,1:30:end])
+plot!(QXr[2,1:30:end],label="y",color=:red,errorbar=QSr[2,1:30:end])
+
 function filter(f,x0::AbstractVector,Σ0, θ,dt,Ne,Y,σy)
   XX=[]
   SS=[]
   Σ=Σ0
   x=x0
   for t=1:length(Y)
     @show x
-    @show eigen(Σ)
+    @show diag(Σ)
     Xt,St,Xp=solve(f,x,Σ,θ,dt,Ne) # prediction
-    Yp = Xp[1,:] # measure only the first variable
-    mYp = mean(Yp)
-    mXp = mean(Xp,dims=2)
-    SYp = cov(Yp)+σy
-    C = mean([(Xp[:,i].-mXp)*(Yp[i].-mYp) for i=1:6])
-    G = C*inv(SYp)
-    x = vec(Xt[:,end] + G*(Y[t]-mYp))
-    Σ = cov(Xp,dims=2) - G*SYp*G'
+    if false
+      @show Xp
+      Yp = Xp[1,:] # measure only the first variable
+      @show Yp
+      mYp = mean(Yp)
+      mXp = mean(Xp,dims=2)
+      SYp = cov(Yp)+σy
+      @show SYp
+      C = mean([(Xp[:,i].-mXp)*(Yp[i].-mYp) for i=1:6])
+      @show C
+      G = C*inv(SYp)
+      @show G
+      x = vec(Xt[:,end] + G*(Y[t]-mYp))
+      Σ = Matrix((Xp.-mXp)*(Xp.-mXp)'/n2) - G*SYp*G'
+
+      # Σ = cov(Xp,dims=2) - G*SYp*G'
+    else
+        x = vec(mean(Xp,dims=2))
+        Σ = Matrix((Xp.-x)*(Xp.-x)'/6.0)
+      
+    end
+
     push!(XX,Xt)
     push!(SS,St)
   end
@@ -110,11 +162,11 @@ end
 
 
 Y = X[1,100:100:end]
-Xh,Sh=filter(f,[1.0,1.0,0.1],diagm([0.1,0.1,0.01].^2),θ0,0.1,100,Y,0.01)
+Xh,Sh=filter(f,[1.0,1.0,0.1],diagm([0.1,0.1,0.01].^2),θ0,0.1,100,Y[1:end],1e8)
 XF=hcat(Xh...)
 SF=hcat(Sh...)
 
-step=10
+step=1
 plot([1:step:size(XF,2)],XF[1,1:step:end],label="x",color=:blue,errorbar=SF[1,1:step:end])
 plot!([1:step:size(XF,2)],XF[2,1:step:end],label="y",color=:red,errorbar=SF[2,1:step:end])
 scatter!([100:100:size(XF,2)],Y,label="measured",marker=:xcross)