Code Monkey home page Code Monkey logo

kissmcmc.jl's Introduction

KissMCMC

Build Status Build Status Coverage

Got a probability density function you want to draw samples from? Don't want to learn all the fancy stuff of the fancy sampler packages? The KissMCMC (Keep it simple, stupid, MCMC) package intends to provide a few simple MCMC samplers.

using KissMCMC
# the distribution to sample from,
logpdf(x::T) where {T} = x<0 ? -convert(T,Inf) : -x
# initial point of walker
theta0 = 0.5

# Metropolis MCMC sampler:
sample_prop_normal(theta) = 1.5*randn() + theta # samples the proposal (or jump) distribution
thetas, accept_ratio = metropolis(logpdf, sample_prop_normal, theta0, niter=10^5)
println("Accept ratio Metropolis: $accept_ratio")

# emcee MCMC sampler:
thetase, accept_ratioe = emcee(logpdf, make_theta0s(theta0, 0.1, logpdf, 100), niter=10^5)
# check convergence using integrated autocorrelation
thetase, accept_ratioe = squash_walkers(thetase, accept_ratioe) # puts all walkers into one
println("Accept ratio emcee: $accept_ratio")

using Plots
histogram(thetas, normalize=true, fillalpha=0.4)
histogram!(thetase, normalize=true, fillalpha=0.1)
plot!(0:0.01:5, map(x->exp(logpdf(x)[1]), 0:0.01:5), lw=3)

outputs:

MCMC samplers:

  • Metropolis (serial) metropolis
  • Affine invariant MCMC, aka emcee emcee (threaded)

References

Other, probably better Julia MCMC packages:

The (original) emcee python package: https://github.com/dfm/emcee

kissmcmc.jl's People

Contributors

github-actions[bot] avatar juliatagbot avatar mauro3 avatar sefffal avatar

Stargazers

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

Watchers

avatar avatar avatar avatar

Forkers

sefffal standardgalactic szabo137

kissmcmc.jl's Issues

Test sometimes fail

It seems the test:

metropolis: Test Failed at /home/pkgeval/.julia/packages/KissMCMC/z6iaT/test/runtests.jl:37
  Expression: all(abs.(std(thetas) - std_) .< abs.(std_ * tol))
Stacktrace:
 [1] test_mean_std(::Array{Array{Float64,1},1}, ::ATest, ::Float64) at /home/pkgeval/.julia/packages/KissMCMC/z6iaT/test/runtests.jl:37
 [2] macro expansion at /home/pkgeval/.julia/packages/KissMCMC/z6iaT/test/metro.jl:16 [inlined]
 [3] macro expansion at /workspace/srcdir/usr/share/julia/stdlib/v1.5/Test/src/Test.jl:1116 [inlined]
 [4] top-level scope at /home/pkgeval/.julia/packages/KissMCMC/z6iaT/test/metro.jl:3
Test Summary: | Pass  Fail  Total
metropolis    |   27     1     28

sometimes randomly fail (perhaps the tolerance is too tight or something).

Crash with BigFloats

I have a crash which I've narrowed down to the following... no idea whether it's this package's fault though:

using LinearAlgebra, KissMCMC

f(θs, ts) = [-t*exp(θ-exp(θ)*t)/length(θs) for t in ts, θ in θs]

g(v, ts) = begin
    j = f(v, ts)
    x = j' * j
    # s, ld = LinearAlgebra.logabsdet(x)
    # ifelse(s>0, Float64(ld), -Inf)
    d = det(x)
    d>0 ? Float64(log(d)) : -Inf
end

KissMCMC.emcee(v -> g(v, 1:3), [rand(3) for _ in 1:10]; niter=10^4) # fine!

versioninfo()
KissMCMC.emcee(v -> g(v, big.(1:3.0)), [big.(rand(3)) for _ in 1:10]; niter=10^4)
julia> versioninfo()
Julia Version 1.5.0
Commit 96786e22cc (2020-08-01 23:44 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i5-7360U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 4
  JULIA_PKG_SERVER = https://geo.pkg.julialang.org

julia> KissMCMC.emcee(v -> g(v, big.(1:3.0)), [big.(rand(3)) for _ in 1:10]; niter=10^4)
julia(303,0x70000e9e8000) malloc: *** error for object 0x7fb6d3ea8d70: pointer being freed was not allocated
julia(303,0x108ff35c0) malloc: *** error for object 0x7fb6d306b460: pointer being freed was not allocated
julia(303,0x70000e9e8000) malloc: *** set a breakpoint in malloc_error_break to debug
julia(303,0x108ff35c0) malloc: *** set a breakpoint in malloc_error_break to debug

signal (6): Abort trap: 6
in expression starting at REPL[5]:1
in expression starting at Rin expression starting at REPL[5]:1
__pthread_kill at /usr/lib/system/libsystem_kernel.dylib (unknown line)
Allocations: 34375251 (Pool: 34366235; Big: 9016); GC: 38
Abort trap: 6

(And the same on Julia 1.4.2)

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.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Image storage

This issue can be used to store images, used in, e.g. the README.

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.