Add uniform_thinning!

docs/src/library.md

@@ -31,6 +31,7 @@ Pages = ["/wiener.jl"]
 ## Stochastic differential equations
 
 ```@docs
+Bridge.a
 sample
 sample!
 quvar
@@ -54,7 +55,8 @@ GammaProcess
 GammaBridge
 Bridge.ExpCounting
 Bridge.CompoundPoisson
-Bridge.nu 
+Bridge.nu
+Bridge.uniform_thinning!
 ```
 
 ## Miscellaneous

src/levy.jl

@@ -77,8 +77,7 @@ function uniform_thinning!(X, P::GammaProcess, γᵒ)
     y = yy[1]
     for i in 2:length(tt)
         dt = tt[i] - tt[i-1]
-        y = yy[i]
-        yy[i] = y + (yy[i] - y)*rand(Beta(dt* γᵒ,  dt*γ))
+        y, yy[i] = yy[i], yy[i-1] + (yy[i] - y)*rand(Beta(dt*γ,  dt*γᵒ))
     end
     X
 end

test/gamma.jl

@@ -1,11 +1,11 @@
 using Base.Test, Bridge
 
 # test the tricky bridge sampling when endpoint is or is not in grid
-srand(10)
+#srand(10)
 G =  GammaProcess(10.0, 1.5)
 GB = GammaBridge(1.0, 2.0, G)
 n = 1000
-@test abs(mean(sample([0.0, 1.0, 3.0],  G).yy[end] for i in 1:1000) -  mean(Bridge.increment(3, G))) < 3*std(Bridge.increment(0.5, G))/sqrt(n)
+@test abs(mean(sample([0.0, 1.0, 3.0],  G).yy[end] for i in 1:1000) -  mean(Bridge.increment(3, G))) < 4*std(Bridge.increment(0.5, G))/sqrt(n)
 
 X1 = sample([0.0, 0.5, 1.0], GB, 0.2)
 X2 = sample([0.0, 0.5, 2.0], GB, 0.2)
@@ -22,11 +22,15 @@ X3 = sample([0.0, 0.5], GB, 0.2)
 @test  X3.yy[2] <= 2.0
 
 γᵒ, γ = 0.75, 1.2
-la = 0.5
-tt = linspace(0, γ, 1000)
-P = GammaProcess(γ, la)
+t = 10000.
+α = 16.0
+n = 10000
+tt = linspace(0, t, n)
+P = GammaProcess(γ, α)
 X = sample(tt, P)
-var(diff(X.yy))
 
 Y = copy(X)
-Bridge.uniform_thinning!(Y, P, γᵒ)
+Bridge.uniform_thinning!(Y, P, γᵒ)
+
+@test abs(γ*t/n/α - mean(diff(X.yy))) < 0.2/sqrt(n)
+@test abs(γᵒ*t/n/α - mean(diff(Y.yy))) < 0.2/sqrt(n)

test/runtests.jl

@@ -11,6 +11,6 @@ include("guip.jl")
 include("linpro.jl")
 include("timechange.jl")
 include("uniformscaling.jl")
-include("gamma.jl")
+include("gamma.jl") 
 include("gaussian.jl")
-include("with_srand.jl") # run last
+include("with_srand.jl") # run last