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 ______________________________________________ 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.
