Some more skilled folks can help with the curve fitting, but the general
answer is yes -- R will handle this quite ably.
Consider the following code:
<<-------------------------------------->>
n = 1e5
length = 1e5
R = matrix(sample(c(-1,1),length*n,replace=T),nrow=n)
R = apply(R,1,cumsum) ## this applies cumsum `row-wise' to R and will make
your life INFINITELY better
R = cbind(rep(0,n),R) ## Now each row is a random walk as you desired.
<<---------------------------------------->>
There are actually even faster ways to do what you are asking for, but this
introduces you to some useful R architecture, above all the apply function.
To see how long the longest stretch of negatives in each row is, the
following is a little sneaky but works pretty well:
countNegative = apply(R,1,function(x){which.max(table(cumsum(x>=0))})
then you can study these random numbers to do whatever with them.
The gist of this code is that it counts how many positive number have been
seen up to each point: for any stretch this doesn't increase, you must be
negative, so this identifies the longest such stretch on each row and
records the length. (It may be off by one so check it on a smaller R matrix.
An empirical confidence interval might be given by this
mu = mean(countNegative)
sig = sd(countNegative)/sqrt(length(countNegative))
So all together:
<<-------------------------------------->>
n = 1e3
length = 1e3
R = matrix(sample(c(-1,1),length*n,replace=T),nrow=n)
R = apply(R,1,cumsum) ## this applies cumsum `row-wise' to R and will make
your life INFINITELY better
R = cbind(rep(0,n),R) ## Now each row is a random walk as you desired.
fTemp <- function(x) {
return(max(table(cumsum(x>=0))))
}
countNegative = apply(R,1,fTemp)
mu = mean(countNegative)
sig = sd(countNegative)/sqrt(length(countNegative))
<<---------------------------------------->>
This runs pretty fast on my laptop, but you'll need to look into the
memory.limit() function if you want to increase the simulation parameters.
There are much faster ways to handle the simulation as well, but this should
get you off to a nice start with R.
Hope this helps,
Michael
On Wed, Jul 27, 2011 at 7:36 PM, Paul Menzel <
paulepan...@users.sourceforge.net> wrote:
> Dear R folks,
>
>
> I need to simulate first passage times for iterated partial sums. The
> related papers are for example [1][2].
>
> As a start I want to simulate how long a simple random walk stays
> negative, which should result that it behaves like n^(-½). My code looks
> like this.
>
> -------- 8< -------- code -------- >8 --------
> n = 100000 # number of simulations
>
> length = 100000 # length of iterated sum
> z = rep(0, times=length + 1)
>
> for (i in 1:n) {
> x = c(0, sign(rnorm(length)))
> s = cumsum(x)
> for (i in 1:length) {
> if (s[i] < 0 && s[i + 1] >= 0) {
> z[i] = z[i] + 1
> }
> }
> }
> plot(1:length(z), z/n)
> curve(x**(-0.5), add = TRUE)
> -------- 8< -------- code -------- >8 --------
>
> This code already runs for over half an hour on my system¹.
>
> Reading about the for loop [3] it says to try to avoid loops and I
> probably should use a matrix where every row is a sample.
>
> Now my first problem is that there is no matrix equivalent for
> `cumsum()`. Can I use matrices to avoid the for loop?
>
> My second question is, is R the right choice for such simulations? It
> would be great when R can also give me a confidence interval(?) and also
> try to fit a curve through the result and give me the rule of
> correspondence(?) [4]. Do you have pointers for those?
>
> I glanced at simFrame [5] and read `?simulate` but could not understand
> it right away and think that this might be overkill.
>
> Do you have any suggestions?
>
>
> Thanks,
>
> Paul
>
>
> ¹ AMD Athlon(tm) X2 Dual Core Processor BE-2350, 2,1 GHz
>
>
> [1] http://www-stat.stanford.edu/~amir/preprints/irw.ps
> [2] http://arxiv.org/abs/0911.5456
> [3]
> http://cran.r-project.org/doc/manuals/R-intro.html#Repetitive-execution
> [4] https://secure.wikimedia.org/wikipedia/en/wiki/Function_(mathematics)
> [5]
> http://finzi.psych.upenn.edu/R/library/simFrame/html/runSimulation.html
>
> ______________________________________________
> R-help@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.
>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
