> For the complete documentation index, see [llms.txt](https://mc.hohoweiya.xyz/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://mc.hohoweiya.xyz/mcintegration/rao-blackwell.md).

# Rao-Blackwellization

## Recall: Rao-Blackwell Theorem

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5MAmYLmHctczsH%2Frao-blackwell-thm.png?generation=1574778107726132\&alt=media)

(Lehmann & Casella, 1998)

**Rao-Blackwellization** is that

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5OdILLPK3W9jfA%2Frao-blackwellization.png?generation=1574778107792770\&alt=media)

(Robert & Casella, 2005)

## Example: Student's $$t$$ expectation

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5QAKeGRJC_U6-w%2Ft-dist.png?generation=1574778107621495\&alt=media)

(Robert & Casella, 2005)

Here Gamma is in (α, β) form.

Alternative way to look at such decomposition

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5S8KjRB3VNqxxc%2Ft-dist2.png?generation=1574778107630441\&alt=media)

(Robert & Casella, 2010)

Compare the performance with two cases,

```julia
using Distributions
using Plots

function rt_dickey(n, μ, ν, σ)
    # sample y firstly
    dist_y = InverseGamma(ν/2, ν/2)
    ys = rand(dist_y, n)
    # sample x
    xs = σ * randn(n) .* sqrt.(ys) .+ μ
    return xs, ys
end

function cmp_res(n, μ, ν, σ)
    xs, ys = rt_dickey(n, μ, ν, σ)
    cum_δm = cumsum( exp.(- xs .^2) )
    cum_δm_star = cumsum( 1 ./ sqrt.(2 * σ^2 .* ys .+ 1) .* exp.(-μ^2 ./ (1 .+ 2*σ^2 .* ys)) )
    δm = cum_δm ./ (1:n)
    δm_star = cum_δm_star ./ (1:n)
    p = plot(δm, label = "MC")
    plot!(p, δm_star, label = "RB")
    return p
end

using Random
Random.seed!(123)
p1 = cmp_res(10000, 0, 4.6, 1)
p2 = cmp_res(10000, 3, 5, 0.5)
plot(p1, p2)
savefig("two-situations.svg")
```

![](https://666993855-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-LJfsESZOIJn_3uGIecs%2F-Luc4iIdHmCZZA-edwg8%2F-Luc4j5UP8jd25MwIYaz%2Ftwo-situations.svg?generation=1574778107808226\&alt=media)