Git Product home page Git Product logo

mschauer / bridge.jl Goto Github PK

View Code? Open in Web Editor NEW
111.0 10.0 19.0 10.44 MB

A statistical toolbox for diffusion processes and stochastic differential equations. Named after the Brownian Bridge.

License: Other

Julia 15.01% R 0.34% Jupyter Notebook 84.65%
diffusion julia stochastic-differential-equations mcmc bayesian-inference sde brownian-motion diffusion-processes simulating-diffusion-bridges gamma-process

bridge.jl's Introduction

Build Status Latest Latest

Logo

Bridge.jl

Statistics and stochastic calculus for Markov processes in continuous time, include univariate and multivariate stochastic processes such as stochastic differential equations or diffusions (SDE's) or Levy processes.

I am personally interested in doing Bayesian inference on discretely observed diffusion processes, but this package is written to be of general use and contributions are welcome. Specifically for our code for parameter inference for diffusion processes from discrete data or passage times, check out the dependent package BridgeSDEInference.jl. The statistical method relies a lot on simulating conditional diffusions (so called "difffusion bridges"). See ./example/tutorial.jl for a more general introduction into working with this package.

  • Define and simulate diffusion processes in one or more dimension
  • Continuous and discrete likelihood using Girsanovs theorem and transition densities
  • Monte Carlo sample diffusion bridges, diffusion processes conditioned to hit a point v at a prescribed time T
  • Brownian motion in one and more dimensions
  • Ornstein-Uhlenbeck processes and Ornstein-Uhlenbeck bridges
  • Bessel processes
  • Gamma processes
  • Inhomogenous poisson process
  • Basic stochastic calculus functionality (Ito integral, quadratic variation)
  • Euler-Scheme and implicit methods (Runge-Kutta)
  • Levy-driven SDEs
  • Continuous-discrete filtering for partially observed diffusion processes

The layout/api was originally written to be compatible with Simon Danisch's package FixedSizeArrays.jl. It was refactored to be compatible with StaticArrays.jl by Dan Getz. Some SDE and ODE solvers in Bridge are accessible with the JuliaDiffEq common interface via BridgeDiffEq.jl.

The example programs in the example/ directory have additional dependencies: ConjugatePriors and a plotting library.

Introduction

The key objects introduced are the abstract type ContinuousTimeProcess{T} parametrised by the state space of the path, for example T == Float64 and various structs suptyping it, for example Wiener{Float64} for a real Brownian motion. These play roughly a similar role as types subtyping Distribution in the Distributions.jl package.

Secondly, the struct

struct SamplePath{T}
    tt::Vector{Float64}
    yy::Vector{T}
    SamplePath{T}(tt, yy) where {T} = new(tt, yy)
end

serves as container for sample path returned by direct and approximate samplers (sample, euler, ...). tt is the vector of the grid points of the simulation and yy the corresponding vector of states.

Help is available at the REPL:

help?> GammaProcess
search: GammaProcess LocalGammaProcess VarianceGammaProcess
GammaProcess

A GammaProcess with jump rate γ and inverse jump size λ has increments Gamma(t*γ, 1/λ) and Levy measure

ν(x) = γ x⁻¹ exp(-λ x),

Here Gamma(α,θ) is the Gamma distribution in julia's parametrization with shape parameter α and scale θ.

Examples

julia> sample(linspace(0.0, 1.0), GammaProcess(1.0, 0.5))

Pre-defined processes defined are Wiener, WienerBridge, Gamma, LinPro (linear diffusion/generalized Ornstein-Uhlenbeck) and others.

It is also quite transparent how to add a new process:

using Bridge
using Plots

# Define a diffusion process
struct OrnsteinUhlenbeck  <: ContinuousTimeProcess{Float64}
    β::Float64 # drift parameter (also known as inverse relaxation time)
    σ::Float64 # diffusion parameter
end

# define drift and diffusion coefficient of OrnsteinUhlenbeck
Bridge.b(t, x, P::OrnsteinUhlenbeck) = -P.β*x
Bridge.σ(t, x, P::OrnsteinUhlenbeck) = P.σ

# simulate OrnsteinUhlenbeck using Euler scheme
W = sample(0:0.01:10, Wiener())
X = solve(EulerMaruyama(), 0.1, W, OrnsteinUhlenbeck(2.0, 1.0))
plot(X, label="X")

OrnsteinUhlenbeck

# Levy (Difference-Gamma process) driven OrnsteinUhlenbeck
Z = sample(0:0.01:10, GammaProcess(100.0,10.0))
Z.yy .-= sample(0:0.01:10, GammaProcess(100.0,10.0)).yy
Y = solve(EulerMaruyama(), 0.1, Z, OrnsteinUhlenbeck(2.0, 1.0))
plot(Y, label="Y")

Levy OrnsteinUhlenbeck

Feedback and Contributing

See the documentation for more functionality and issue #12 (Feedback and Contribution) for coordination of the development. Bridge is free software under the MIT licence. If you use Bridge.jl in a closed environment I’d be happy to hear about your use case in a mail to [email protected] and able to give some support.

Literature

F. v. d. Meulen, M. Schauer: Bayesian estimation of discretely observed multi-dimensional diffusion processes using guided proposals. Electronic Journal of Statistics 11 (1), 2017, doi:10.1214/17-EJS1290.

M. Schauer, F. v. d. Meulen, H. v. Zanten: Guided proposals for simulating multi-dimensional diffusion bridges. Bernoulli 23 (4A), 2017, doi:10.3150/16-BEJ833.

bridge.jl's People

Contributors

femtocleaner[bot] avatar getzdan avatar github-actions[bot] avatar jeffbezanson avatar macorstanje avatar mortenpi avatar mschauer avatar simondanisch avatar tkelman avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

bridge.jl's Issues

SDE solver interface overhaul

I plan to bring some oder into the zoo of integrators by introducing the following interface. Currently I am thinking of

For ODEs

solve!(method::ODESolver, X::SamplePath, x0, F::Function)
solve(method::ODESolver, x0, F::Function)

where for example method=BS3() and for SDEs

solve!(method, X::SamplePath, x0, W, P)
solve!(method, X::SamplePath, x0, W, b, σ)
solve(method, x0, W, P)
solve(method, x0, W, b, σ)
bridge!(method, X::SamplePath, W, P)
bridge(method, W, P)

where for example method=Euler() which should cover all cases needed in Bridge.

`bi` interface

I am getting rid of the bi(i, x, P) variants of the generic drift function b(t, x, P) which take an index i instead of a time t. In future SDE solver will call
b((i,t), x, P) with fallback b((i,t), x, P) = b(t, x, P) and models which need the corresponding index and not only the time point t to compute the drift should use dispatch to overwrite that fallback.

Failed to precompile Bridge

I am getting the following error when trying to precompile Bridge (Julia v1.1 on Mac)

[ Info: Precompiling Bridge [2d3116d5-4b8f-5680-861c-71f149790274]
ERROR: LoadError: LoadError: UndefVarError: Zip2 not defined
Stacktrace:
 [1] getproperty(::Module, ::Symbol) at ./sysimg.jl:13
 [2] top-level scope at none:0
 [3] include at ./boot.jl:326 [inlined]
 [4] include_relative(::Module, ::String) at ./loading.jl:1038
 [5] include at ./sysimg.jl:29 [inlined]
 [6] include(::String) at /Users/antonio/.julia/packages/Bridge/s3VsI/src/Bridge.jl:1
 [7] top-level scope at none:0
 [8] include at ./boot.jl:326 [inlined]
 [9] include_relative(::Module, ::String) at ./loading.jl:1038
 [10] include(::Module, ::String) at ./sysimg.jl:29
 [11] top-level scope at none:2
 [12] eval at ./boot.jl:328 [inlined]
 [13] eval(::Expr) at ./client.jl:404
 [14] top-level scope at ./none:3
in expression starting at /Users/antonio/.julia/packages/Bridge/s3VsI/src/types.jl:127
in expression starting at /Users/antonio/.julia/packages/Bridge/s3VsI/src/Bridge.jl:101
ERROR: Failed to precompile Bridge [2d3116d5-4b8f-5680-861c-71f149790274] to /Users/antonio/.julia/compiled/v1.1/Bridge/cfWTr.ji.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1197
 [3] _require(::Base.PkgId) at ./loading.jl:960
 [4] require(::Base.PkgId) at ./loading.jl:858
 [5] require(::Module, ::Symbol) at ./loading.jl:853```

Initial release

  • Travis activated
  • Attobot activated
  • .travis.yml
  • LICENSE.md README.md REQUIRE
  • Tag Release

Address t. kelmans comments

  • various instances of 'type piracy', e.g. chol
  • remove randn(T) is defined in base
  • propose logdet(::Number) for base and add to compat
  • dot is defined for numbers in base
  • avoid dubious dot product for linear scalings
  • use ellipsisnotation.jl
  • remove Compat dependency

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

Feedback and Contributing

This issue is held open for general feedback, feature requests and to coordinate contributions to the package; or just to say "hello".

3DExamples a few issues with upgrade to Julia 1.2

Trying to run the examples on Julia 1.2, I encouter linspace deprecated and replaced by range and some problem with SArray here:

S = 9.0:0.05:11.0
Ρ = 26.0:0.05:30.0
Β = 2.0:0.02:4.0
θref = s0, Ρ[end÷2], Β[end÷2]
@show θref;

θref = (10.0, 27.95, 2.98)

llsurface = [loglikelihood((s0, ρ, β), θref, X) for ρ in Ρ, β in Β];

gives:

The size of type SArray{Tuple{S},T,1,S} where T where S is not known.

If you were trying to construct (or convert to) a StaticArray you
may need to add the size explicitly as a type parameter so its size is
inferrable to the Julia compiler (or performance would be terrible). For
example, you might try

m = zeros(3,3)
SMatrix(m)      # this error
SMatrix{3,3}(m) # correct - size is inferrable

Stacktrace:
[1] error(::String) at .\error.jl:33
[2] missing_size_error(::Type{SArray{Tuple{S},T,1,S} where T where S}) at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\traits.jl:72
[3] Size(::Type{SArray{Tuple{S},T,1,S} where T where S}) at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\traits.jl:88
[4] length(::Type{SArray{Tuple{S},T,1,S} where T where S}) at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\abstractarray.jl:2
[5] convert at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\convert.jl:22 [inlined]
[6] Type at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\convert.jl:7 [inlined]
[7] Diagonal{T,SArray{Tuple{N},T,1,N}} where T where N(::Diagonal{Float64,SArray{Tuple{3},Float64,1,3}}) at C:\Users\Denis.julia\packages\StaticArrays\DBECI\src\SDiagonal.jl:7
[8] σ at C:\Users\Denis.julia\packages\Bridge\vPofv\src\Models.jl:57 [inlined]
[9] a(::Float64, ::SArray{Tuple{3},Float64,1,3}, ::Lorenz) at C:\Users\Denis.julia\packages\Bridge\vPofv\src\types.jl:32
[10] Γ(::Float64, ::SArray{Tuple{3},Float64,1,3}, ::Lorenz) at C:\Users\Denis.julia\packages\Bridge\vPofv\src\types.jl:33
[11] girsanov(::SamplePath{SArray{Tuple{3},Float64,1,3}}, ::Lorenz, ::Lorenz) at C:\Users\Denis.julia\packages\Bridge\vPofv\src\diffusion.jl:119
[12] loglikelihood(::Tuple{Float64,Float64,Float64}, ::Tuple{Float64,Float64,Float64}, ::SamplePath{SArray{Tuple{3},Float64,1,3}}) at .\In[36]:4
[13] (::getfield(Main, Symbol("##17#18")))(::Tuple{Float64,Float64}) at .\none:0
[14] iterate at .\generator.jl:47 [inlined]
[15] collect(::Base.Generator{Base.Iterators.ProductIterator{Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}},getfield(Main, Symbol("##17#18"))}) at .\array.jl:606
[16] top-level scope at In[38]:1

Implement OU bridge

@ChrisRackauckas

Real OU bridge was implemented in 7397f01

For the vector valued OU bridge some matrix algebra is missing, but the one-dimensional formulas give a good idea how the answer should look like.

Bridge.expint inappropriately returns NaN

For exponential integrals with n>1, expint(n, 0.0) returns NaN, instead of 1/(n-1). Compare the behavior to scipy.special.expn, which returns the appropriate value.

Usefulness of LabelledArrays

Hi, Thanks for all your work in making Bridge.jl available.
I'm new to Julia, and feeling my way in terms of leading/best-of-breed packages and idioms.

It seems to me that Bridge.jl is the best of breed diffusion simulating package - closely followed by the DiffEqNoiseProcess.jl, to my mind the BridgeDiffEq,jl supports that conclusion.

I am planning a package that likely would take as input an array such as that generated by Bridge.jl.
In that context I hoped you might share your thoughts on the usefulness of adopting JuliaDiffEq/LabelledArrays.jl.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.