Skip to content

Running my own log-likelihood function #98

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
urirolls5987 opened this issue Jul 27, 2023 · 1 comment
Open

Running my own log-likelihood function #98

urirolls5987 opened this issue Jul 27, 2023 · 1 comment

Comments

@urirolls5987
Copy link

urirolls5987 commented Jul 27, 2023

Hi!
Thank you for an incredible package, this is really quite amazing.
I was faced with an issue when it comes to running my own log-likelihood function, and sampling that similar to how I would do in Dynesty.jl:

loglikelihood(sample_params) = somefunction(sample_params) smplr = NestedSampler(ndim, nlive=500) res = dysample(loglikelihood, identity, smplr; dlogz=0.5)

In Pigeons there is an example on how to do something like this (general_target.jl) but it's unclear how to define the Prior.

For example, using a uniform prior distribution (essentially the identity).

Thanks!

@miguelbiron
Copy link
Collaborator

miguelbiron commented Jul 28, 2023

Hi -- sorry for the delay. I've been working on a PR that exposes a shortcut interface to pigeons in the case where your loglikelihood is given by a blackbox function, and your reference is determined by a Distribution.jl object.

In the meantime, if you want to check it out, you can install it via

]add git@github.com:Julia-Tempering/Pigeons.jl.git#bread-crumbs-api

Note: this interface is still experimental and subject to changes in the future

To replicate the black-box version of the unidentifiable model example in the docs using this simplified interface, simply run

    using Pigeons, Distributions, MCMCChains

    # define the target loglikelihood 
    function unid_log_potential(x; n_trials=100, n_successes=50) 
        p = prod(x)
        return n_successes*log(p) + (n_trials-n_successes)*log1p(-p)
    end
    ref_dist = product_distribution(Uniform(), Uniform()) # define the reference distribution 
    pt = pigeons(
        BreadCrumbs(unid_log_potential, ref_dist),
        n_rounds = 12,
        record = [traces]
    )

    # collect the statistics and convert to MCMCChains' Chains
    samples = Chains(sample_array(pt), variable_names(pt))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants