## student t random number generator
## using data augmentation method

n<-10000

df = 5
mu =0
sig = 1

w<-rep(NA, n)
t.y<-rep(NA, n)

for (i in 1:n){
	#w[i]<-rchisq(1,df)
	w[i]<-rgamma(1, df/2, df/2)
	t.y[i]<-rnorm(1, mu, sig/sqrt(w[i]))
	}
	
## Q-Q plot
prob<-seq(.01, .99, .01)

plot(qt(prob, df), quantile(t.y, prob))

qqplot(t.y, rt(n, df))