Dirichlet Process

[gcgibson]

Hi All, I really like Julia so far and Mamba is a great addition. I have been trying my hand at a dirichlet process implementation using the stick breaking construction.

#!/Applications/Julia-0.5.app/Contents/Resources/julia/bin/julia
using Mamba


## Data
line = Dict{Symbol, Any}(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]


function stick_breaking(beta)
	cum_prod = 1 - beta[1]
	weights = ones(5)
	for i = 2:5
		weights[i] = beta[i]*cum_prod
		cum_prod = cum_prod*(1-beta[i])
	end
	weights
end

model = Model(

  y = Stochastic(1,
    (mu) ->  Normal(mean(mu), 1),
    false
  ),

  mu = Logical(5,
  	(beta) -> stick_breaking(beta)
    ,
    true
  ),

  beta = Stochastic(5,
    () -> Beta(1, 2),
    true
  )

)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => randn(5)*1
  )
]


## Sampling Scheme
scheme = [NUTS(:beta)]

setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, line, inits, 10000, burnin=2500, thin=2, chains=1)
describe(sim)


## Posterior Predictive Distribution
ppd = predict(sim, :x)
describe(ppd)

[gcgibson]

Not sure why the code block formatting isn't working, however I get the following error

ERROR: LoadError: MethodError: Cannot convert an object of type Array{Float64,1} to an object of type Array{Float64,5}
This may have arisen from a call to the constructor Array{Float64,5}(...),
since type constructors fall back to convert methods.

So I assume there is something up with the dimensions of the initial values, but I don't see what is wrong.

[bdeonovic]

The error is because of the way you defined stochastic nodes.

Stochastic(n, ...) defines an n-dimensional node. Your mu and beta are 1-dimensional arrays not 5-dimensional arrays (i.e. matrices are 2-dimensional, etc).

---

You can do code blocks by using triple back-tick (i.e. ``` to open block and ``` to close). If you put ```julia as open block you will get julia syntax highlighting.

single backtick will give you inline code code . So does double backtick, but I'm not sure if that is any different.

[gcgibson]

Thanks, that is awesome! Here is my updated code for DP, I feel like I am very close but I'm still missing something, it seems is doesn't like me multiplying the weight by the Normal
ERROR: LoadError: MethodError: no method matching *(::Float64, ::Distributions.Normal{Mamba.ScalarStochastic})

so how can I take a vector of weights and produce a Gaussian mixture likelihood?

#!/Applications/Julia-0.5.app/Contents/Resources/julia/bin/julia
using Mamba


## Data
data = Dict{Symbol, Any}(
  :y => [1, 3, 3, 3, 5]
)

number_of_components = 5


function stick_breaking(beta)
  cum_prod = 1 - beta[1]
  weights = ones(5)
  weights[1] = beta[1]
  for i = 2:number_of_components
    weights[i] = beta[i]*cum_prod
    cum_prod = cum_prod*(1-beta[i])
  end
  weights
end


model = Model(

  y = Stochastic(1,
    (weights,mu,tau) -> UnivariateDistribution[weights[i]*Normal(mu, tau) for i in 1:number_of_components]

  ),
  weights = Logical(1,
    (beta) -> stick_breaking(beta),
    false
  ),

  beta = Stochastic(1,
    () -> Beta(1, 2)
    ),

  mu = Stochastic(
    () -> Normal(0,5)
  ),

  tau = Stochastic(
    () -> Gamma(1,1)
  )

)

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => data[:y],
    :beta =>rand(number_of_components)*1,
    :mix_dists => rand(number_of_components),
    :mu => 1,
    :tau => 1
  )
]


## Sampling Scheme
scheme = [NUTS(:beta)]

setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, data, inits, 10000, burnin=2500, thin=2, chains=1)
describe(sim)

[bdeonovic]

A Stochastic node has to be a take a function that returns a Distribution. In your function you are trying to multiply a distribution by a weight which is not defined i.e.

julia> 5 * Normal()
ERROR: MethodError: no method matching *(::Int64, ::Distributions.Normal{Float64})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:424
  *(::Real, ::Complex{Bool}) at complex.jl:254
  *(::Real, ::Complex) at complex.jl:266
  ...

What you want is a mixture distribution which is defined in Distributions.jl ie

julia> mu = collect(1:5);
       sigma = collect(1:5);
       w = rand(Dirichlet(ones(5)));
       MixtureModel( [Normal(mu[i], sigma[i]) for i in 1:5], w)
MixtureModel{Distributions.Normal{Float64}}(K = 5)
components[1] (prior = 0.1456): Distributions.Normal{Float64}(μ=1.0, σ=1.0)
components[2] (prior = 0.2764): Distributions.Normal{Float64}(μ=2.0, σ=2.0)
components[3] (prior = 0.4236): Distributions.Normal{Float64}(μ=3.0, σ=3.0)
components[4] (prior = 0.0463): Distributions.Normal{Float64}(μ=4.0, σ=4.0)
components[5] (prior = 0.1081): Distributions.Normal{Float64}(μ=5.0, σ=5.0)

[gcgibson]

Thanks!

[bdeonovic]

Also just want to point out that until recently you couldn't use MixtureModel in Mamba because there were some missing functions in Distributions.jl which was fixed recently, so I am not sure if you need to checkout Distributions.jl master or if they have already put out another release. (see #122)

[gcgibson]

Thanks! when I finish this up should I submit a PR to the examples? Would other people be interested in this?

[bdeonovic]

Feel free to contribute an example, try to match the format of the others.