eloceanografo / marginallogdensities.jl Goto Github PK
View Code? Open in Web Editor NEWMarginalized log-probability functions in Julia
License: MIT License
Marginalized log-probability functions in Julia
License: MIT License
The sparse Hessian coloring/decompression code in this package has been improved and incorporated into SparseDiffTools.jl (JuliaDiff/SparseDiffTools.jl#190). Should replace it with that implementation.
SparsityDetection.jl is no longer maintained, and the sparsity detection functionality in Symbolics.jl still seems somewhat brittle. Currently this package just uses ForwardDiff, but that is a suboptimal solution (very slow in high dimensions).
Hi! Why not in regisry?
Hi @ElOceanografo, thanks for this very nice (albeit in-progress) package!
I am trying to implement the "urchins" example from Simon Wood's "Core Statistics" to test the package on a somewhat less typical example (see here for the example in R).
I've implemented it as below. When calling the callable MarginalLogDensity
with both the random and fixed effects as input, it returns the log likelihood correctly. However, when calling it with only the fixed effects to marginalize over the random effects I get an error about dimension mismatch. I tracked it down to this particular line in the source code https://github.com/ElOceanografo/MarginalLogDensities.jl/blob/master/src/MarginalLogDensities.jl#L241 but I'm not sure where in my code I've implemented things incorrectly to get such an error. If you are able to take a look, I'd greatly appreciate it.
using Plots
using MarginalLogDensities
using Distributions
using Random
using Optim
function urchin_V(ω, p, g, a)
ω = exp(ω)
p = exp(p)
g = exp(g)
am = log(p/(g*ω))/g
if a < am
return ω*exp(g*a)
else
return p/g + p*(a-am)
end
end
# number of data
const n = 100
# true parameters
# fixed effects
θ_true = zeros(206)
θ_true[1] = -4.0
θ_true[2] = -0.5
θ_true[3] = log(0.2)
θ_true[4] = 1
θ_true[5] = log(0.1)
θ_true[6] = log(0.1)
# random effects
θ_true[7:106] = rand(Normal(θ_true[2], exp(θ_true[3])), n) #g
θ_true[107:206] = rand(Normal(θ_true[4], exp(θ_true[5])), n) #p
# simulate true data
# urchin ages
a = Float64.(sample(1:30, n))
# urchin volumes
v = urchin_V.(repeat([θ_true[1]],n), θ_true[107:206], θ_true[7:106], a)
v_samp = rand.(Normal.(sqrt.(v), exp(θ_true[6]))) # data + measurement error for likelihood computation
# plot the "data"
scatter(a,v,legend=false)
scatter(a,v_samp,legend=false)
# likelihood function for model with fixed and random effects
function loglik_urchin(θ::Vector{T}) where {T<:Real}
# parameters
log_ω = θ[1]
log_g = θ[7:106]
log_p = θ[107:206]
μ_g = θ[2]
log_σ_g = θ[3]
μ_p = θ[4]
log_σ_p = θ[5]
log_σ = θ[6]
# estimated volumes of urchins conditional on current θ
v_est = urchin_V.(repeat([log_ω],n), log_p, log_g, a)
# data likelihood
data_lik = sum(logpdf.(Normal.(sqrt.(v_est), exp(log_σ)), v_samp))
# random effect (growth rates) likelihood
g_lik = loglikelihood(Normal(μ_g, exp(log_σ_g)), log_g)
p_lik = loglikelihood(Normal(μ_p, exp(log_σ_p)), log_p)
return data_lik + g_lik + p_lik
end
loglik_urchin(θ_true)
n_θ = length(θ_true)
marginal_ix = collect(7:206)
marginalloglik_urchin = MarginalLogDensity(loglik_urchin, n_θ, marginal_ix)
# check it works at the true values
marginalloglik_urchin(θ_true[7:end], θ_true[1:6])
# just fixed effects
marginalloglik_urchin(θ_true[1:6])
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!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.