Thanks for that! Both your solutions are MUCH quicker!
afaik overwork_ratio = 4/pi # what a hard way to calculate this :-)
The real problem I'm actually working on is a little more complicated :-) - 6 'random' variables, and 18 dependant variables (at last count)... currently only 4 computed terms in the 'if() c<-c+1' statement... I just created the x,y area problem as a simple example of the style I was using.
In the past I've worked on other 'accident' related problems and one difficulty was in getting enough final 'z' values to be statistically significant as there was an enormous rejection ratio - but it was running on a 486 written in compiled turbo pascal and it would take days of work...
I'll have a go at recoding using your techniques.
Many Thanks, Sean
Patrick Burns wrote:
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 [EMAIL PROTECTED] +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);
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
