[Rd] rhyper has too high variance (PR#7314)

[Bob Wheeler]

If one has phyper(x,m,n,k) then pghyper(x,k,m,m+n).

oehl_list at gmx.de 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
> 
> 
> 
> 
> 
> 
> --
> 
> ______________________________________________
> R-devel at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 


-- 
Bob Wheeler --- http://www.bobwheeler.com/
         ECHIP, Inc. ---
Randomness comes in bunches.