df=5
x2=rgamma(10000, shape=df/2, rate=df/2)
x1=rnorm(10000, mean = 0, sd=sqrt(1/x2))

library(MASS)
truehist(x1, col="white")
curve(dt(x, df=5), from=-5, to=5, add=T, col=2, lwd=2)