Hi all,
Sorry I have too many questions. I could not think of a way to fix
problem. Can anyone give some suggestions on fixing this?
Hannah
2010/5/28 li li <hannah....@gmail.com>
> Hi Sarah,
> Thanks for your kind help. I now know where the problem is.
> Hannah
>
> 2010/5/28 Sarah Goslee <sarah.gos...@gmail.com>
>
>> My initial guess appears to be right: you're working with something
>>
>> exceedingly sensitive to floating point precision. You may have to
>> reconsider your methods.
>>
>> Your problem is:
>>
>> rho.f(rho = 0.3)
>> gives a different answer than
>> rho.f(seq(0, 1, by=.1)[4])
>>
>> even though
>> all.equal(0.3, seq(0, 1, by=.1)[4]) == TRUE
>>
>>
>> The only thing rho.f does with rho is passes it to rdata1.mnorm
>> so the rest of the function is irrelevant for this question.
>>
>> The only thing rdata1.mnorm does with rho is passes it to var.f
>>
>> So we've already ruled out most of your code and narrowed the problem
>> down to one small section.
>>
>> Take a look at this example and run it for yourself:
>>
>> # starting parameters
>> m=30;n=10;mzero=5;mu0=0;mu1=2
>> mean1 <- c(rep(mu0, mzero), rep(mu1, m-mzero))
>>
>> ## make two different var.f() results
>>
>> varf1 <- var.f(0.3, m, rep(1, m))
>> varf2 <- var.f(seq(0, 1, by=0.1)[4], m, rep(1, m))
>>
>> ## compare them
>>
>> all.equal(varf1, varf2)
>> all(varf1 == varf2)
>>
>> ## here's where the problem is
>> ## The function you're calling is extremely sensitive, and 0.3
>> ## cannot be represented exactly.
>>
>> set.seed(103)
>> mvrnorm(n, mean1, varf1)
>>
>> set.seed(103)
>> mvrnorm(n, mean1, varf2)
>>
>> ## Here's a check on that idea:
>>
>> set.seed(103)
>> mvrnorm(n, mean1, round(varf1, 1))
>>
>> set.seed(103)
>> mvrnorm(n, mean1, round(varf2, 1))
>>
>>
>> #### and compare to this
>>
>>
>> ## make two different var.f() results
>>
>> varf1 <- var.f(0.2, m, rep(1, m))
>> varf2 <- var.f(seq(0, 1, by=0.1)[3], m, rep(1, m))
>>
>> ## compare them
>>
>> all.equal(varf1, varf2)
>> all(varf1 == varf2)
>>
>> ## Look: 0.2 can be represented exactly.
>>
>> set.seed(103)
>> mvrnorm(n, mean1, varf1)
>>
>> set.seed(103)
>> mvrnorm(n, mean1, varf2)
>>
>>
>> --
>> Sarah Goslee
>> http://www.functionaldiversity.org
>>
>
>
[[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.
