Towards v0.2.0

[xukai92]

This PR will address the rest issues in Plan for v0.2 #64

• More robust and modular way of detecting divergence More robust and modular way of detecting divergence  #16
• Return more information for each step Return more information for each step #59
• Implement Base.show for all user side types Implement Base.show for all user side types #36
• RFC Adapter types RFC Adapter types #28
• Instantiation of new AbstractMetrics Instantiation of new AbstractMetrics #45
• Simplifying Metric types by using AbstractPDMat Simplifying Metric types by using AbstractPDMat #34
• Adding noise to step size? Adding noise to step size? #10

[xukai92]

@yebai Can you remind me why the change of metric length trick is necessay? It it for potential dynamic space, or it's also necessary for static HMC?

[xukai92]

Was a bug here: we forgot to add the negative sign to neg_energy(z). @yebai

[xukai92]

All issues are resolved except from #34 - please check my comment in #34 (comment).

For #54, @cpfiffer agrees with that we can stay with returning statistic within single transition using NamedTuple. This would not introduce any breaking changes to Turing.jl for now.

[xukai92]

Looks like this

samples = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; progress=true);
Sampling100%|███████████████████████████████| Time: 0:00:04
  iteration:           100000
  acceptance_rate:     0.6724052307089154
  n_steps:             4
  hamiltonian_energy:  13.170723025896388
  numerical_error:     false
  is_accept:           true
  precondition:        DiagEuclideanMetric([0.479825, 0.53532, 0.75214 ...])
  step_size:           1.1207818211894387
  log_density:         -11.689584321232719
  tree_depth:          2
┌ Info: Finished 100000 sampling steps in 4.238515317 (s)
│   h = Hamiltonian(metric=DiagEuclideanMetric([0.479825, 0.53532, 0.75214 ...]))
│   τ = NUTS{Multinomial}(integrator=Leapfrog(ϵ=1.12), max_depth=10), Δ_max=1000.0)
│   EBFMI(Hs) = 112748.64011387681
└   mean(αs) = 0.8000706158165503

or this

samples = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; verbose=true);
┌ Info: Finished 2000 adapation steps
│   adaptor = StanHMCAdaptor(n_adapts=2000, pc=DiagPreconditioner, ssa=NesterovDualAveraging(γ=0.05, t_0=10.0, κ=0.75, δ=0.8, state.ϵ=1.1116666464941625), init_buffer=75, term_buffer=50)
│   τ.integrator = Leapfrog(ϵ=1.05)
└   h.metric = DiagEuclideanMetric([0.479825, 0.53532, 0.75214 ...])
┌ Info: Finished 100000 sampling steps in 3.190849454 (s)
│   h = Hamiltonian(metric=DiagEuclideanMetric([0.479825, 0.53532, 0.75214 ...]))
│   τ = NUTS{Multinomial}(integrator=Leapfrog(ϵ=1.11), max_depth=10), Δ_max=1000.0)
│   EBFMI(Hs) = 114838.007040791
└   mean(αs) = 0.8035625060120144

in AHMC ATM.

Looks like this in Turing.jl with MCMCChains.jl ATM.

julia> describe(chn, sections=:internals)
2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters         │ mean     │ std         │ naive_se    │ mcse        │ ess     │ r_hat    │
│     │ Symbol             │ Float64  │ Float64     │ Float64     │ Float64     │ Any     │ Any      │
├─────┼────────────────────┼──────────┼─────────────┼─────────────┼─────────────┼─────────┼──────────┤
│ 1   │ acceptance_rate    │ 0.996775 │ 0.00717163  │ 0.000226787 │ 0.000216056 │ 1000.0  │ 1.00124  │
│ 2   │ eval_num           │ 7002.0   │ 0.0         │ 0.0         │ 0.0         │ NaN     │ NaN      │
│ 3   │ hamiltonian_energy │ 6.7005   │ 1.45964     │ 0.0461578   │ 0.102329    │ 309.238 │ 1.00569  │
│ 4   │ is_accept          │ 0.998    │ 0.044699    │ 0.00141351  │ 0.00133333  │ 1000.0  │ 0.998999 │
│ 5   │ log_density        │ -5.67199 │ 1.00641     │ 0.0318255   │ 0.0895646   │ 133.668 │ 1.00688  │
│ 6   │ lp                 │ -4.82565 │ 2.66587e-15 │ 8.43022e-17 │ 0.0         │ 4.01606 │ 5.73882  │
│ 7   │ n_steps            │ 5.0      │ 0.0         │ 0.0         │ 0.0         │ NaN     │ NaN      │
│ 8   │ step_size          │ 0.1      │ 1.38847e-16 │ 4.39074e-18 │ 0.0         │ 4.01606 │ 5.0447   │

Quantiles

│ Row │ parameters         │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%    │
│     │ Symbol             │ Float64  │ Float64  │ Float64  │ Float64  │ Float64  │
├─────┼────────────────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1   │ acceptance_rate    │ 0.978634 │ 0.996369 │ 0.999995 │ 1.0      │ 1.0      │
│ 2   │ eval_num           │ 7002.0   │ 7002.0   │ 7002.0   │ 7002.0   │ 7002.0   │
│ 3   │ hamiltonian_energy │ 4.82941  │ 5.6462   │ 6.41803  │ 7.38837  │ 10.4049  │
│ 4   │ is_accept          │ 1.0      │ 1.0      │ 1.0      │ 1.0      │ 1.0      │
│ 5   │ log_density        │ -8.27467 │ -6.15494 │ -5.35321 │ -4.94637 │ -4.62537 │
│ 6   │ lp                 │ -4.82565 │ -4.82565 │ -4.82565 │ -4.82565 │ -4.82565 │
│ 7   │ n_steps            │ 5.0      │ 5.0      │ 5.0      │ 5.0      │ 5.0      │
│ 8   │ step_size          │ 0.1      │ 0.1      │ 0.1      │ 0.1      │ 0.1      │

[yebai]

@yebai Can you remind me why the change of metric length trick is necessay? It it for potential dynamic space, or it's also necessary for static HMC?

When AHMC metric is created, we might not know the model dimension yet. This trick removes the dependency of AHMC algorithm configuration on the model dimension in Turing.

[xukai92]

@yebai Can you remind me why the change of metric length trick is necessay? It it for potential dynamic space, or it's also necessary for static HMC?

When AHMC metric is created, we might not know the model dimension yet. This trick removes the dependency of AHMC algorithm configuration on the model dimension in Turing.

I see. Thanks!

[yebai]

Thanks, @xukai92. This PR looks correct to me. I made some minor style-related comments below.

[xukai92]

TODOs

• wait for response from Towards v0.2.0 #81 (comment)
• rename Divergence to Termination
• reimplement or replace renew by reconstruct
• update named tuple signature

[yebai]

Hi @xukai92, I'm ready to merge this PR once you're happy.

[xukai92]

Great. I will merge it in 30 mins.

[xukai92]

@yebai I didn't compute the new sqrtM⁻¹ based on the new M⁻¹.