Hi Alex,
Michael is correct that vectorizing this is the way to go. Since this
is for class and vectorizing requires a bit of thought, I am not going
to demonstrate it with your function; however, here are some very
minor changes you can make that give fairly substantial performance
increases. It can be hard when you are learning create super
efficient code, but as a general rule of thumb (particularly when
loops are involved), byte compiling cannot hurt and may help quite a
bit, and it does not really require any extra thought or effort on
your part.
Cheers,
Josh (examples and benchmarks below)
## what you had just changing phat <- 0L
## and indenting the code and using <- instead of =
buffonslow <- function(n) {
x <- NULL
theta <- NULL
a <- NULL
phat <- 0L
i <- 1L
while (i <= n) {
x[i] <- runif(1, 0, 1/2)
theta[i] <- runif(1, 0, pi/2)
a[i] <- 1/2 * (sin(theta[i])) - x[i]
if (a[i] <= 0) phat <- phat + 1
else phat <- phat
i <- i + 1L
}
return(phat)
}
## basically what you had
## but instead of incrementally growing
## x, theta, and a (which requires R to make copies at each expansion
## instatiate them as n length vectors
buffon <- function(n) {
x <- vector("numeric", n)
theta <- vector("numeric", n)
a <- vector("numeric", n)
phat <- 0L
i <- 1L
while (i <= n) {
x[i] <- runif(1, 0, 1/2)
theta[i] <- runif(1, 0, pi/2)
a[i] <- 1/2 * (sin(theta[i])) - x[i]
if (a[i] <= 0) phat <- phat + 1
i <- i + 1L
}
return(phat)
}
## take the trivial step of byte compiling the function
require(compiler)
cbuffon <- cmpfun(buffon)
## benchmark slow and uncompiled versus compiled
## on 500 iterations
set.seed(1)
system.time(for (i in 1:500) buffonslow(i))
set.seed(1)
system.time(for (i in 1:500) buffon(i))
set.seed(1)
system.time(for (i in 1:500) cbuffon(i))
## but maybe you will not be doing many iterations
## just one large number
set.seed(1)
system.time(buffonslow(50000))
set.seed(1)
system.time(buffon(50000))
set.seed(1)
system.time(cbuffon(50000))
On my system, here are the times. Notice that just instatiating a
full sized vector yields a *huge* speedup when buffon is called for
large N. Compiling helps further.
> ## benchmark slow and uncompiled versus compiled
> ## on 500 iterations
> set.seed(1)
> system.time(for (i in 1:500) buffonslow(i))
user system elapsed
9.82 0.00 9.96
>
> set.seed(1)
> system.time(for (i in 1:500) buffon(i))
user system elapsed
8.06 0.01 8.38
>
> set.seed(1)
> system.time(for (i in 1:500) cbuffon(i))
user system elapsed
3.35 0.00 3.44
>
> ## but maybe you will not be doing many iterations
> ## just one large number
> set.seed(1)
> system.time(buffonslow(50000))
user system elapsed
89.63 0.29 93.70
> set.seed(1)
> system.time(buffon(50000))
user system elapsed
3.21 0.00 3.31
> set.seed(1)
> system.time(cbuffon(50000))
user system elapsed
1.36 0.00 1.39
On Sat, Nov 19, 2011 at 9:08 PM, R. Michael Weylandt
<michael.weyla...@gmail.com> <michael.weyla...@gmail.com> wrote:
> Also, vectorize this puppy - that's how you really get the best performance
> out of R, but in this case it will easily get rid of your problem.
>
> Michael
>
> On Nov 19, 2011, at 11:31 PM, Peter Langfelder <peter.langfel...@gmail.com>
> wrote:
>
>> Well, you assign numeric(0) or numeric(0)+1, which is still
>> numeric(0). No wonder the return value is always numeric(0).
>>
>> You probably need to replace numeric(0) simply by 0. Numeric(0) does
>> not mean 0, it means a numeric vector of length zero (i.e., empty).
>>
>> HTH,
>>
>> Peter
>>
>> On Sat, Nov 19, 2011 at 6:52 PM, alex_janssen <janss...@uoguelph.ca> wrote:
>>> Hi,
>>>
>>> I am trying to code buffons needle in R for a class
>>>
>>> This is my code w/ output from R, if anyone could tell me why this is
>>> happening it would be amazing,
>>>
>>> I can generate correct results without putting the steps into a function but
>>> alas that is the assignment.
>>>> buffon = function(n){
>>> + x = NULL
>>> + theta = NULL
>>> + a = NULL
>>> + phat = numeric(0)
>>> + i = 1
>>> +
>>> + while ( i <= n){
>>> + x[i] <- runif(1,0,1/2)
>>> + theta[i] <- runif(1,0,pi/2)
>>> + a[i] <- 1/2 * (sin(theta[i])) - x[i]
>>> +
>>> +
>>> + if (a[i] <= 0)phat <- phat + 1
>>> + else phat<-phat
>>> + i <- i + 1
>>> +
>>> + }
>>> + return(phat)
>>> + }
>>>> buffon(10)
>>> numeric(0)
>>>>
>>> Thanks Again
>>>
>>> the response should be the number of times a[i] is less than or equal to
>>> zero, which i am calling phat
>>>
>>> --
>>> View this message in context:
>>> http://r.789695.n4.nabble.com/Function-gives-numeric-0-for-every-input-tp4087834p4087834.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>>
>>> ______________________________________________
>>> 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.
>>>
>>
>> ______________________________________________
>> 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.
>
> ______________________________________________
> 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.
>
--
Joshua Wiley
Ph.D. Student, Health Psychology
Programmer Analyst II, ATS Statistical Consulting Group
University of California, Los Angeles
https://joshuawiley.com/
______________________________________________
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.
