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.
