Saturday, May 15, 2010

Simulating from a T1EV (Gumbel) Distribution

I have recently starting writing my own functions in R, and this has been very helpful in getting some of my coding projects done. Because R allows you to write functions, you can use this feature to greatly augment R's already-powerful capabilities.

For example, I recently needed to simulate 1.5 million independent draws from a Type 1 Extreme Value distribution, which is known on Wikipedia as the Gumbel distribution. I tried using rgumbel(), but no such function exists in the standard packages. I also conducted a quick search of R help files... to no avail.

What's an R user to do? Write a function!

Fortunately, I remembered this really handy theorem that transforming any continuous random variable by its CDF F() always produces a random variable that is uniformly distributed on the [0,1] interval. As a practical matter, this means that you can simulate n random draws from a distribution by:
1. Computing n random draws from a Uniform distribution on [0,1]
2. Applying the inverse function of the CDF to the uniformly draw random variables.
I used this algorithm and the inverse formula for the Gumbel CDF to write the following function.

rt1ev = function(n){ ## Function that returns n
U = runif(n) ## independent draws from T1EV
t1ev = -log(-log(U))
t1ev
}

Because IO economists are fond of using the T1EV distribution and I am studying IO now, I can see myself using this function repeatedly.

Here's how it works. Suppose you want a vector x to have 100 draws from a T1EV distribution. After loading my code, all you need to do to use this function is type:

x = rt1ev(100)

Not only is this a useful function if you want to simulate from a Gumbel distribution, but it is a great illustration of why knowing how to write functions in R can make your life easier.