10
$\begingroup$

How can I generate random numbers that follow a student-t distribution? From several sources I understand that this can be done using a random sample of size $n$ drawn from a normally distributed population, as follows:

$$ t = \frac {(x - m)} {(s/\sqrt n)} $$

Where $x$ is the sample mean, $m$ is the mean of the normal distribution (I assume you can just use the standard normal distribution, so $m=0$?), and $s$ is the sample standard deviation.

The degrees of freedom of the student-t distribution will then be $n-1$.

Do I understand correctly that in order to generate a random student-t value with $f$ degrees of freedom, I should first generate $f+1$ normally distributed values (i.e. standard normal), and then calculate the mean ($x$) and standard deviation ($s$) of these and apply the above formula? And if I do this repeatedly many times, the resulting random values will approach a student-t distribution with $f$ degrees of freedom?

I tried this in Excel using a macro that uses the above formula and another macro that generates random Gaussians (which works, I tested it) but the resulting random values do not seem to be completely student-t distributed. For instance with 6 degrees of freedom, the variance of 10,000 random values is about $1.7$ while it should be $6/(6-2) = 1.5$.

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ See Bailey, R.W. (1994). Polar generation of random variates with the t-distribution. Mathematics of Computation, 62(206), 779-781. $\endgroup$
    – Ray Koopman
    Sep 28, 2013 at 7:13

3 Answers 3

Reset to default
8
$\begingroup$

I have an answer to the practical part of your question, though not quite the theoretical one.

There is a function called TINV that directly does this. Except that it conly returns positive random t variates. You can get around that limitation with the following formula:

=TINV(RAND(),6)*(RANDBETWEEN(0,1)*2-1)

...you can replace 6 with whatever value you want for the DF, and the rand() can be replaced with any number between 0 and 1. The rest of it simply guarantees equal probability of negative and positive values.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks. This seems to work. I was only wondering if something like =TINV(RAND(),6)*IF(RAND()<0.5,1,-1) wouldn't be better. But either way the variance is 1.5. $\endgroup$
    – nominator
    Sep 17, 2013 at 17:45
  • 4
    $\begingroup$ This is a good Excel solution. I haven't checked recent upgrades of Excel, but previous versions had terrible implementations of TINV. They break down in the tails; it's really noticeable beyond (approximately) the 99.99 and 0.01 percentiles. Not a problem for quick-and-dirty work, but when you start using t-variates in simulations and for other calculations, the errors quickly become serious. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 18:06
  • $\begingroup$ Good to know. Hopefully everyone realizes that Excel is not statistical software and should only be used for quick and dirty work in the first place. :) $\endgroup$
    – f1r3br4nd
    Sep 18, 2013 at 10:54
  • $\begingroup$ @nominator , the two formulas should give the same results. My intuition tells me the first version would run faster, but I haven't actually tested this intuition. Speaking of runtimes, this thing will bog down hard if used as part of a large spreadsheet because every single instance of the formula will update whenever you make an edit to anything. So you might want it to be a standalone sheet from which you copy and paste as static values into you working sheet. ...or use R and presumably get better accuracy with rt(). $\endgroup$
    – f1r3br4nd
    Sep 18, 2013 at 11:02
10
$\begingroup$

Given a generator of i.i.d. standard gaussian random variates, you can generate $t_k$ distributed random variates (with any positive integer degree of freedom $k$) by using the relation:

$$Y=\frac{X_{k+1}}{\sqrt{k^{-1}\sum_{i=1}^k X_i^2}}$$

where $Y\sim t_k$ and $X_i\sim\text{i.i.d. }\mathcal{N}(0,1),i=1,\ldots,k+1.$

Improve this answer
$\endgroup$
10
  • $\begingroup$ I'm not sure what Xk+1 is supposed to be. Is it just one more standard gaussian random variate (which isn't one of the Xi in the denominator) ? or is it a mean? $\endgroup$
    – nominator
    Sep 17, 2013 at 17:51
  • 3
    $\begingroup$ You can do much better than this by generating one variate from a $\chi^2$ distribution for the denominator. Thus the t-variate can be realized by the ratio of a Normal to an independent $\chi^2$ value. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 18:04
  • $\begingroup$ @whuber But does generating one chi-squared variate not involve generating k gaussian variates? Or is there a faster way? $\endgroup$
    – nominator
    Sep 17, 2013 at 19:00
  • 4
    $\begingroup$ @nominator There are faster ways, including inverting the CDF and rejection sampling. Wikipedia remarks there is no universally fast way, but specialized methods have been developed for various ranges of the parameter. (Note that a chi-squared variate is a scaled Gamma variate.) For some of these methods visit the thread at stats.stackexchange.com/questions/7969. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 19:12
  • 1
    $\begingroup$ @user603 I suppose for proper scaling - to get Var(Y) = k/(k-2) - the summed squares in the denominator still need to be divided by k. $\endgroup$
    – nominator
    Sep 26, 2013 at 11:56
7
$\begingroup$

A fast way of generating a t variate, faster than the gaussian-only approach for all but the smallest degrees of freedom, is to use the fact that a t distribution is a mixture of Normals, with the mixture distribution being an inverted gamma distribution on the variance. Here's an example in R, where we generate 1,000,000 t(10) variates this way and compare to the theoretical distribution using the Kolmogorov-Smirnov test ("proof" by large experiment!):

> df <- 10
> s2 <- 1/rgamma(1000000, df/2, df/2)
> tv <- rnorm(1000000,0,sqrt(s2))
> 
> ks.test(tv, pt, df=df)

    One-sample Kolmogorov-Smirnov test

data:  tv 
D = 6e-04, p-value = 0.8826
alternative hypothesis: two-sided

This approach also works for non-integer degrees of freedom, which the gaussian-only approach does not.

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ +1 (clever use of a mixture indeed). But what "Gaussian-only approach" do you refer to? The one that divides a Gaussian by an independent chi-squared variate also works for nonintegral DF and is actually a little simpler to implement on many platforms. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 18:09
  • $\begingroup$ @whuber - the Gaussian-only approach was the method in the earlier answer, which used a $\chi^2$ variate generated by summing squared Gaussian variates. Of course, if you generate the $\chi^2$ variate using any of lots of other methods, you avoid the integral DF problem... now I'm wondering if there's a simple (fast) way of transforming a $\chi^2(1)$ variate to a $\chi^2(k<1)$ variate (w/o going through the probability integral transform)... $\endgroup$
    – jbowman
    Sep 17, 2013 at 19:29
  • 1
    $\begingroup$ I found this in my searches. This thread may be applicable, too: stats.stackexchange.com/questions/7969. $\endgroup$
    – whuber ♦
    Sep 17, 2013 at 19:32

Your Answer

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

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