full precision may not have been achieved in '

问题:

I am running R version 2.14.0 on a PC which uses Windows 7 Ultimate (Intel Core i5-2400 3GHz processor with 8.00GB ram). Let me know if other specs needed.

I am trying to simulate correlated beta distributed data. The method I am using is an extension of what is written in this paper:

http://onlinelibrary.wiley.com/doi/10.1002/asmb.901/pdf

1. Basically, I start by simulating multivariate normal data (using mvrnorm() function from MASS).
2. Then I use pnorm() to apply the probit transform to these data such that my new data vector(s) live on (0,1). And are still correlated according to the previous statement.
3. Then given these probit transformed data I apply the qbeta() function with certain shape1 and shape2 parameters, to get back correlated beta data with certain mean and dispersion properties.

I know other methods for generating correlated beta data exist. I am interested in why qbeta() causes this method to fail for certain "seeds". Below is the error message I get.

Warning message:
In qbeta(probit_y0, shape1 = a0, shape2 = b0) :
  full precision may not have been achieved in 'qbeta'

What does this mean? How can it be avoided? When it does occur within the context of a larger simulation what is the best way to ensure that this problem does not terminate the entire sourced (using source()) simulation code?

I ran the following code for integer seeds from 1:1000. Seed=899 was the only value which gave me problems. Though if its problematic here, it inevitably will be problematic for other seeds too.

library(MASS)
set.seed(899)
n0 <- 25  
n1 <- 25    
a0 <- 0.25    
b0 <- 4.75    
a1 <- 0.25    
b1 <- 4.75    
varcov_mat <- matrix(rep(0.25,n0*n0),ncol=n0)
diag(varcov_mat) <- 1
y0 <- mvrnorm(1,mu=rep(0,n0),Sigma=varcov_mat)
y1 <- mvrnorm(1,mu=rep(0,n1),Sigma=varcov_mat)
probit_y0 <- pnorm(y0)
probit_y1 <- pnorm(y1)
beta_y0 <- qbeta(probit_y0, shape1=a0, shape2=b0)
beta_y1 <- qbeta(probit_y1, shape1=a1, shape2=b1)

The above code is a fragment of a piece of a larger simulation project. But the qbeta() warning message is what is giving me a headache now.

Any help the group could provide would be greatly appreciated.

Cheers Chris

回答1:

The reason for the error is that the algorithm used to calculate qbeta did not converge for those values of the parameters.

R uses AS 109 to calculate qbeta (Cran, G. W., K. J. Martin and G. E. Thomas (1977). Remark AS R19 and Algorithm AS 109, Applied Statistics, 26, 111–114, and subsequent remarks (AS83 and correction).). R tries to calculate the value in 1000 iterations. If it cannot in the 1000 iterations, you get the error message you saw.

Here is the code for qbeta.