```julia
using Random
using Distributions
using Statistics
using Printf
using Plots
```


```julia
# for readability, these are good settings to use
default(xtickfontsize=14,  ytickfontsize=14, ztickfontsize=14,
    guidefontsize=14, legendfontsize=12, lw=2,ms=8)
```

## Example 6.5
Taken from Rubinstein and Kroese.  Sample from the following density in $\mathbb{R}^2$ using the Gibbs sampler:
$$
f(x,y) = c_1(y) \exp \left(-\frac{1+y^2}{2} \left(x - \frac{4}{1+y^2}\right)^2\right)
$$
After some algebra, it can be shown that
\begin{align}
X\mid Y = y &\sim N\left(\tfrac{4}{1+y^2}, \tfrac{1}{1+y^2}\right)\\
Y\mid X = x &\sim N\left(\tfrac{4}{1+x^2}, \tfrac{1}{1+x^2}\right)
\end{align}
allowing us to sample basic Gaussians at each step.



```julia
X0 = ones(2);
T = 10^6; # total number of iterations
```


```julia
X = deepcopy(X0);
X_samples = [deepcopy(X)];
Random.seed!(100);
for _ in 1:T-1
    x = rand(Normal(4/(1+X[2]^2), 1/sqrt(1+X[2]^2)));
    y = rand(Normal(4/(1+x^2), 1/sqrt(1+x^2)));
    @. X  = [x,y];
    push!(X_samples, deepcopy(X))    
end
```


```julia
histogram2d([X_[1] for X_ in X_samples],[X_[2] for X_ in X_samples], norm =:pdf)
xlabel!("x")
ylabel!("y")
```


```julia
@show mean([X_[1] for X_ in X_samples]);
```


```julia
plot(1:T, cumsum([X_[1] for X_ in X_samples])./(1:T), label="Running Avg.", legend=:bottomright, xscale=:log10)
xlabel!("Iterate")
ylabel!("Estimate")
```

121feb633e9b0a64822051114aecf4a5634ca7fe
ipynb
Jupyter Notebook
notebooks/Gibbs Sampler.ipynb
gideonsimpson/Math_483_2022
9f6d5c7411971ad0d6f74d4db9f270a1809d5069
MIT
notebooks/Gibbs Sampler.ipynb
gideonsimpson/Math_483_2022
9f6d5c7411971ad0d6f74d4db9f270a1809d5069
MIT
notebooks/Gibbs Sampler.ipynb
gideonsimpson/Math_483_2022
9f6d5c7411971ad0d6f74d4db9f270a1809d5069
MIT
Qwen/Qwen-72B
1. YES
2. YES