Hi Folks,
I'm trying to learn R. One of my intentions is to do some Monte-Carlo
type modelling of road "accidents".
Below, to simplify things, I've appended a little program which does a
'monte-carlo' type simulation. However, it is written in a way which
seems a bit un-natural in R. Could someone help me make this a bit more
R-ish please?
Or is there a completely different approach I should be taking?
Many thanks in advance,
Sean O'Riordain
seanpor AT acm.org
--------------------------------------------
n <- 900; # number of valid items required...
x <- numeric(n);
y <- numeric(n);
z <- numeric(n);
c <- 1; # current 'array' pointer
tc <- 0; # total items actually looked at...
while (c <= n) {
x[c] = runif(1, 0, 1);
y[c] = runif(1, 0, 1);
z[c] = sqrt(x[c]^2 + y[c]^2);
if (z[c] < 1)
c <- c + 1;
tc <- tc + 1;
}
print("'overwork' ratio");
print(tc/(c-1));
plot(x,y);
1. I suggest you avoid using "c" as a loop index, as it conflicts
with the name of a function. R is smart enough to figure out the
difference in some cases but not in all.
2. How about the following:
n <- 900
x <- runif(n)
y <- runif(n)
z <- sqrt(x^2+y^2)
print(list(Overwork.ratio = n/sum(z<1)))
plot(x, y)
theta <- seq(0, pi/2, length=31)
lines(sin(theta), cos(theta))
######################
This won't give you n numbers z < 1. If you need exactly that,
what about first generating, say, 2*n, then either throwing away the
excess or generating another 2*n if you don't have enough? You can use
a recursive function for this, which would almost never recurse with 2*n
but might with 1.1*n.
hope this helps. spencer graves
Sean O'Riordain wrote:
> Hi Folks,
>
> I'm trying to learn R. One of my intentions is to do some Monte-Carlo
> type modelling of road "accidents".
>
> Below, to simplify things, I've appended a little program which does a
> 'monte-carlo' type simulation. However, it is written in a way
which
> seems a bit un-natural in R. Could someone help me make this a bit
> more R-ish please?
>
> Or is there a completely different approach I should be taking?
>
> Many thanks in advance,
> Sean O'Riordain
> seanpor AT acm.org
>
>
> --------------------------------------------
> n <- 900; # number of valid items required...
>
> x <- numeric(n);
> y <- numeric(n);
> z <- numeric(n);
>
> c <- 1; # current 'array' pointer
> tc <- 0; # total items actually looked at...
>
> while (c <= n) {
> x[c] = runif(1, 0, 1);
> y[c] = runif(1, 0, 1);
>
> z[c] = sqrt(x[c]^2 + y[c]^2);
> if (z[c] < 1)
> c <- c + 1;
> tc <- tc + 1;
> }
>
> print("'overwork' ratio");
> print(tc/(c-1));
> plot(x,y);
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
For this particular problem you can probably use
polar coordinates.
But something similar to your code could be:
x <- runif(900)
y <- runif(900)
z <- sqrt(x^2 + y^2)
okay <- z < 1
while(any(!okay)) {
n.bad <- sum(!okay)
x[!okay] <- runif(n.bad)
y[!okay] <- runif(n.bad)
z <- sqrt(x^2 + y^2) # restricting to !okay may or may not be useful
okay <- z < 1
}
Patrick Burns
Burns Statistics
patrick at burns-stat.com
+44 (0)20 8525 0696
http://www.burns-stat.com
(home of S Poetry and "A Guide for the Unwilling S User")
Sean O'Riordain wrote:
> Hi Folks,
>
> I'm trying to learn R. One of my intentions is to do some Monte-Carlo
> type modelling of road "accidents".
>
> Below, to simplify things, I've appended a little program which does a
> 'monte-carlo' type simulation. However, it is written in a way
which
> seems a bit un-natural in R. Could someone help me make this a bit
> more R-ish please?
>
> Or is there a completely different approach I should be taking?
>
> Many thanks in advance,
> Sean O'Riordain
> seanpor AT acm.org
>
>
> --------------------------------------------
> n <- 900; # number of valid items required...
>
> x <- numeric(n);
> y <- numeric(n);
> z <- numeric(n);
>
> c <- 1; # current 'array' pointer
> tc <- 0; # total items actually looked at...
>
> while (c <= n) {
> x[c] = runif(1, 0, 1);
> y[c] = runif(1, 0, 1);
>
> z[c] = sqrt(x[c]^2 + y[c]^2);
> if (z[c] < 1)
> c <- c + 1;
> tc <- tc + 1;
> }
>
> print("'overwork' ratio");
> print(tc/(c-1));
> plot(x,y);
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>
>
On Wed, Oct 15, 2003 at 10:06:36AM +0100, Sean O'Riordain wrote:> n <- 900; # number of valid items required... > > x <- numeric(n); > y <- numeric(n); > z <- numeric(n); > > c <- 1; # current 'array' pointer > tc <- 0; # total items actually looked at... > > while (c <= n) { > x[c] = runif(1, 0, 1); > y[c] = runif(1, 0, 1); > > z[c] = sqrt(x[c]^2 + y[c]^2); > if (z[c] < 1) > c <- c + 1; > tc <- tc + 1; > } > > print("'overwork' ratio"); > print(tc/(c-1)); > plot(x,y);If I read this correctly you want to find the frequency of the event "sqrt(x^2 + y^2) < 1" where 0 <= x, y <= 1 right? How about this: n <- 1000 z <- sqrt(runif(n,0,1)^2 + runif(n,0,1)^2) ratio <- (length(z)-1) / (length( z[z<1] )) The main difference of course is that I uses a fixed number of "total items" rather than "valid items". If that's a problem and you need exactly 900 "valid items" things get a bit more complicated... cu Philipp -- Dr. Philipp Pagel Tel. +49-89-3187-3675 Institute for Bioinformatics / MIPS Fax. +49-89-3187-3585 GSF - National Research Center for Environment and Health Ingolstaedter Landstrasse 1 85764 Neuherberg, Germany