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