\begin{Schunk}
\begin{Sinput}
 par(mfrow=c(2,2))
 x <- seq(0.00001, 6, length=1500)
 y1 <- dgamma(x, shape = 1/3, scale = 3)
 plot(x, y1, main = "shape = 1/3", type = "l", ylim = c(0,2), 
       xlab = expression(delta[i]),
       ylab = expression(density(delta[i])))
 y2 <- dgamma(x, shape = 2/3, scale = 3/2)
 plot(x, y2, main = "shape = 2/3", type = "l", ylim = c(0,2), 
       xlab = expression(delta[i]),
       ylab = expression(density(delta[i])))
 y3 <- dgamma(x, shape = 1/1, scale = 1)
 plot(x, y3, main="shape = 1", type = "l", ylim = c(0,2), 
       xlab = expression(delta[i]),
       ylab = expression(density(delta[i])))
 y4 <- dgamma(x, shape = 3/2, scale = 2/3)
 plot(x, y4, main="shape = 1.5", type = "l", ylim = c(0,2), 
       xlab = expression(delta[i]),
       ylab = expression(density(delta[i])))
 par(mfrow=c(1,1))
\end{Sinput}
\end{Schunk}
