Merge 30e7b96 into bed6655

.travis.yml

@@ -11,6 +11,7 @@ notifications:
   email: false
 matrix:
   allow_failures:
+    - julia: 1.0                # FIXME remove when MCMCDiagnostics catches up
     - julia: nightly
 # uncomment the following lines for debugging Documenter.jl
 # env:

Project.toml

@@ -0,0 +1,27 @@
+name = "DynamicHMC"
+uuid = "bbc10e6e-7c05-544b-b16e-64fede858acb"
+authors = ["Tamas K. Papp <tkpapp@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+
+[extras]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+MCMCDiagnostics = "6e857e4b-079a-58c4-aeab-bc2670384359"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Distributions", "Documenter", "ForwardDiff", "MCMCDiagnostics", "StatsBase", "Suppressor", "Test"]

src/DynamicHMC.jl

@@ -1,4 +1,3 @@
-__precompile__()
 module DynamicHMC
 
 import Base: rand, length, show

src/hamiltonian.jl

@@ -159,7 +159,7 @@ phasepoint_in(H::Hamiltonian, q, p) = PhasePoint(q, p, H.ℓ(q))
 Extend a position `q` to a phasepoint with a random momentum according to the
 kinetic energy of `H`.
 """
-rand_phasepoint(rng, H, q) = phasepoint_in(H, q, rand(rng, H.κ))
+rand_phasepoint(rng::AbstractRNG, H, q) = phasepoint_in(H, q, rand(rng, H.κ))
 
 """
     $SIGNATURES

test/runtests.jl

@@ -1,6 +1,75 @@
-include("setup-and-utilities.jl")
-include("test-utilities.jl")
-include("test-basics.jl")
+using DynamicHMC
+
+using DynamicHMC:
+    # Hamiltonian
+    GaussianKE, Hamiltonian, PhasePoint, neg_energy, phasepoint_in,
+    rand_phasepoint, leapfrog, move, is_valid_ℓq,
+    # building blocks
+    rand_bool, TurnStatistic, combine_turnstats, Proposal,
+    combined_logprob_logweight, combine_proposals,
+    DivergenceStatistic, combine_divstats, divergence_statistic, isdivergent,
+    get_acceptance_rate, isturning, adjacent_tree, sample_trajectory,
+    # stepsize
+    InitialStepsizeSearch, find_initial_stepsize,
+    # transitions and tuning
+    NUTS_transition, NUTS_init, StepsizeTuner, StepsizeCovTuner, tune,
+    TunerSequence, bracketed_doubling_tuner
+
+using Test
+
+using ArgCheck: @argcheck
+using DataStructures
+import DiffResults
+using DiffResults: MutableDiffResult
+using Distributions
+using DocStringExtensions: SIGNATURES
+import ForwardDiff
+using LinearAlgebra
+using MCMCDiagnostics: effective_sample_size, potential_scale_reduction
+using Parameters
+import Random
+using Random: randn, rand
+using StatsBase: mean_and_cov, mean_and_std
+using StatsFuns: logaddexp
+using Statistics: mean, quantile
+using Suppressor
+
+
+# general test environment
+
+const RNG = Random.GLOBAL_RNG   # shorthand
+Random.seed!(RNG, UInt32[0x23ef614d, 0x8332e05c, 0x3c574111, 0x121aa2f4])
+
+"Tolerant testing in a CI environment."
+const RELAX = (k = "CONTINUOUS_INTEGRATION"; haskey(ENV, k) && ENV[k] == "true")
+
+"""
+    $SIGNATURES
+
+Random positive definite matrix of size `n` x `n` (for testing).
+"""
+function rand_Σ(::Type{Symmetric}, n)
+    A = randn(n, n)
+    Symmetric(A'*A .+ 0.01)
+end
+
+rand_Σ(::Type{Diagonal}, n) = Diagonal(randn(n).^2 .+ 0.01)
+
+rand_Σ(n::Int) = rand_Σ(Symmetric, n)
+
+## use MvNormal as a test distribution
+(ℓ::MvNormal)(p) = DiffResults.DiffResult(logpdf(ℓ, p), (gradlogpdf(ℓ, p), ))
+
+"Random Hamiltonian `H` with phasepoint `z`, with dimension `K`."
+function rand_Hz(K)
+    μ = randn(K)
+    Σ = rand_Σ(K)
+    κ = GaussianKE(inv(rand_Σ(Diagonal, K)))
+    H = Hamiltonian(MvNormal(μ, Matrix(Σ)), κ)
+    z = rand_phasepoint(RNG, H, μ)
+    H, z
+end
+
 include("test-Hamiltonian-leapfrog.jl")
 include("test-buildingblocks.jl")
 include("test-stepsize.jl")

test/test-Hamiltonian-leapfrog.jl

@@ -1,21 +1,71 @@
-import DynamicHMC:
-    GaussianKE, Hamiltonian, PhasePoint, neg_energy, phasepoint_in,
-    rand_phasepoint, leapfrog, move, is_valid_ℓq, InitialStepsizeSearch,
-    find_initial_stepsize
+
+# utility functions
 
-import DiffResults: MutableDiffResult
+"Test `loggradient` vs autodiff `neg_energy`."
+function test_loggradient(ℓ, x)
+    ∇ = DynamicHMC.loggradient(ℓ, x)
+    ∇_AD = ForwardDiff.gradient(x->neg_energy(ℓ,x), x)
+    @test ∇ ≈ ∇_AD
+end
+
+"""
+    $SIGNATURES
+
+Return a reasonable estimate for the largest stable stepsize (which may not be
+stable, but is a good starting point for finding that).
+
+`q` is assumed to be normal with variance `Σ`. `κ` is the kinetic energy.
+
+Using the transformation ``p̃ = W⁻¹ p``, the kinetic energy is
+
+``p'M⁻¹p = p'W⁻¹'W⁻¹p/2=p̃'p̃/2``
+
+Transforming to ``q̃=W'q``, the variance of which becomes ``W' Σ W``. Return the
+square root of its smallest eigenvalue, following Neal (2011, p 136).
+
+When ``Σ⁻¹=M=WW'``, this the variance of `q̃` is ``W' Σ W=W' W'⁻¹W⁻¹W=I``, and
+thus decorrelates the density perfectly.
+"""
+find_stable_ϵ(κ::GaussianKE, Σ) = √eigmin(κ.W'*Σ*κ.W)
+
+function find_stable_ϵ(H::Hamiltonian{Tℓ, Tκ}) where
+    {Tℓ <: Distribution{Multivariate,Continuous}, Tκ}
+    find_stable_ϵ(H.κ, cov(H.ℓ))
+end
+
+"""
+    $SIGNATURES
+
+Simulated mean and covariance of `N` values from `f()`.
+"""
+function simulated_meancov(f, N)
+    s = f()
+    K = length(s)
+    x = similar(s, (N, K))
+    for i in 1:N
+        x[i, :] = f()
+    end
+    m, C = mean_and_cov(x, 1)
+    vec(m), C
+end
+
+@testset "simulated meancov" begin
+    d = MvNormal([2,2], [2.0,1.0])
+    m, C = simulated_meancov(()->rand(d), 10000)
+    @test m ≈ mean(d) atol = 0.05 rtol = 0.1
+    @test C ≈ cov(d) atol = 0.05 rtol = 0.1
+end
 
-######################################################################
-# Hamiltonian and leapfrog
-######################################################################
+
+# testsets
 
 @testset "Gaussian KE full" begin
     for _ in 1:100
         K = rand(2:10)
         Σ = rand_Σ(Symmetric, K)
         κ = GaussianKE(inv(Σ))
         @test κ.Minv * κ.W * κ.W' ≈ Diagonal(ones(K))
-        m, C = simulated_meancov(()->rand(RNG, κ), 10000)
+        m, C = simulated_meancov(()->rand(κ), 10000)
         @test Matrix(Σ) ≈ C rtol = 0.1
         test_loggradient(κ, randn(K))
     end
@@ -28,7 +78,7 @@ end
         κ = GaussianKE(inv(Σ))
         # FIXME workaround for https://github.com/JuliaLang/julia/issues/28869
         @test κ.Minv * Matrix(κ.W) * Matrix(κ.W') ≈ Diagonal(ones(K))
-        m, C = simulated_meancov(()->rand(RNG, κ), 10000)
+        m, C = simulated_meancov(()->rand(κ), 10000)
         @test Matrix(Σ) ≈ C rtol = 0.1
         test_loggradient(κ, randn(K))
     end