Merge pull request #16 from mschauer/nine

Towards release 0.8

docs/src/library.md

@@ -5,6 +5,7 @@
 ```@docs
 ContinuousTimeProcess{T}
 SamplePath{T}
+Bridge.GSamplePath
 valtype
 Bridge.outertype
 ```
@@ -50,6 +51,13 @@ StochasticHeun
 Bridge.NoDrift
 ```
 
+## In place solvers
+```@docs
+Bridge.R3!
+Bridge.σ!
+Bridge.b!
+```
+
 ## Levy processes
 ```@docs
 GammaProcess
@@ -60,6 +68,12 @@ Bridge.nu
 Bridge.uniform_thinning!
 ```
 
+## Poisson processes
+```@docs
+Bridge.ThinningAlg
+Bridge.InhomogPoisson
+```
+
 ## Miscellaneous
 
 ```@docs
@@ -139,7 +153,7 @@ Bridge.dotVs
 Bridge.SDESolver
 Bridge.Increments
 Bridge.sizedtype
+Bridge.piecewise
+Bridge.BridgePre!
 ```
-
-
-
+

project/gamma3.jl

@@ -24,12 +24,10 @@ if sim == :fire
     #b = [0.4, 1.2, 3.6]
     b = [0.75, 2.5, 3.6]
     b = [0.5, 1.5, 3]
-    b = [0.5, 1, 1.5, 2, 2.5, 3., 3.5]
-    b = [0.75, 2, 4]
-    b = [0.75, 2.5, 3.6]
-    b = [1.0, 2.0, 4.0]
-    #b = [1.0, 1.5, 3.5, 5]
-#    b = [3.0, 5.0]
+    b = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5]
+    b = [0.5, 1.0, 1.5, 2.0, 3.0, 4.0 ]
+
+ #   b = [1.0, 2.0, 4.0]
     N = length(b)
 end
 
@@ -43,7 +41,7 @@ end
 tt = linspace(0.0, T, n + 1)
 
 
-iterations = 200_000
+iterations = [100_000, 200_000][1]
 
 
 # two gamma processes
@@ -83,7 +81,7 @@ elseif sim == :fire
     T = tt_[end]
     tt = tt_
 
-    beta0 = 1/dt*1.812 # ml 
+    beta0 = 1/dt*1.812 # Maximum likelihood estimate from the R script
     alpha0 = 0.588942 
 
 end
@@ -178,11 +176,12 @@ else
     X = SamplePath(tt, Bridge.cumsum0(dxx))
 end
 
+alphahat = mean(dxx)/var(dxx) # use estimate for alpha
 
 # grid points
 # b = quantile(increment(dt, GammaProcess(beta0,alpha0)),(1:(N))/(N+1)) # theoretical
 if simid != 3 && N >= 1
-    b = cumsum(0.3:0.4:2.0)[1:N]
+    b = cumsum(0.5:0.4:2.0)[1:N]
     #b = quantile(diff(X.yy), (N:2N-1)/(2N)) # first bin resembles 50% of emperical increment distributions.
 elseif N == 0
     b = Float64[]
@@ -246,12 +245,14 @@ beps = 0.0
 #var(yy[yy.< b[1]/4])
 
 
-alpha = alpha0
+#alpha = alpha0
+alpha = alphahat
+
 beta = beta0
 
 alpha = alpha
 if transdim
-    beta = beta
+    beta = 1.2*beta
 end
 
 #c = beta*(T/(n*m*alpha)) *(1-exp(-alpha*beps)) # compensator for small jumps
@@ -262,8 +263,10 @@ rho = zeros(N)
 
 # prior chain step std
 alphasigma = 0.15
-thsigma = 0.05
-rhosigma = 0.05*ones(N) # variance of rho1 = 0
+thsigma = 0.02*ones(N)
+#thsigma[1] = 0.02 #Fixme
+
+rhosigma = 0.02*ones(N) # variance of rho1 = 0
 if N > 0 && sim != :fire
     rhosigma[1] = 0.0
 end
@@ -356,9 +359,9 @@ for iter in 1:iterations
 
     # sample parameters
     # update theta and rho
-    if iter % 5 != 2 # remember to update formula for acceptane rates
+    if iter % 5 != 2 # remember to update formula for acceptance rates
         P0 = GammaProcess(beta, alpha)
-        thetaº = theta + thsigma*randn(length(theta))
+        thetaº = theta + thsigma.*randn(length(theta))
         rhoº = rho + rhosigma.*randn(length(rho))
 
         if N == 0 || thetaº[end] + alpha < eps() 
@@ -382,7 +385,7 @@ for iter in 1:iterations
         end
     end
 
-    if transdim && iter % 5 == 4
+    if transdim && iter % 5 == 2
         betaº = beta + betasigma*randn()
         if betaº <= 0
             #reject

src/Bridge.jl

@@ -1,7 +1,7 @@
 __precompile__()
 
 module Bridge
-export ContinuousTimeProcess, SamplePath, VSamplePath
+export ContinuousTimeProcess, SamplePath, VSamplePath, GSamplePath
 export stack
 
 export LinPro, Wiener, WienerBridge, CSpline
@@ -50,6 +50,35 @@ function mcsvd3
 end
 function visualize_uncertainty
 end
+function B!
+end
+function a!
+end
+
+"""
+    b!(t, y, tmp1, P)
+
+Compute drift ``b`` in `y` (without factor ``Δt``, modifying `tmp1`.
+"""
+function b!
+end
+function bi!
+end
+function ri!
+end
+"""
+    σ!(t, y, Δw, tmp2, P)
+
+Compute stochastic increment at `y`, ``σ Δw``, modifying `tmp2`.
+"""
+function σ!
+end
+
+function P
+end
+function Pt
+end
+
 
 hasbi(::Any) = false
 hasai(::Any) = false
@@ -67,10 +96,13 @@ include("cspline.jl")
 include("wiener.jl")
 include("ellipse.jl")
 include("ode.jl")
+include("ode!.jl")
 include("diffusion.jl")
 include("linpro.jl")
 include("guip.jl")
 include("euler.jl")
+
+include("sde!.jl")
 include("levy.jl")
 include("poisson.jl")
 

src/ellipse.jl

@@ -3,7 +3,7 @@ const .. = Val{:...}
 
 setindex!(A::AbstractArray{T,1}, x, ::Type{Val{:...}}, n) where {T} = A[n] = x
 setindex!(A::AbstractArray{T,2}, x, ::Type{Val{:...}}, n) where {T} = A[ :, n] = x
-setindex!(A::AbstractArray{T,3}, x, ::Type{Val{:...}}, n) where {T} = A[ :, :, n] =x
+setindex!(A::AbstractArray{T,3}, x, ::Type{Val{:...}}, n) where {T} = A[ :, :, n] = x
 
 getindex(A::AbstractArray{T,1}, ::Type{Val{:...}}, n) where {T} = A[n]
 getindex(A::AbstractArray{T,2}, ::Type{Val{:...}}, n) where {T} = A[ :, n]

src/euler.jl

@@ -29,6 +29,9 @@ Precomputed Euler-Maruyama scheme for bridges using `bi`.
 struct BridgePre <: SDESolver
 end
 
+
+
+
 """
     StochasticHeun() <: SDESolver
 
@@ -86,19 +89,37 @@ end
 
 
 """
-    solve!(method::SDESolver, Y, u, W::SamplePath, P) -> X
+    solve(method::SDESolver, u, W::SamplePath, P) -> X
+    solve(method::SDESolver, u, W::SamplePath, (b, σ)) -> X
   
 Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t`` 
 using `method` in place.
+
+# Example
+
+```
+solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), ((t,x)->-x, (t,x)->I))
+```
+
+```
+import Bridge: b, σ
+struct OU <: ContinuousTimeProcess{Float64}
+    μ::Float64
+end
+Bridge.b(s, x, P::OU) = -P.μ*x
+Bridge.σ(s, x, P::OU) = I
+
+solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), OU(1.4))
+```
 """
 solve(method::SDESolver, u::T, W::SamplePath, P::ProcessOrCoefficients) where {T} =
     solve!(method, SamplePath{T}(W.tt, T[zero(u) for t in W.tt]), u, W, P)
 
 """
-    solve!(method::SDESolver, Y, u, W::VSamplePath, P) -> X
+    solve(method::SDESolver, u, W::VSamplePath, P) -> X
   
 Solve stochastic differential equation ``dX_t = b(t,X_t)dt + σ(t,X_t)dW_t`` 
-using `method` in place.
+using `method`.
 """
 solve(method::SDESolver, u, W::VSamplePath{T}, P::ProcessOrCoefficients) where {T} =
     solve!(method, VSamplePath(W.tt, zeros(T, size(u)..., length(W.tt))), u, W, P)
@@ -147,38 +168,8 @@ function solve!(::EulerMaruyama, Y, u::T, W::AbstractPath, P::ProcessOrCoefficie
     yy[.., N] = y
     Y
 end
-function solve!(::EulerMaruyama!, Y, u::T, W::AbstractPath, P::ProcessOrCoefficients) where {T}
-    N = length(W)
-    N != length(Y) && error("Y and W differ in length.")
 
-    tt = Y.tt
-    tt[:] = W.tt
-    yy = Y.yy
-    y::T = copy(u)
 
-    assert(size(Y.yy) == (length(y), N))
-    assert(size(W.yy) == (length(y), N))
-    tmp1 = copy(y)
-    tmp2 = copy(y)
-    dw = W.yy[.., 1]
-    for i in 1:N-1
-        t¯ = tt[i]
-        dt = tt[i+1] - t¯ 
-        for k in eachindex(tmp1)
-            @inbounds yy[k, i] = y[k]
-        end
-        for k in eachindex(dw)
-            @inbounds dw[k] = W.yy[k, i+1] - W.yy[k, i]
-        end
-        b!(t¯, y, tmp1, P)
-        σ!(t¯, y, dw, tmp2, P)
-        for k in eachindex(y)
-            @inbounds y[k] = y[k] + tmp1[k]*dt + tmp2[k]
-        end
-    end
-    yy[.., N] = y
-    Y
-end
 
 """
     bridge(method, W, P) -> Y
@@ -218,6 +209,9 @@ function bridge!(::BridgePre, Y, W::SamplePath, P::ContinuousTimeProcess{T}) whe
     Y
 end
 
+
+#### Guided Bridges
+
 bridge(u::T, W::SamplePath, P::GuidedBridge{T}) where {T} = bridge!(samplepath(W.tt, u), W, P)
 """
     bridge!(Y, u, W, P::GuidedBridge) -> v
@@ -249,6 +243,9 @@ function bridge!(Y, u, W::SamplePath, P::GuidedBridge{T}) where {T}
     yy[.., N]
 end
 
+####
+
+
 """
     bridge!(method, Y, W, P) -> Y
 

src/gaussian.jl

@@ -52,7 +52,7 @@ whiten(Σ::UniformScaling, z) = z/sqrt(Σ.λ)
 sqmahal(P::Gaussian, x) = norm_sqr(whiten(P.Σ,x - P.μ))
 
 rand(P::Gaussian) = P.μ + chol(P.Σ)'*randn(typeof(P.μ))
-rand(P::Gaussian{Vector}) = P.μ + chol(P.Σ)'*randn(T, length(P.μ))
+rand(P::Gaussian{Vector}) = P.μ + chol(P.Σ)'*randn(eltype(P.μ), length(P.μ))
 
 logpdf(P::Gaussian, x) = -(sqmahal(P,x) + _logdet(P.Σ, dim(P)) + dim(P)*log(2pi))/2    
 pdf(P::Gaussian, x) = exp(logpdf(P::Gaussian, x))

src/guip.jl

@@ -11,6 +11,11 @@ end
 tau(s, t, T) = t + (s.-t).*(2-(s-t)/(T-t))
 tau(ss::Vector) = tau(ss, ss[1], ss[end])
 
+
+# fallback
+btilde(t, x, Po) = b(t, x, P(Po))
+btilde!(t, x, out, Po) = b!(t, x, out, P(Po))
+
 #####################
 
 
@@ -217,7 +222,7 @@ function gpupdate(H♢, V, L, Σ, v)
         V_ = (L' * inv(Σ) * L)\(L' * inv(Σ) *  v)
         H♢_, V_
     else
-        Z = I - H♢*L'*inv(Σ + L*H♢*L')*L
+        Z = I - H♢*L'*inv(Σ*I + L*H♢*L')*L
         Z*H♢, Z*H♢*L'*inv(Σ)*v + Z*V
     end
 end
@@ -248,6 +253,8 @@ function logdetU(GP1, GP2, L, Σ)
     logdet(inv(K) + L'*inv(Σ)*L + inv(H)) + logdet(Σ) + logdet(H) + logdet(K) + 2logdet(PhiTS)
 end
 
+#################################################
+abstract type GuidedBridge!{T} <: ContinuousTimeProcess{T} end
 
 #################################################