# [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

rho.f(rho = 0.3)
rho.f(seq(0, 1, by=.1))

even though
all.equal(0.3, seq(0, 1, by=.1)) == 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), 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), 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

```