Skip to content

Instantly share code, notes, and snippets.

@tpoisot

tpoisot/turing.jl

Last active February 13, 2020 02:49
Download ZIP

Revisions

  1. tpoisot revised this gist Feb 13, 2020. 1 changed file with 29 additions and 14 deletions.
    43 changes: 29 additions & 14 deletions turing.jl
    View file
    Original file line number Diff line number Diff line change
    @@ -3,20 +3,21 @@ Pkg.activate(".")

    using Turing
    using StatsPlots
    using StatsBase
    import CSV

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

    @model model(x,y) = begin
    @model flexible_links(x,y) = begin
    N = length(x)
    # Transformed data
    F = x.*x .- (x.-1)
    # Parameters
    ϕ ~ Normal(3.0, 0.5)
    ϕ ~ Normal(3.0, 0.5)
    μ ~ Beta(3.0, 7.0)
    for i in 1:N
    y[i] ~ BetaBinomial(F[i], μ*exp(ϕ), (1-μ)*exp(ϕ))
    y[i] ~ BetaBinomial(F[i], μ*exp(ϕ), (1-μ)*exp(ϕ))
    end
    return μ, ϕ
    end
    @@ -25,7 +26,7 @@ s = vec(d[:,:S])
    l = vec(d[:,:L])
    m = l .- (s.-1)

    tchain = mapreduce(c -> sample(model(s,m), HMC(0.01,10), 2000), chainscat, 1:8);
    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,:,:])
    @@ -34,16 +35,30 @@ function prediction(chain, S)
    return rand(BetaBinomial(S^2-(S-1), μ*exp(ϕ), (1-μ)*exp(ϕ)))+(S-1)
    end

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

    scatter(S, L, lab="", alpha=0.05)
    plot!(s, s.-1, lw=1, c=:black, lab="")
    plot!(s, s.^2, lw=2, c=:black, lab="")
    xaxis!(:log, (3,750), "Richness")
    yaxis!(:log, (1,10000), "Links")
    savefig("prediction.png")
    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

    scatter(S, L./(S.^2))
    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")
  2. tpoisot created this gist Feb 12, 2020.
    49 changes: 49 additions & 0 deletions turing.jl
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,49 @@
    import Pkg
    Pkg.activate(".")

    using Turing
    using StatsPlots
    import CSV

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

    @model model(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)

    tchain = mapreduce(c -> sample(model(s,m), HMC(0.01,10), 2000), chainscat, 1:8);

    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:750
    S = rand(3:750, 5000)
    L = [prediction(tchain, s) for s in S]

    scatter(S, L, lab="", alpha=0.05)
    plot!(s, s.-1, lw=1, c=:black, lab="")
    plot!(s, s.^2, lw=2, c=:black, lab="")
    xaxis!(:log, (3,750), "Richness")
    yaxis!(:log, (1,10000), "Links")
    savefig("prediction.png")


    scatter(S, L./(S.^2))