Hi all,
    Sorry I have too many questions. I could not think of a way to fix
problem.   Can anyone give some suggestions on fixing this?
                    Hannah

2010/5/28 li li <hannah.hlx@gmail.com>

> Hi Sarah,
>     Thanks for your kind help. I now know where the problem is.
>                         Hannah
>
>  2010/5/28 Sarah Goslee <sarah.goslee@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]]

