In [1]:
using Turing, Distributions
using PyPlot, PyCall
using Mamba: describe

WARNING: Environment variable CMDSTAN_HOME not found. Use set_cmdstan_home!.

[Turing]: AD chunk size is set as 40

In [2]:
y = [ 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 1.0, 1.0 ];
N = length(y);  K = 3;

In [3]:
plot(y)

Out[3]:
1-element Array{PyCall.PyObject,1}:
PyObject <matplotlib.lines.Line2D object at 0x7f9d37776490>
In [5]:
@model BayesHmm(y) = begin
s = tzeros(Int, N)
m = Vector{Real}(K)
T = Vector{Vector{Real}}(K)
for i = 1:K
T[i] ~ Dirichlet(ones(K)/K)
# m[i] ~ Normal(1, 0.1) # Defining m this way causes label-switching problem.
m[i] ~ Normal(i, 0.01)
end
s[1] ~ Categorical(ones(Float64, K)/K)
for i = 2:N
s[i] ~ Categorical(vec(T[s[i-1]]))
y[i] ~ Normal(m[s[i]], 0.01)
end
return(s, m)
end

Out[5]:
BayesHmm (generic function with 2 methods)
In [11]:
g = Gibbs(500, HMC(1, 0.2, 3, :m, :T), PG(50, 1, :s))
c = sample(BayesHmm(y), g);

[Turing]:  Assume - T is a parameter
@~(::ANY, ::ANY) at compiler.jl:76
[Turing]:  Assume - m is a parameter (ignoring m found in global scope)
@~(::ANY, ::ANY) at compiler.jl:76
[Turing]:  Assume - s is a parameter (ignoring s found in global scope)
@~(::ANY, ::ANY) at compiler.jl:76
[Turing]:  Observe - y is an observation
@~(::ANY, ::ANY) at compiler.jl:57

[Gibbs] Sampling... 99%  ETA: 0:00:01
[Gibbs] Finished with
Running time    = 53.564168918000036;

[Gibbs] Sampling...100% Time: 0:00:54

In [12]:
describe(c)

Iterations = 1:500
Thinning interval = 1
Chains = 1
Samples per chain = 500

Empirical Posterior Estimates:
Mean                     SD                           Naive SE                        MCSE                 ESS
T[2][1]  4.792786378×10⁻¹    0.0000000000000015003021363   0.000000000000000067095551270   0.000000000000000000000000 500.000000
T[2][2]  2.784222581×10⁻¹    0.0000000000000001667002374   0.000000000000000007455061252   0.000000000000000000000000 500.000000
T[2][3]  2.422991041×10⁻¹    0.0000000000000003611838476   0.000000000000000016152632713   0.000000000000000013877788 500.000000
lf_num     2.9940000×10⁰    0.1341640786499871118575555   0.005999999999999987981835758   0.005999999999999961093622 500.000000
s[2]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[4]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[9]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[1]     2.0400000×10⁰    0.5127629585914010856839695   0.022931456635085653572581066   0.121202310209005501007162  17.898288
s[6]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[14]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
T[1][1]  8.065717592×10⁻¹    0.0000000000000011113349158   0.000000000000000049700408348   0.000000000000000000000000 500.000000
T[1][2]  1.795871242×10⁻¹    0.0000000000000001389168645   0.000000000000000006212551044   0.000000000000000000000000 500.000000
T[1][3] 1.3841116588×10⁻²    0.0000000000000000017364608   0.000000000000000000077656888   0.000000000000000000000000 500.000000
elapsed  1.071283378×10⁻¹    0.0597890418667208489722498   0.002673847238471375054730261   0.003990129674620669468499 224.527307
s[12]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[8]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[11]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
T[3][1]  6.169442702×10⁻¹    0.0000000000000004445339663   0.000000000000000019880163339   0.000000000000000000000000 500.000000
T[3][2]  3.444090352×10⁻¹    0.0000000000000004445339663   0.000000000000000019880163339   0.000000000000000000000000 500.000000
T[3][3] 3.8646694554×10⁻²    0.0000000000000000138916864   0.000000000000000000621255104   0.000000000000000000000000 500.000000
epsilon    2.0000000×10⁻¹    0.0000000000000004723173392   0.000000000000000021122673548   0.000000000000000013877788 500.000000
s[7]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[10]     1.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
m[1]    2.71155628×10⁰    0.0000000000000008890679326   0.000000000000000039760326679   0.000000000000000000000000 500.000000
s[3]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
m[2]    1.89621886×10⁰    0.0000000000000015558688821   0.000000000000000069580571688   0.000000000000000111022302 196.392786
s[5]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[15]     3.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
s[13]     2.0000000×10⁰    0.0000000000000000000000000   0.000000000000000000000000000   0.000000000000000000000000        NaN
lp   -4.46911188×10⁴ 3656.1798791090955091931391507 163.509334953098033338392269798 163.360946325820577840204351 500.000000
m[3]  7.111201435×10⁻¹    0.0000000000000010002014242   0.000000000000000044730367514   0.000000000000000000000000 500.000000

Quantiles:
2.5%             25.0%             50.0%             75.0%             97.5%
T[2][1]  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹  4.792786378×10⁻¹
T[2][2]  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹  2.784222581×10⁻¹
T[2][3]  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹  2.422991041×10⁻¹
lf_num     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
s[2]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
s[4]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
s[9]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰
s[1]     1.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     3.0000000×10⁰
s[6]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
s[14]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
T[1][1]  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹  8.065717592×10⁻¹
T[1][2]  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹  1.795871242×10⁻¹
T[1][3] 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻² 1.3841116588×10⁻²
elapsed 8.2674258675×10⁻²   9.60775315×10⁻²   1.06401758×10⁻¹  1.098029557×10⁻¹  1.302064362×10⁻¹
s[12]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
s[8]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰
s[11]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
T[3][1]  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹  6.169442702×10⁻¹
T[3][2]  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹  3.444090352×10⁻¹
T[3][3] 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻² 3.8646694554×10⁻²
epsilon    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹    2.0000000×10⁻¹
s[7]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
s[10]     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰     1.0000000×10⁰
m[1]    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰    2.71155628×10⁰
s[3]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
m[2]    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰    1.89621886×10⁰
s[5]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
s[15]     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰     3.0000000×10⁰
s[13]     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰     2.0000000×10⁰
lp   -4.45298711×10⁴   -4.45288443×10⁴   -4.45270086×10⁴   -4.45270086×10⁴   -4.45270086×10⁴
m[3]  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹  7.111201435×10⁻¹


In [14]:
m = c[:m][111];
s = c[:s][111];
PyPlot.plot(y, linestyle="None", marker="+", color = "r")
PyPlot.plot(m[s], linestyle="-", marker="*", color = "b")

Out[14]:
1-element Array{PyCall.PyObject,1}:
PyObject <matplotlib.lines.Line2D object at 0x7f9d33551050>