turing.jl

import Pkg
Pkg.activate(".")

using Turing
using StatsPlots
using StatsBase
import CSV

d_orig = CSV.read("ls.csv")
d = d_orig[d_orig.P .> 0,:]

@model flexible_links(x,y) = begin
    N = length(x)
    # Transformed data
    F = x.*x .- (x.-1)
    # Parameters
    ϕ ~ Normal(3.0, 0.5)
    μ ~ Beta(3.0, 7.0)
    for i in 1:N
        y[i] ~ BetaBinomial(F[i], μ*exp(ϕ), (1-μ)*exp(ϕ))
    end
    return μ, ϕ
end

s = vec(d[:,:S])
l = vec(d[:,:L])
m = l .- (s.-1)

fl_chain = mapreduce(c -> sample(flexible_links(s,m), HMC(0.01,10), 2000), chainscat, 1:2);

function prediction(chain, S)
    p = get_params(chain[200:end,:,:])
    i = rand(1:length(p.μ))
    μ, ϕ = p.μ[i], p.ϕ[i]
    return rand(BetaBinomial(S^2-(S-1), μ*exp(ϕ), (1-μ)*exp(ϕ)))+(S-1)
end

s = 3:1:100
S = rand(minimum(s):maximum(s), 5000)
L = [prediction(fl_chain, s) for s in S]

function quantiles(v, q)
    lims = zeros(eltype(v), length(q))
    eq = ecdf(v)(v)
    for (i,s) in enumerate(q)
        lims[i] = v[last(findmin(abs.(eq.-s)))]
    end
    return lims
end

q = [0.07, 0.13, 0.5, 1.0-0.13, 1.0-0.07]
Qs = zeros(Int64, (length(q), length(s)))
@progress for i in eachindex(s)
    Li = [prediction(fl_chain, s[i]) for rep in 1:2000]
    Qs[:,i] = quantiles(Li, q)
end

plot(Qs[5,:], fillrange=Qs[1,:], fillalpha=0.2, c=:lightgrey, lw=0, lc=:transparent, lab="", frame=:box)
plot!(Qs[4,:], fillrange=Qs[2,:], fillalpha=0.2, c=:lightgrey, lw=0, lc=:transparent, lab="")
plot!(Qs[3,:], c=:black, lw=1.5, lab="", ls=:dash)
plot!(s, s.-1, lw=1, c=:black, lab="")
plot!(s, s.^2, lw=2, c=:black, lab="")
xaxis!(:log, (3,100), "Richness")
yaxis!(:log, (1,5000), "Links")