4
$\begingroup$

I trying to sample using Gibbs from a proportional distribucion $f_{Z}(z)$: \begin{align*} f_{Z}(z) \propto e^{-z}\left(1-e^{-z}\right)^4, \quad z >0 \end{align*} using the joint $f_{Z,\textbf{U}}(z,\textbf{u})$ $$ f_{Z,\textbf{U}}(z,\textbf{u}) \propto e^{- z} \prod_{i=1}^4 I_{\left(e^{-z}, 1\right)}\left(u_i\right) $$ with $U \sim U(0,1)$. I found that $p(U_i|Z,U_{-i}) \sim U(e^{-z},1)$, but i cant find $p(Z|\textbf{U})$.How should I sample from $p(Z|\textbf{U})$? it would depend on all the $u_i$ values? I got stuck on this.

$\endgroup$
1
  • 2
    $\begingroup$ This is a slice sampler, the whole point of uniforms is to replace challenging parts of the target with indicators. In the current case, a single uniform would have been enough. $\endgroup$
    – Xi'an
    Dec 4, 2023 at 3:08

1 Answer 1

6
$\begingroup$

If $u_i$ has to be between $e^{-z}$ and $1$, then $e^{-z}$ needs to be between $0$ and $u_i$. Thus the distribution of $z$ given $\mathbf u$ is effectively a standard exponential RV with truncation:

\begin{align} p(z|\mathbf u) &\propto e^{-z}\prod_{i=1}^4I_{(0, u_i)}(e^{-z}) \\ &= e^{-z}\prod_{i=1}^4I_{(\log(u_i), \infty)}(z) \\ &= e^{-z}\prod_{i=1}^4I(z > -\log u_i) \\ &= e^{-z}I(z > -\log u_{(1)}) \end{align}

where $u_{(1)} = \min (u_1, u_2, u_3, u_4)$. The easiest way to sample from this distribution is with a rejection sampler, with steps:

  1. Sample $z^\star \sim Exp(1)$
  2. If $z^\star > -\log u_{(1)}$, then accept. Otherwise return to step 1.

Verification of algorithm:

enter image description here


R Code:

# GIBBS SAMPLER
# Number of samples
M <- 10000 
# Initialize values
z <- rep(NA, M)
u <- matrix(NA, nrow=M, ncol=4)
z[1]  <- rexp(1)
u[1,] <- rexp(4)
# Begin sampler
for(m in 2:M){
  umin <- min(u[m-1,])
  # Sample z | u
  while(is.na(z[m])){
    zstar <- rexp(1)
    if(zstar > -log(umin))
      z[m] <- zstar
  }
  # Sample u_i | u_-i, z
  for(i in 1:4){
    u[m,i] <- runif(1, exp(-z[m]), 1)
  }
}

# MAKE FIGURE
fz <- function(z) exp(-z)*(1-exp(-z))^4
const <- integrate(fz, lower=0, upper=100)$value

hist(z, freq=FALSE, breaks=30)
curve(fz(x)/const, add=TRUE, lwd=2, col='orange')

$\endgroup$
1
  • $\begingroup$ ooh I couldn't see that condition to accept, thank you very much! :D $\endgroup$
    – daniel
    Dec 4, 2023 at 1:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.