Untitled

```
np.random.seed(410)
n = 30
x = stats.norm(0, 1).rvs(n)

C, G = 2, 1000
theta = np.zeros((C, G))
tau   = np.ones((C, G))  # initialize tau to 1

for c in range(C):
    for g in range(1, G):
        tau_prop = stats.gamma(a=4, scale=0.5).rvs()
        log_like_prop = np.sum(stats.norm(loc=theta[c, g-1],
                                          scale=tau_prop).logpdf(x))
        log_like_curr = np.sum(stats.norm(loc=theta[c, g-1],
                                          scale=tau[c, g-1]).logpdf(x))
        
        # Prior: half-normal with scale 3 for tau:
        log_prior_prop = stats.halfnorm(scale=3).logpdf(tau_prop)
        log_prior_curr = stats.halfnorm(scale=3).logpdf(tau[c, g-1])
        
        # Proposal density evaluations: note these are for i.i.d. samples so
        # q(·) = Gamma(4, scale=0.5).pdf(·)
        log_q_prop = stats.gamma(a=4, scale=0.5).logpdf(tau_prop)
        log_q_curr = stats.gamma(a=4, scale=0.5).logpdf(tau[c, g-1])
        
        log_r_tau = (log_like_prop + log_prior_prop + log_q_curr) - (log_like_curr + log_prior_curr + log_q_prop)
        
        if np.log(np.random.rand()) < min(0, log_r_tau):
            tau[c, g] = tau_prop
        else:
            tau[c, g] = tau[c, g-1]
        
        theta_prop = stats.norm(loc=theta[c, g-1], scale=0.25).rvs()
        
        # Compute the log likelihoods with the updated tau[c, g]:
        log_like_prop = np.sum(stats.norm(loc=theta_prop,
                                          scale=tau[c, g]).logpdf(x))
        log_like_curr = np.sum(stats.norm(loc=theta[c, g-1],
                                          scale=tau[c, g]).logpdf(x))
        
        # Prior for theta is a t-distribution with 2 degrees of freedom:
        log_prior_prop = stats.t(df=2).logpdf(theta_prop)
        log_prior_curr = stats.t(df=2).logpdf(theta[c, g-1])
        
        # For theta the proposal is symmetric, so the q's cancel:
        log_r_theta = (log_like_prop + log_prior_prop) - (log_like_curr + log_prior_curr)
        
        if np.log(np.random.rand()) < min(0, log_r_theta):
            theta[c, g] = theta_prop
        else:
            theta[c, g] = theta[c, g-1]
```