Merge branch 'breaking' into mhauru/dppl-0.35

Author: penelopeysm

.github/workflows/Tests.yml

@@ -14,7 +14,7 @@ concurrency:
 jobs:
   test:
     # Use matrix.test.name here to avoid it taking up the entire window width
-    name: test ${{matrix.test.name}} (${{ matrix.runner.os }}, ${{ matrix.runner.version }}, ${{ matrix.runner.arch }}, ${{ matrix.runner.num_threads }})
+    name: test ${{matrix.test.name}} (${{ matrix.runner.version }}, ${{ matrix.runner.os }}, ${{ matrix.runner.num_threads }})
     runs-on: ${{ matrix.runner.os }}
     continue-on-error: ${{ matrix.runner.version == 'pre' }}
 
@@ -40,49 +40,36 @@ jobs:
           # Default
           - version: '1'
             os: ubuntu-latest
-            arch: x64
-            num_threads: 1
-          # x86
-          - version: '1'
-            os: ubuntu-latest
-            arch: x86
             num_threads: 1
           # Multithreaded
           - version: '1'
             os: ubuntu-latest
-            arch: x64
             num_threads: 2
           # Windows
           - version: '1'
             os: windows-latest
-            arch: x64
             num_threads: 1
           # macOS
           - version: '1'
             os: macos-latest
-            arch: aarch64
             num_threads: 1
           # Minimum supported Julia version
           - version: 'min'
             os: ubuntu-latest
-            arch: x64
             num_threads: 1
           # Minimum supported Julia version, multithreaded
           - version: 'min'
             os: ubuntu-latest
-            arch: x64
             num_threads: 2
           # Pre-release Julia version
           - version: 'pre'
             os: ubuntu-latest
-            arch: x64
             num_threads: 1
 
     steps:
       - name: Print matrix variables
         run: |
           echo "OS: ${{ matrix.runner.os }}"
-          echo "Architecture: ${{ matrix.runner.arch }}"
           echo "Julia version: ${{ matrix.runner.version }}"
           echo "Number of threads: ${{ matrix.runner.num_threads }}"
           echo "Test arguments: ${{ matrix.test.args }}"
@@ -93,7 +80,6 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '${{ matrix.runner.version }}'
-          arch: ${{ matrix.runner.arch }}
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       # TODO: Use julia-actions/julia-runtest when test_args are supported

docs/Project.toml

@@ -3,4 +3,5 @@ Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
+DynamicPPL = "366bfd00-2699-11ea-058f-f148b4cae6d8"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"

docs/make.jl

@@ -4,6 +4,7 @@ using Turing
 # from those packages.
 using Distributions
 using Bijectors
+using DynamicPPL
 
 using DocumenterInterLinks
 

docs/src/api.md

@@ -127,7 +127,7 @@ Distributions.BernoulliLogit
 ### Predictions
 
 ```@docs
-predict
+DynamicPPL.predict
 ```
 
 ### Querying model probabilities and quantities

src/optimisation/Optimisation.jl

@@ -15,6 +15,7 @@ using Printf: Printf
 using ForwardDiff: ForwardDiff
 using StatsAPI: StatsAPI
 using Statistics: Statistics
+using LinearAlgebra: LinearAlgebra
 
 export maximum_a_posteriori, maximum_likelihood
 # The MAP and MLE exports are only needed for the Optim.jl interface.
@@ -235,11 +236,61 @@ end
 
 # Various StatsBase methods for ModeResult
 
-function StatsBase.coeftable(m::ModeResult; level::Real=0.95)
+"""
+    StatsBase.coeftable(m::ModeResult; level::Real=0.95, numerrors_warnonly::Bool=true)
+
+
+Return a table with coefficients and related statistics of the model. level determines the
+level for confidence intervals (by default, 95%).
+
+In case the `numerrors_warnonly` argument is true (the default) numerical errors encountered
+during the computation of the standard errors will be caught and reported in an extra
+"Error notes" column.
+"""
+function StatsBase.coeftable(m::ModeResult; level::Real=0.95, numerrors_warnonly::Bool=true)
     # Get columns for coeftable.
     terms = string.(StatsBase.coefnames(m))
     estimates = m.values.array[:, 1]
