I'm trying a simple Monte Carlo example which led to some confusion about ways of generating random Geometrically distributed values in R.
"A supercomputer is shared by 250 independent subscribers. Each day, each subscriber uses the facility with probability 0.3. The number of tasks sent by each active user has Geometric distribution with parameter 0.15, and each task takes a Gamma(10, 3) distributed computer time (in minutes). Tasks are processed consecutively. What is the probability that all the tasks will be processed, that is, the total requested computer time is less than 24 hours? Estimate this probability, attaining the margin of error ±0.01 with probability 0.99."
The solution is given in pure MATLAB and I will spare you everything unrelated to my question. The Geometric distribution is calculated as
Y=ceil( log(1-rand(X,1))/log(1-q) );
X is a Binomial variable (the number of active subscribers a given day) and
q is the parameter to the Geometric distribution. I translated this into R as
Y <- ceiling(log(1-runif(X))/log(1-q))
and got the expected result (~0.17) when I ran the simulation. However, as I understand it, R, unlike pure MATLAB, can generate Geometrically distributed random variables directly, like so
Y <- rgeom(X, q)
and when I tried that instead the simulation suddenly started returning values (~0.55).
I'm probably missing something obvious but can anyone explain what's going on here? I tried plotting the different distributions given by the expected value of X and they look similar enough but the results are obviously not equal. If any necessary details are missing tell me, and I will add them.