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>