[R] problem with a function
Sarah Goslee
sarah.goslee at gmail.com
Fri May 28 19:58:40 CEST 2010
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
More information about the R-help
mailing list