Added Banana shaped distribution example

[burtonjosh]

as mentioned in #10, this PR adds an example of using the package on a 2D banana shaped distribution.

[mschauer]

If this is interesting, for context https://discourse.julialang.org/t/mcmc-benchmark/63013

[burtonjosh]

I like @mschauer's suggestion about the n-dimensional Rosenbrock distribution. Maybe it is a little complicated/involved for the quickstart page? I have an implementation for this if you are interested to see. Not sure if I should put it here or make a new PR?

[sethaxen]

I like @mschauer's suggestion about the n-dimensional Rosenbrock distribution. Maybe it is a little complicated/involved for the quickstart page? I have an implementation for this if you are interested to see. Not sure if I should put it here or make a new PR?

It's a good suggestion. I agree, it would be nice to have. RE the Quickstart, I think we only want to include it if it's not too complicated and multi-path Pathfinder more-or-less (or almost) works (we already have the funnel for which that isn't the case). Since you have an implementation, could you generate the plot of draws and the animation of fitting and drop them here in a comment? Then we can discuss whether we think they should go into the quickstart.

[burtonjosh]

here is the animation,

and here are two of the marginals for the multipath plot (the contours are done by sampling directly then plotting a kde)

[sethaxen]

Thanks! Hm, yeah these don't look terribly interesting. I'm surprised at how poor the samples look. It looks like it's initializing very close to the mode, while the scale is quite large. Could you use init_scale=10 perhaps, and rerun? Also, could you share the code so I can check a few things?

[burtonjosh]

Sure!

using LinearAlgebra, Pathfinder, Printf, StatsPlots, Random, Distributions, KernelDensity
Random.seed!(23)

struct HybridRosenbrock{Tμ,Ta,Tb} <: ContinuousMultivariateDistribution
    μ::Tμ
    a::Ta
    b::Tb
end

function Distributions.logpdf(D::HybridRosenbrock, x::AbstractVector{})
    n2, n1 = size(D.b)
    n = length(x)
    X = reshape(@view(x[2:end]), n1, n2)
    y = -D.a*(x[1] - D.μ)^2 
    for j in 1:n2 
        y -= D.b[j,1]*(X[j,1] - x[1]^2)^2
        for i in 2:n1
            y -= D.b[j,i]*(X[j,i] - X[j,i-1]^2)^2
        end
    end
    C = 0.5log(D.a) + 0.5sum(log.(D.b)) - n/2*log(pi)
    y + C
end

function Distributions.rand(rng::Random.AbstractRNG, D::HybridRosenbrock)
    x = zeros(length(D.b) + 1)
    n2, n1 = size(D.b)
    X = reshape(@view(x[2:end]), n2, n1)
    x[1] = D.μ + randn(rng)/(2*D.a)
    for j in 1:n2 
        X[j,1] = x[1]^2 + randn(rng)/(2D.b[j,1])
        for i in 2:n1
            X[j,i] = X[j,i-1]^2 + randn(rng)/(2D.b[j,i])
        end
    end
    x
end

μ = 1
a = 1/20
b = 100/20*ones(2,2)
n2, n1 = size(b)

function logp_5drosenbrock(x)
    n = length(x)
    X = reshape(@view(x[2:end]), n1, n2)
    y = -a*(x[1] - μ)^2 
    for j in 1:n2 
        y -= b[j,1]*(X[j,1] - x[1]^2)^2
        for i in 2:n1
            y -= b[j,i]*(X[j,i] - X[j,i-1]^2)^2
        end
    end
    C = 0.5log(a) + 0.5sum(log.(b)) - n/2*log(pi)
    y + C
end

result = pathfinder(logp_5drosenbrock; dim=5, init_scale=10)

iterations = length(result.optim_trace) - 1
trace_points = result.optim_trace.points
trace_dists = result.fit_distributions

xrange = -20:20
yrange = -45:350

# sample directly and replace contour with kde of the samples
D = HybridRosenbrock(μ, a, b)
x_direct = hcat([rand(D) for k in 1:20000]...)'
x_contour = kde(x_direct[:,1:2])

anim = @animate for i in 1:iterations
    plot(x_contour,label="")
    trace = trace_points[1:(i + 1)]
    dist = trace_dists[i + 1]
    plot!(first.(trace), last.(trace); seriestype=:scatterpath, color=:black, msw=0, label="trace")
    covellipse!(dist.μ[1:2], dist.Σ[1:2, 1:2]; n_std=2.45, alpha=0.7, color=1, linecolor=1, label="MvNormal 95% ellipsoid")
    title = "Iteration $i"
    plot!(; xlims=extrema(xrange), ylims=extrema(yrange), xlabel="x₁", ylabel="x₂", legend=:top, title)
end
gif(anim, fps=5)

[sethaxen]

Ah, okay, I see what the issue is. The marginal distribution of logp_5drosenbrock is

$$x[1] ~ Normal(μ, inv(2a)) x[2] ~ Normal(x[1]^2, inv(2b[1,1]))$$

If we plot the contours of this density, it looks very different than the KDE, which overaggressively smooths it. It basically looks like a parabola. This is a very hard distribution to fit, and understandably, Pathfinder fails to produce draws over the whole distribution, though it does produce draws from the high-density part. Since we already have an example like this in the quickstart (failure due to high curvature/correlations), I think we should stick with the original banana function.

[sethaxen]

Just a few minor tweaks. Then I think it's good to go.

[codecov]

Codecov Report

Merging #68 (49d8e9b) into main (2723aff) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##             main      #68   +/-   ##
=======================================
  Coverage   92.89%   92.89%           
=======================================
  Files          13       13           
  Lines         521      521           
=======================================
  Hits          484      484           
  Misses         37       37           

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2723aff...49d8e9b. Read the comment docs.

[burtonjosh]

accepted all the changes :)

[sethaxen]

Thanks @burtonjosh for the excellent contribution!