Add DIC

Project.toml

@@ -23,6 +23,7 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsModelComparisons = "854dedd9-9477-4a25-907d-7fd989bfdd01"
 TableTraits = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
@@ -40,6 +41,7 @@ PrettyTables = "0.9, 0.10, 0.11"
 RecipesBase = "0.7, 0.8, 1.0"
 SpecialFunctions = "^0.8, 0.9, 0.10, 1.0"
 StatsBase = "0.32, 0.33"
+StatsModelComparisons = "0.1.1"
 TableTraits = "0.4, 1"
 Tables = "1"
 julia = "1"

src/MCMCChains.jl

@@ -24,6 +24,8 @@ import Random
 import Serialization
 import Statistics: std, cor, mean, var, mean!
 
+import StatsModelComparisons
+
 export Chains, chains, chainscat
 export setrange, resetrange
 export set_section, get_params, sections, sort_sections, setinfo

src/modelstats.jl

@@ -3,47 +3,13 @@ export dic
 #################### Posterior Statistics ####################
 
 """
-    dic(chain::Chains, logpdf::Function) -> (DIC, pD)
+    dic(chain::Chains, loglik::Symbol)
 
-Compute the deviance information criterion (DIC).
-(Smaller is better)
+Compute the deviance information criterion (DIC) from `chain` on posterior log likelihood samples specified by parameter name `loglik`.
 
 Note: DIC assumes that the posterior distribution is approx. multivariate Gaussian and tends to select overfitted models.
-
-## Returns:
-* `DIC`: The calculated deviance information criterion
-* `pD`: The effective number of parameters
-
-## Usage:
-
-```
-chn ... # sampling results
-lpfun = function f(chain::Chains) # function to compute the logpdf values
-    niter, nparams, nchains = size(chain)
-    lp = zeros(niter + nchains) # resulting logpdf values
-    for i = 1:nparams
-        lp += map(p -> logpdf( ... , x), Array(chain[:,i,:]))
-    end
-    return lp
-end
-
-DIC, pD = dic(chn, lpfun)
-```
-
 """
-function dic(chain::Chains, logpdf::Function)
-
-    # expectation of each parameter
-    Eθ = reshape(mean(Array(chain), dims = [1,3]), 1,:,1)
-    Echain = Chains(Eθ)
-    EθD = -2*mean(logpdf(Echain))
-
-    D = -2*logpdf(chain)
-    ED = mean(D)
-
-    pD = 2*(ED - EθD)
-
-    DIC = EθD + pD
-
-    return DIC, pD
+function StatsModelComparisons.dic(chain::Chains, loglik::Symbol)
+    lps = Array(chain[:, loglik, :])
+    return dic(vec(lps))
 end