brian-j-smith / Mamba.jl

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

Sample from mixture of VonMises distribution.(Ungiven Float64) #159

Open sosuts opened 4 years ago

sosuts commented 4 years ago

Hello, I’m a student majoring in biology in Japan. I usually use Python so I'm very new to Julia...

I am trying to fit my experimental data(~400,000 angles, -pi ~ +pi bimodal) with mixture of von mises distribution. I want to predict the parameters(mu, kappa, weights) and the number of components.

mydata
# Dict{Symbol,Any} with 1 entry:
  :angle => ....

Below is the model I want to describe.

kappa1 ~ Gamma(alpha=1., beta=1.0)
kappa2 ~ Gamma(alpha=1., beta=1.0)
mu1 ~ VonMises(mu=0, kappa=0.5)
mu2 ~ VonMises(mu=pi, kappa=0.5)
w ~ Dirichlet([0.5, 0.5])
y ~ w[1]*VonMIses(mu=mu1, kappa=kappa1) + w[2]*VonMises(mu=mu2, kappa=kappa2)

What I want is a posterior distribution of kappa, mu, and weights. From the previous post : Incompatible Initial Value?, it looks like I'm specifying a wrong dimension to those parameters. However, I can't figure out which code I must fix...

Any suggestions are welcomed. I'm sorry for very basic qustion. Thank you in advance.

The following is my code.

using Mamba
using Distributions

model = Model(
    kappa1 = Stochastic(1,
        () -> Gamma(1, 1), 
    ),
    kappa2 = Stochastic(1,
        () -> Gamma(1, 1),
    ), 
    mu1 = Stochastic(1,
        () -> VonMises(0, 0.5), 
    ),
    mu2 = Stochastic(1,
        () -> VonMises(π, 0.5),
    ), 
    w = Stochastic(1, 
        () -> Dirichlet(ones(2))
    ),
    y = Stochastic(1,
        (mu1, mu2, kappa1, kappa2) -> 
        MixtureModel([VonMises(mu1, kappa1), VonMises(mu2, kappa2)], w),
    false),
)
scheme = [NUTS([:y, :kappa1, :kappa2, :mu1, :mu2, :w])]
setsamplers!(model, scheme)

inits = [
  Dict{Symbol, Any}(
        :y => data[:angle],
        :kappa1 => rand(Gamma(1, 1)),
        :kappa2 => rand(Gamma(1, 1)),
        :mu1 => rand(VonMises(0, 0.5)), 
        :mu2 => rand(VonMises(π, 0.5)), 
        :w => [0.5, 0.5]
  )
]
## MCMC Simulations
nuts_inference = mcmc(model, data, inits, 5000, burnin=1000, thin=2, chains=1)

returns

ArgumentError: incompatible initial value for node : kappa2
sosuts commented 4 years ago

I've now encountered another error. Details are in MCMC of mixture model using Mamba (Julia discourse).

bdeonovic commented 4 years ago

Sorry for the late response. This appears to be a problem with Distributions package (the VonMises distribution) and parameter checking. I will look into it today and file a pull request if I can fix it.

sosuts commented 4 years ago

Thank you so much for your reply. I'm sorry for the late response. I'm glad if it's fixed.