Skip to content
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

Beta Hierarchy #157

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Beta Hierarchy #157

wants to merge 6 commits into from

Conversation

mberaha
Copy link
Contributor

@mberaha mberaha commented Mar 12, 2024

This PR implements a Beta mixture kernel where the parameters of the Beta distribution are given independent Gamma priors. This allows using models of the kind

$$\begin{align*} y_1, \dots, y_n &\sim \int \text{Beta} ( \cdot \mid \alpha, \gamma) P(\text{d}\alpha \, \text{d}\beta ) \\[3pt] P &\sim DP(\theta G_0) \end{align*}$$

(or any variation including Pitman-Yor process, MFM, truncated stickbreaking..) where

$$G_0(\text{d}\alpha \, \text{d}\beta) = \text{Gamma}(\text{d}\alpha | a_1, b_1) \text{Gamma}(\text{d}\gamma | a_2, b_2)$$

The resulting hierarchy is not conjugate and is sampled using a simple Random Walk Updater, which works well in my experience.

Upon merging, this closes #156.

TODO:

  • Write docstrings for most of the classes
  • Check why some I deleted some of the tests (LOL)

@jdtuck
Copy link

jdtuck commented Mar 12, 2024

This is what I am looking for. You can also use the prior form Kottas (2006) which is

$$U(d\alpha|[0,T])Inv-Gamma(d\gamma|a_0,b_0)$$

It would be interesting to compare as he suggests $a_0=2$ and $b_0=8$. You can also put a prior on $b_0$ of $Gamma(1,0.125)$ as well. I will test it out.

@jdtuck
Copy link

jdtuck commented Mar 12, 2024

Also do you have an example for running this?

@jdtuck
Copy link

jdtuck commented Mar 13, 2024

Just tried to build this PR and does not build

*** Building BayesMix executable ***
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/blocked_gibbs_algorithm.cc.o
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/conditional_algorithm.cc.o
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/base_algorithm.cc.o
[  5%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/marginal_algorithm.cc.o
[  7%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/neal2_algorithm.cc.o
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/blocked_gibbs_algorithm.cc:1:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/blocked_gibbs_algorithm.h:5:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/conditional_algorithm.h:7:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/base_algorithm.h:14:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/base_hierarchy.h:11:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/abstract_hierarchy.h:15:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/priors/abstract_prior_model.h:11:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/likelihoods/states/includes.h:6:
/Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/likelihoods/states/shape_rate_state.h:68:20: error: no member named 'sr_state' in 'bayesmix::AlgorithmState_ClusterState'
    shape = state_.sr_state().shape();
            ~~~~~~ ^

@TeoGiane
Copy link
Contributor

TeoGiane commented Mar 22, 2024

Just tried to build this PR and does not build

*** Building BayesMix executable ***
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/blocked_gibbs_algorithm.cc.o
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/conditional_algorithm.cc.o
[  4%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/base_algorithm.cc.o
[  5%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/marginal_algorithm.cc.o
[  7%] Building CXX object CMakeFiles/bayesmix.dir/src/algorithms/neal2_algorithm.cc.o
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/blocked_gibbs_algorithm.cc:1:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/blocked_gibbs_algorithm.h:5:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/conditional_algorithm.h:7:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/algorithms/base_algorithm.h:14:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/base_hierarchy.h:11:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/abstract_hierarchy.h:15:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/priors/abstract_prior_model.h:11:
In file included from /Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/likelihoods/states/includes.h:6:
/Users/jdtuck/Documents/Research/bayesmix/src/hierarchies/likelihoods/states/shape_rate_state.h:68:20: error: no member named 'sr_state' in 'bayesmix::AlgorithmState_ClusterState'
    shape = state_.sr_state().shape();
            ~~~~~~ ^

@jdtuck Thank you for checking this out. Some protobuf messages were missing. Now I'm able to build the betaGG hierarchy. 👍

@jdtuck
Copy link

jdtuck commented Mar 26, 2024

Awesome, I will check it out later this week.

@jdtuck
Copy link

jdtuck commented Mar 27, 2024

built on my end, will try to test it soon.

@jdtuck
Copy link

jdtuck commented Mar 28, 2024

@mberaha do you have a test example that I can compare with. I am having trouble

terminate called after throwing an instance of 'std::domain_error'
  what():  categorical_rng: Probabilities parameter is not a valid simplex. sum(Probabilities parameter) = -nan, but should be 1
Aborted (core dumped)
Error in run_mcmc("BetaGG", "DP", data1, g0_params_allprior, dp_params[2],  : 
  Something went wrong: 'run_mcmc()' exit with status 134

@mberaha
Copy link
Contributor Author

mberaha commented Apr 7, 2024

sorry for the latency, will try to work on it asap

@mberaha
Copy link
Contributor Author

mberaha commented Apr 19, 2024

Hey @jdtuck sorry for the delay.

I ran some tests and it looks fine to me. See the snippet below for a minimal working example through the Python interface

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta

from bayesmixpy import run_mcmc, build_bayesmix

build_bayesmix(4)

data = np.concatenate([beta(3, 10).rvs(100), beta(10, 3).rvs(100)])

dp = """
fixed_value {
    totalmass: 1.0
}
"""

g0 = """
fixed_values {
  a_rate: 2.0
  a_shape: 2.0
  b_rate: 2.0
  b_shape: 2.0
}
"""

algo = """
algo_id: "Neal8"
rng_seed: 20201124
iterations: 10000
burnin: 5000
init_num_clusters: 3
neal8_n_aux: 3
"""

dens_grid = np.linspace(0, 1, 500)

eval_dens, _, _, best_clus, _= run_mcmc(
        "BetaGG", "DP", data, g0, dp, algo, dens_grid,
        return_clusters=False, return_num_clusters=False,
        return_best_clus=True)

plt.plot(dens_grid, np.exp(np.mean(eval_dens, axis=0)), label="estimated_dens")
plt.plot(dens_grid, 0.5 * beta(3, 10).pdf(dens_grid) + 0.5 * beta(10, 3).pdf(dens_grid), label="true dens", color="red")
uniq_clus = np.unique(best_clus)

for u in uniq_clus:
    wh = np.where(best_clus == u)[0]
    plt.scatter(data[wh], np.zeros_like(wh))

plt.legend()

@mberaha mberaha marked this pull request as ready for review April 19, 2024 15:29
@mberaha
Copy link
Contributor Author

mberaha commented Apr 19, 2024

I pushed some docs, @TeoGiane I think this is ready for review.

@jdtuck any input on the code is most welcome as well!

@mberaha mberaha requested a review from TeoGiane April 19, 2024 15:29
@TeoGiane
Copy link
Contributor

I pushed some docs, @TeoGiane I think this is ready for review.

Great @mberaha, thanks a lot! I'll dig in a couple of days!

@jdtuck any input on the code is most welcome as well!

Totally agree, any ideas or suggesion is welcomed!

Copy link
Contributor

@TeoGiane TeoGiane left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems all good with proper DOCS. @mberaha, please merge at will.

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

Successfully merging this pull request may close these issues.

Support for Beta Mixture Model
3 participants