trappmartin/adst.jl

## adst.jl
using Turing

# This is a crappy implementation of Isabel's model.
# For simplicity, I assume all data to be positiv-real valued.
# Isabel's paper: http://proceedings.mlr.press/v70/valera17a/valera17a.pdf

@model function adst(y, ::Type{TV}=Vector{Float64}) where {TV}

  N,D = size(y)

  π ~ Dirichlet(3, 0.1)
  z ~ filldist(Categorical(π), D)

  μ = TV(undef, D)
  β = TV(undef, D)
  λ = TV(undef, D)

  for d in 1:D
    if z[d] == 1
      μ[d] ~ Normal(0.0, 1.0)
      for i in 1:N
        y[i] ~ Normal(μ[d], 1.0)
      end
    elseif z[d] == 2
      β[d] ~ Gamma(1.0, 1.0)
      for i in 1:N
        y[i] ~ Gamma(1.0, β[d])
      end
    else
      λ[d] ~ Gamma(1.0, 1.0)
      for i in 1:N
        y[i] ~ Exponential(λ[d])
      end
    end
  end
end

x = abs.(rand(100, 3)  + randn(100,3) .* randn(3)')
model = adst(x)

@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=0)

@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=1)

@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=2)
	using Turing

	# This is a crappy implementation of Isabel's model.
	# For simplicity, I assume all data to be positiv-real valued.
	# Isabel's paper: http://proceedings.mlr.press/v70/valera17a/valera17a.pdf

	@model function adst(y, ::Type{TV}=Vector{Float64}) where {TV}

	N,D = size(y)

	π ~ Dirichlet(3, 0.1)
	z ~ filldist(Categorical(π), D)

	μ = TV(undef, D)
	β = TV(undef, D)
	λ = TV(undef, D)

	for d in 1:D
	if z[d] == 1
	μ[d] ~ Normal(0.0, 1.0)
	for i in 1:N
	y[i] ~ Normal(μ[d], 1.0)
	end
	elseif z[d] == 2
	β[d] ~ Gamma(1.0, 1.0)
	for i in 1:N
	y[i] ~ Gamma(1.0, β[d])
	end
	else
	λ[d] ~ Gamma(1.0, 1.0)
	for i in 1:N
	y[i] ~ Exponential(λ[d])
	end
	end
	end
	end

	x = abs.(rand(100, 3) + randn(100,3) .* randn(3)')
	model = adst(x)

	@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=0)

	@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=1)

	@time sample(model, Gibbs(PG(10,:z), HMC(0.01, 5, :π, :μ, :β, :λ)), 100, chain_type=Any, specialize_after=2)