R help - Aug 2008 - Vectorization of duration of the game in the gambler ruin's problem

Hey fellas:

In the context of the gambler's ruin problem, the following R code obtains
the mean duration of the game, in turns:

# total.capital is a constant, an arbitrary positive integer
# initial.capital is a constant, an arbitrary positive integer between, and not
including
# 0 and total.capital
# p is the probability of winning 1$ on each turn
# 1-p is the probability of loosing 1$
# N is a large integer representing the number of times to simulate
# dur is a vector containing the simulated game durations


T <- total.capital
dur <- NULL
for (n in 1:N) {
    x <- initial.capital
    d <- 0
    while ((x!=0)&(x!=T)) {
       x <- x+sample(c(-1,1),1,replace=TRUE,c(1-p,p))
       d <- d+1
    }
   dur <- c(dur,d)
}
mean(dur) #returns the mean duration of the game

The problem with this code is that, using the traditional control structures
(while, for, etc.) it is rather slow. Does anyone know of a way i could
vectorize the "while" and the "for" to produce a faster
code?

And while I'm at it, does anyone know of a discrete-event simulation package
in R such as the "SimPy" for Python?


Thanks in advance


	[[alternative HTML version deleted]]

Hi Jose,

If you are only interested in the expected duration, the problem can be solved
analytically - no simulation is needed.
Let P be the probability to get total.capital (and then 1-P is the probability
to loose all the money) when starting with initial.capital. This probability P
is well known (I do not remember it now but I can derive the formula if you need
- let me know). Let X(i) be the gain at game i and let D be the duration. Let
S(n) = X(1)+...+X(n).
Since EX(i) = p - (1-p) = 2p-1, S(n) - n*(2p-1) is a martingale, and since D is
a stopping time we get that E(S(D) - (2p-1)*D) = 0, so that (2p-1)*E(D) =
E(S(D)) = P*(total.capital-initial.capital) + (1-P)*(-initial.capital), and so
E(D) can be computed provided that p != 1/2.
If p = 1/2 then S(n) is a martingale and then by Wald's Lemma, E(S(D)^2) =
E(D)*E(X^2) = E(D). Since E(S(D)^2) = P*(total.capital-initial.capital)^2 +
(1-P)*(-initial.capital)^2, we can compute E(D).

Regards,

Moshe.

--- On Fri, 15/8/08, jose romero <jlaurentum at yahoo.com> wrote:
> From: jose romero <jlaurentum at yahoo.com>
> Subject: [R] Vectorization of duration of the game in the gambler
ruin's problem
> To: r-help at r-project.org
> Received: Friday, 15 August, 2008, 2:26 PM
> Hey fellas:
> 
> In the context of the gambler's ruin problem, the
> following R code obtains the mean duration of the game, in
> turns:
> 
> # total.capital is a constant, an arbitrary positive
> integer
> # initial.capital is a constant, an arbitrary positive
> integer between, and not including
> # 0 and total.capital
> # p is the probability of winning 1$ on each turn
> # 1-p is the probability of loosing 1$
> # N is a large integer representing the number of times to
> simulate
> # dur is a vector containing the simulated game durations
> 
> 
> T <- total.capital
> dur <- NULL
> for (n in 1:N) {
> ??? x <- initial.capital
> ??? d <- 0
> ??? while ((x!=0)&(x!=T)) {
> ?????? x <-
> x+sample(c(-1,1),1,replace=TRUE,c(1-p,p))
> ?????? d <- d+1
> ??? }
> ?? dur <- c(dur,d)
> }
> mean(dur) #returns the mean duration of the game
> 
> The problem with this code is that, using the traditional
> control structures (while, for, etc.) it is rather slow.
> Does anyone know of a way i could vectorize the
> "while" and the "for" to produce a
> faster code?
> 
> And while I'm at it, does anyone know of a
> discrete-event simulation package in R such as the
> "SimPy" for Python?
> 
> 
> Thanks in advance
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained,
> reproducible code.