-    stderrors = StatsBase.stderror(m)
+    # If numerrors_warnonly is true, and if either the information matrix is singular or has
+    # negative entries on its diagonal, then `notes` will be a list of strings for each
+    # value in `m.values`, explaining why the standard error is NaN.
+    notes = nothing
+    local stderrors
+    if numerrors_warnonly
+        infmat = StatsBase.informationmatrix(m)
+        local vcov
+        try
+            vcov = inv(infmat)
+        catch e
+            if isa(e, LinearAlgebra.SingularException)
+                stderrors = fill(NaN, length(m.values))
+                notes = fill("Information matrix is singular", length(m.values))
+            else
+                rethrow(e)
+            end
+        else
+            vars = LinearAlgebra.diag(vcov)
+            stderrors = eltype(vars)[]
+            if any(x -> x < 0, vars)
+                notes = []
+            end
+            for var in vars
+                if var >= 0
+                    push!(stderrors, sqrt(var))
+                    if notes !== nothing
+                        push!(notes, "")
+                    end
+                else
+                    push!(stderrors, NaN)
+                    if notes !== nothing
+                        push!(notes, "Negative variance")
+                    end
+                end
+            end
+        end
+    else
+        stderrors = StatsBase.stderror(m)
+    end
     zscore = estimates ./ stderrors
     p = map(z -> StatsAPI.pvalue(Distributions.Normal(), z; tail=:both), zscore)
 
@@ -251,7 +302,7 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95)
     level_ = 100 * level
     level_percentage = isinteger(level_) ? Int(level_) : level_
 
-    cols = [estimates, stderrors, zscore, p, ci_low, ci_high]
+    cols = Vector[estimates, stderrors, zscore, p, ci_low, ci_high]
     colnms = [
         "Coef.",
         "Std. Error",
@@ -260,6 +311,10 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95)
         "Lower $(level_percentage)%",
         "Upper $(level_percentage)%",
     ]
+    if notes !== nothing
+        push!(cols, notes)
+        push!(colnms, "Error notes")
+    end
     return StatsBase.CoefTable(cols, colnms, terms)
 end
 

test/Project.toml

@@ -59,7 +59,7 @@ LinearAlgebra = "1"
 LogDensityProblems = "2"
 LogDensityProblemsAD = "1.4"
 MCMCChains = "5, 6"
-Mooncake = "0.4.61"
+Mooncake = "0.4.95"
 NamedArrays = "0.9.4, 0.10"
 Optim = "1"
 Optimization = "3, 4"

test/optimisation/Optimisation.jl

@@ -627,6 +627,65 @@ using Turing
         maximum_likelihood(m; adtype=adbackend)
         maximum_a_posteriori(m; adtype=adbackend)
     end
+
+    @testset "Collinear coeftable" begin
+        xs = [-1.0, 0.0, 1.0]
+        ys = [0.0, 0.0, 0.0]
+
+        @model function collinear(x, y)
+            a ~ Normal(0, 1)
+            b ~ Normal(0, 1)
+            return y ~ MvNormal(a .* x .+ b .* x, 1)
+        end
+
+        model = collinear(xs, ys)
+        mle_estimate = Turing.Optimisation.estimate_mode(model, MLE())
+        tab = coeftable(mle_estimate)
+        @assert isnan(tab.cols[2][1])
+        @assert tab.colnms[end] == "Error notes"
+        @assert occursin("singular", tab.cols[end][1])
+    end
+
+    @testset "Negative variance" begin
+        # A model for which the likelihood has a saddle point at x=0, y=0.
+        # Creating an optimisation result for this model at the x=0, y=0 results in negative
+        # variance for one of the variables, because the variance is calculated as the
+        # diagonal of the inverse of the Hessian.
+        @model function saddle_model()
+            x ~ Normal(0, 1)
+            y ~ Normal(x, 1)
+            Turing.@addlogprob! x^2 - y^2
+            return nothing
+        end
+        m = saddle_model()
+        ctx = Turing.Optimisation.OptimizationContext(DynamicPPL.LikelihoodContext())
+        optim_ld = Turing.Optimisation.OptimLogDensity(m, ctx)
+        vals = Turing.Optimisation.NamedArrays.NamedArray([0.0, 0.0])
+        m = Turing.Optimisation.ModeResult(vals, nothing, 0.0, optim_ld)
+        ct = coeftable(m)
+        @assert isnan(ct.cols[2][1])
+        @assert ct.colnms[end] == "Error notes"
+        @assert occursin("Negative variance", ct.cols[end][1])
+    end
+
+    @testset "Same coeftable with/without numerrors_warnonly" begin
+        xs = [0.0, 1.0, 2.0]
+
+        @model function extranormal(x)
+            mean ~ Normal(0, 1)
+            return x ~ Normal(mean, 1)
+        end
+
+        model = extranormal(xs)
+        mle_estimate = Turing.Optimisation.estimate_mode(model, MLE())
+        warnonly_coeftable = coeftable(mle_estimate; numerrors_warnonly=true)
+        no_warnonly_coeftable = coeftable(mle_estimate; numerrors_warnonly=false)
+        @assert warnonly_coeftable.cols == no_warnonly_coeftable.cols
+        @assert warnonly_coeftable.colnms == no_warnonly_coeftable.colnms
+        @assert warnonly_coeftable.rownms == no_warnonly_coeftable.rownms
+        @assert warnonly_coeftable.pvalcol == no_warnonly_coeftable.pvalcol
+        @assert warnonly_coeftable.teststatcol == no_warnonly_coeftable.teststatcol
+    end
 end
 
 end