# Compute Deviance Information Criterion for a Bayesian Model

**URL:** <https://discourse.julialang.org/t/compute-deviance-information-criterion-for-a-bayesian-model/132195>\
**Category:** Probabilistic Programming\
**Tags:** question, turing, mcmc\
**Created:** [September 8, 2025, 6:22pm UTC](https://discourse.julialang.org/t/compute-deviance-information-criterion-for-a-bayesian-model/132195 "2025-09-08T18:22:00Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![dhruva-ahuja](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva-ahuja/32/218137_2.png) [@dhruva-ahuja](https://discourse.julialang.org/u/dhruva-ahuja)\
**Post date:** [September 8, 2025, 6:22pm UTC](https://discourse.julialang.org/t/compute-deviance-information-criterion-for-a-bayesian-model/132195/1 "2025-09-08T18:22:00Z")

</div>

Hello,

I am struggling to compute Deviance Information Criterion (DIC) in Julia as Turing/MCMCDiagnosticTools do not have an implementation of DIC as in R `dic.samples(mod2, n.iter=1e3)`.

Backdrop: I am enrolled in a course on Bayesian Statistics on Coursera where these models are created in R. I am new to Julia and found that Turing.jl has a very expressive syntax (chef’s kiss), I decided to the coursework in Julia. In assignments and quizzers, many times I am required to compute DIC, so I thought of implementing it on my own but my results are not close to what I get from R’s `dic.samples(mod, n.iter=1e3` function.

Here is my implementation for a simpler case, any help to figure out what I am missing here would be extremely appreciated:

```julia
using RDatasets, DataFrames
using Turing
using CategoricalArrays
using Distributions
using Statistics

# Load data
df = RDatasets.dataset("datasets", "warpbreaks")
println(first(df, 5))

# Create categorical variable and log-transform
df.cat = levelcode.(categorical(string.(df.Tension, "-", df.Wool)))
df.logBreaks = log.(df.Breaks)

# ANOVA model
@model function anova_t(X, y)
	num_cat = length(unique(X.cat))
	num_obs = length(y)
	μ ~ filldist(Normal(0.0, 1e3), num_cat)
	σ² ~ filldist(InverseGamma(0.5, 0.5), num_cat)
	for i in 1:num_obs
		y[i] ~ Normal(μ[X.cat[i]], sqrt(σ²[X.cat[i]]))
	end
end

# Fit model
model_1 = anova_t(df[:, [:cat]], df.logBreaks)
chains = sample(model_1, NUTS(), 5000)

# Log-likelihood function
function log_likelihood(dat, params)
	μ = params[iter=1, var=[Symbol("μ[$i]") for i in 1:length(unique(dat.cat))]]
	σ = sqrt.(params[iter=1, var=[Symbol("σ²[$i]") for i in 1:length(unique(dat.cat))]])
	sum(logpdf.(Normal.(μ.data[dat.cat], σ[dat.cat]), dat.logBreaks))
end

# Deviance calculations
mean_deviance = -2 * mean([log_likelihood(df, chains.value[iter=[i], var=:, chain=1]) for i in 1:size(chains, 1)])
deviance_at_mean = -2 * log_likelihood(df, mean(chains.value[iter=:, var=:, chain=1], dims=1))

println("mean_deviance: $mean_deviance") # prints mean_deviance: 60.31756252213876
println("p_D: $(mean_deviance - deviance_at_mean)") # prints p_D: 4.72884678655096
println("DIC: $(2 * mean_deviance - deviance_at_mean)") # prints DIC: 65.04640930868972

describe(chains)

```

While my R code prints (Which is correct as per the course)

```julia-auto
Mean deviance: 60.26 
penalty 14.68 
Penalized deviance: 74.94

```

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [September 9, 2025, 9:05am UTC](https://discourse.julialang.org/t/compute-deviance-information-criterion-for-a-bayesian-model/132195/2 "2025-09-09T09:05:46Z")

</div>

I cannot help directly with implementation and therefore I do not want to send you down a rabbit hole, but…

I have implemented AIC and BIC in various languages, comparing them to output from other languages. I can tell you that they all treat the constants different (e.g. \sqrt{\pi}). Some include them, some don’t. This doesn’t change things like minima and maxima, but you cannot compare the numbers per se. If you can find the exact function that R uses in the manual you can easily rule that out.