[R] problem with a function

[Sarah Goslee]

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