[R] A coding question
Uwe Ligges
ligges at statistik.uni-dortmund.de
Sat Jun 3 12:23:10 CEST 2006
xpRt.wannabe wrote:
> Uwe and Ben,
>
> Thank you both for your help.
>
> To me, both sets of code seem to do the job and should produce the
> same results. However, as a test I inserted set.seed( ) as follows.
> Unless I put set.seed( ) in the wrong lines, the results produced by
> both sets of code turn out to be different. I am baffled.
> y <- replicate(10, {
> set.seed(123)
So you reset the RNG 10 times, makes no sense, move set.seed before the
first replicate().
> rp <- rpois(8, 5)
You get 8 (!) different values here (in contrast to the code below).
> mysum <- sapply(rp, function(x) {
> set.seed(123)
And again. Makes no sense either.
> x <- rnorm(x)
> x <- x - max(0, x-15) + max(0, x-90)
> sum(x)
> })
> })
>
> ##
>
> tmpf <- function() {
> set.seed(123)
Here, you reset 80 times, each time you get the same value for rpois().
> x <- rnorm(rpois(1,5))
> sum(ifelse(x>90,x-75,pmin(x,15)))
> }
> replicate(10,replicate(8,tmpf()))
You cannot compare the results of both code fragments with set.seed()
due to the different stucture when which random number is generated.
(even permutation of generation of normal and poisson distributed random
numbers).
Uwe Ligges
> On 6/2/06, Uwe Ligges <ligges at statistik.uni-dortmund.de> wrote:
>
>> xpRt.wannabe wrote:
>>
>> > Dear List:
>> >
>> > I have the follow code:
>> >
>> > y <- replicate(10,replicate(8,sum(rnorm(rpois(1,5)))))
>> >
>> > Now I need to apply the following condition to _every_ randomly
>> generated
>> > Normal number in the code above:
>> >
>> > x - max(0,x-15) + max(0,x-90), where x represents the individual Normal
>> > numbers.
>> >
>> > In other words, the said condition needs to be applied before
>> > replicate(...(replicate(...(sum(...))) takes place.
>> >
>> > Any help would be greatly appreciated.
>>
>>
>> y <- replicate(10, {
>> rp <- rpois(8, 5)
>> mysum <- sapply(rp, function(x) {
>> x <- rnorm(x)
>> x <- x - max(0, x-15) + max(0, x-90)
>> sum(x)
>> })
>> })
More information about the R-help
mailing list