Generating Geometric variates without using inbuilt function in R

user_8484779 0
Viewed 233

I am learning R and have just started simulating from different probability distributions. I was trying to generate a random variate from the Geometric Distribution. I am not willing to use the inbuilt function rgeom in R for doing the same. Among other methods, I decided to use the result that:

The Geometric distribution is the probability distribution of the number of Bernoulli trials needed to get one success.

I wrote the function for generating one Bernoulli variable as in my previous post named as ber, with p being the success probability.

Now to generate a Geometric variable, I can run a loop (I used a repeat loop here) and continuously generate a Bernoulli variable, check if it is a success or not and so on...

I did it like the following:

geom<-function(p)
{
    x<-1
    repeat
    {
        if(ber(p)==1) break
        x<-x+1
    }
    return(x)
}

One output as generated by R is:

geom(0.035)
[1] 9

What I want to do and where am I stuck:

I want to use a vectorized approach instead of this repeat loop method to generate a Geometric variable but I seem to have run out of ideas as to how that can be achieved. Any helps/suggestions and hints are very much welcome in here.

Note 1: As to why I defined a new function geom(p), I originally planned to generate not one, but many (say n) Geometric variates. I figured that defining a function as like this, and then using the replicate function to generate more similar variables with the same success probability is easier.

Note 2: I really apologize if some of my previous posts seem similar (in the sense that asking for suggestion for replacing a loop with some vector methods), but I am a newbie in R and I am still learning about vectorized methods.

r
1 Answers

Without using a loop, you could use the fact that geometric random variables have an invertible cumulative distribution function and use inverse transform sampling. The CDF of your version of the geometric variable (number of first success rather than the version where you count the number of failures before the first success) is:

F(x) = 1 - (1-p)^x,  x = 1, 2, 3, ....

We can solve for x in this:

x = log(1-F(x))/(log(1-p))

To use the with inverse sampling we need to take the ceiling of this to make it discrete:

geom <- function(n,p) {
  ceiling(log(runif(n))/log(1-p))
}

The above code exploits the fact that if U is uniform in [0,1] then so is 1-U.

For example:

> set.seed(42)
> geom(10,1/6)
 [1]  1  1  7  2  3  4  2 11  3  2
> mean(geom(100000,1/6))
[1] 6.00359

On Edit An approach which is more along the lines of what you currently have is to create an entire vector of random variables, initialized all to 0, and keep a list of indices of variables still "active" in the sense of experiments still waiting for a success. At each step, add 1 to these active variables, and with probability 1-p keep some indices active into the next step:

geom2 <- function(n,p) {
  nums <- numeric(n)
  active <- 1:n
  k <- n
  while(k > 0){
    nums[active] <- nums[active]+1
    active <- active[runif(k) <= 1-p]
    k = length(active)
  }
  nums
}

This is still reasonably fast, though once n gets into the millions it visibly lags in comparison with the first definition.

I don't see any way to avoid using either a loop or an inverse transform.

user_49962480
Related