If one has phyper(x,m,n,k) then pghyper(x,k,m,m+n).
[EMAIL PROTECTED] wrote:
Dear all,
it looks like rhyper() gives wrong results compared to theory and compared to sample() and rghyper(SuppDists).
Best regards
Jens
K <- 100 J <- 60 N <- K+J p <- K/N n <- 50
nn <- 100000
urn <- rep(0:1, c(J,K)) x <- sapply(1:nn, function(i){ sum(sample(urn, n)) }) y <- rhyper(nn, K,J, n)
require(SuppDists)
z <- rghyper(nn, a=K, k=n, N=N)
# hypergeometric mean and variance p*n p*(1-p)*n*(N-n)/(N-1)
# check we have parametrized ghyper correctly sghyper(a=K, k=n, N=N)[c("Mean", "Variance")]
var(x) var(y) # wrong var(z)
version
K <- 100 J <- 60 N <- K+J p <- K/N n <- 50
nn <- 100000
urn <- rep(0:1, c(J,K)) x <- sapply(1:nn, function(i){
+ sum(sample(urn, n)) + })
y <- rhyper(nn, K,J, n)
require(SuppDists)
[1] TRUE
z <- rghyper(nn, a=K, k=n, N=N)
# hypergeometric mean and variance p*n
[1] 31.25
p*(1-p)*n*(N-n)/(N-1)
[1] 8.107311
# check we have parametrized ghyper correctly sghyper(a=K, k=n, N=N)[c("Mean", "Variance")]
$Mean [1] 31.25
$Variance [1] 8.107311
var(x)
[1] 8.060095
var(y) # wrong
[1] 8.61289
var(z)
[1] 8.150472
version
_ platform i386-pc-mingw32
arch i386 os mingw32 system i386, mingw32 status major 2 minor 0.0 year 2004 month 10 day 04 language R
--
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
-- Bob Wheeler --- http://www.bobwheeler.com/ ECHIP, Inc. --- Randomness comes in bunches.
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
