brian-j-smith / Mamba.jl

Markov chain Monte Carlo (MCMC) for Bayesian analysis in julia
Other
254 stars 52 forks source link

Mixtures. The Jain-Neal split-merge samplers #106

Open AngelBerihuete opened 8 years ago

AngelBerihuete commented 8 years ago

Dear Mamba users,

I wonder if it is possible to include the Jain-Neal split-merge samplers in order to do a nonparametric Bayesian mixture and/or Bayesian clustering as @jwmi does with BayesianMixtures.jl for mixture of finite mixtures. Maybe the best way to do that it would be to incorporate the BayesianMixtures.jl as sampling function ... but I don't know if that thing is possible.

Here a reproducible example in order to do inference, a three components mixture of 2D Gaussians.

using Mamba
srand(12345)

spread = 5

mu_0 = zeros(2)
mu_1 = spread*ones(2)
mu_2 = -spread*ones(2)
sigma_0 = [1.0 0.0; 0.0 1.0]
sigma_1 = [2.0 0.0; 0.0 1.0]
sigma_2 = [2.0 0.0; 0.0 1.0]

X = MixtureModel(MvNormal, [(mu_0, sigma_0), (mu_1, sigma_1), (mu_2, sigma_2),
 ], [0.2, 0.2, 0.6])

data = Dict{Symbol, Any}(
  :y => rand(X, 100)'
)

data[:N] = size(data[:y])[1]
data[:K] = 3
data[:alpha] = ones(data[:K])
data[:mean] = ones(2)
data[:var] = eye(2)
data[:P0] = [200.0 0.0; 0.0 1.0]

model = Model(

  P = Stochastic(1,
    alpha -> Dirichlet(alpha)
  ),

  T = Stochastic(1,
    (N, P) -> UnivariateDistribution[Categorical(P) for i in 1:N], 
    false
  ),

  mu = Stochastic(2,
    (K, mu0, sigma0) -> MultivariateDistribution[MvNormal(mu0, sigma0) 
    for k in 1:K]
  ),

  mu0 = Stochastic(1,
    (mean, var) -> MvNormal(mean, var),
    false
    ),

  sigma0 = Stochastic(2,
    (P0) -> InverseWishart(2.0, P0),
    false
  ),

  y = Stochastic(2,
    (mu, N, T) -> 
    begin
      MultivariateDistribution[
      begin
        mu_aux = mu[Int(T[i]),:]
        MvNormal(mu_aux, eye(2))
      end
      for i in 1:N
      ]
    end,
    false
  )
)

inits = [
  Dict(:y => data[:y], :T => fill(1, data[:N]), :P => ones(data[:K])/data[:K],
    :mu0 => zeros(2), :sigma0 => data[:P0], :mu => repmat([0 0], data[:K],1)),
  Dict(:y => data[:y], :T => fill(1, data[:N]), :P => ones(data[:K])/data[:K],
    :mu0 => zeros(2), :sigma0 => data[:P0], :mu => repmat([0 0], data[:K], 1))
]

# Sampling Scheme
# ---------------
scheme = [DGS(:T),
          Slice(:mu, 1.0),
          SliceSimplex(:P, scale=0.75)]

setsamplers!(model, scheme)

# MCMC Simulations
# ----------------

sim = mcmc(model, data, inits, 1000, burnin=250, thin=2, chains=2)
bdeonovic commented 8 years ago

After skimming the paper it looks like something reasonable that could be implemented. I won't have time for awhile (December? January?) But feel free to try and take a stab at it yourself! Brian and I can answer any questions you have.

jwmi commented 8 years ago

Hi everyone, It would be great if you want to use the BayesianMixtures package in Mamba. I'd be interested to stay in the loop. Cheers, Jeff

On Mon, Oct 17, 2016 at 8:59 PM, Benjamin Deonovic <notifications@github.com

wrote:

After skimming the paper it looks like something reasonable that could be implemented. I won't have time for awhile (December? January?) But feel free to try and take a stab at it yourself! Brian and I can answer any questions you have.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/brian-j-smith/Mamba.jl/issues/106#issuecomment-254377591, or mute the thread https://github.com/notifications/unsubscribe-auth/AHrTJk7vMrvu6D8UiFe-_Z2uuPzYtqCIks5q1BoGgaJpZM4KYt17 .

AngelBerihuete commented 8 years ago

Dear @bdeonovic and @jwmi I am afraid that I would not be able to code such thing in a short period of time. It is my second week using Julia and I think this task has a lot of specific requirements. My proposal tries to resolve a problem I proposed to @jwmi last week in other context.

Anyway, only for curiosity, do you have any example (documentation) to code a sampling function? What would be the steps to incorporate BaesianMixtures.jl to Mamba? Cheers, Angel

bdeonovic commented 8 years ago

A good example of a basic sampler can be found in the tutorial where a Gibbs sampler is coded up for Bayesian linear regression