如何处理极大数的算术运算？

Question

我试图找到特定类型订单统计的累积分布函数；逐步审查的统一订单统计数据。我生成了以下

R

代码：

#The CDF for rth order progressively censored uniform order statistics
n <- 30 #Total number of experimental units
m <- 15 #Desired numbers of failure
R <- c(rep(0, m - 1), n - m) #Progressive censoring scheme
order <- m #Order of the censored order statistics (Here, maximum order)

gam <- NA
Cr <- NA
for(i in 1 : order)
{
  gam[i] <- m - i + 1 + sum(R[i : m])
}
for(i in 1 : order)
{
  Cr[i] <- prod(gam[1 : i])
}
air <- array(dim = c(order, order))

for(i in 1 : order)
{
  for (j in 1 : order) {
    
    if(i != j)
    {
      air[i, j] <- 1/(gam[j] - gam[i])
    }
  }
}
A <- NA
for(i in 1 : order)
{
  A[i] = prod(na.omit(air[i,]))
}
#CDF of progressively censored uniform order statistics
progU_CDF <- function(u)
{
  CDF = NA
  for(i in 1 : length(u))
  {
    CDF[i] <- 1 - (Cr[order] * sum((A/gam) * ((1 - u[i])^(gam))))
  }
  return(CDF)
}

现在，

progU_CDF(0)

应该给出 0，而

progU_CDF(1)

应该给出 1。尽管在这种情况下，

progU_CDF(0)

产生的数字非常接近 0，而不是完全 0。从数学上来说，这个想法是，当

u = 0

时，

 Cr[order] * sum(A/gam) = 1

。

此外，当我绘制 CDF 时，形状是理想的，即单调非递减。

plot(seq(0, 1, 0.01), progU_CDF(seq(0, 1, 0.01)), type = "l")

阶次统计 CDF 图 n = 30，m = 15

然而，当我服用

n = 50

和

m = 25

时，事情开始变得疯狂。没有任何改变，但当

Cr[order] * sum(A/gam)

时，

u = 0

远不接近 1。 CDF 图如下所示：

阶次统计 CDF 图 n = 50，m = 25

我怀疑这是由于对极大数进行算术运算造成的。但我无法追踪它。

更令人困惑的是

Cr[order] * sum(A/gam)

和

sum(Cr[order] * A/gam)

产生两个不同的数字，这是违反直觉的，因为

Cr[order]

是一个常数。

我的问题是，为什么它对

n = 30, m = 15

有效，但对

n = 50, m = 25

无效？有没有办法处理这么大的数字，使得无论

Cr[order] * sum(A/gam)

和

u = 0

的值如何，每当

n

时，

m

始终接近 1？

Answer 1

用

Rmpfr

进行矢量化：

library(Rmpfr)

fprogU_CDF <- function(n, m, precBits = 128) {
  R <- c(numeric(m - 1), n - m) #Progressive censoring scheme
  order <- m #Order of the censored order statistics (Here, maximum order)
  
  gam <- mpfr(rev(cumsum(R[m:order])) + m:(m - order + 1), precBits)
  Cr <- cumprod(gam[1:order])
  air <- 1/outer(gam, gam, "-")
  diag(air) <- 1
  A <- apply(air, 1, prod)
  diag(air) <- NA
  #CDF of progressively censored uniform order statistics
  function(u) {
    as.numeric(1 - Cr[order]*colSums(A/gam*outer(gam, 1 - u, \(x, y) y^x)))
  }
}

使用

n = 30

和

m = 15

进行测试。

progU_CDF <- fprogU_CDF(30, 15)
progU_CDF(0)
#> [1] 2.15461e-27
progU_CDF(1)
#> [1] 1
curve(progU_CDF(x), 0, 1)

使用

n = 50

和

m = 25

进行测试。

progU_CDF <- fprogU_CDF(50, 25, 1024)
progU_CDF(0)
#> [1] 4.050669e-288
progU_CDF(1)
#> [1] 1
curve(progU_CDF(x), 0, 1)