Kolomogorov-Smirnov检验：C到R翻译问题

Question

Kolomogorov-Smirnov检验：C到R翻译问题

3

我在将一个C语言算法翻译为R语言时遇到了困难，这个算法与Kolmogorov-Smirnov检验有关，更具体地说是KS概率函数。

在《C语言经典数值计算方法》中，'probks'的代码如下：

#include <math.h>
#define EPS1 0.001
#define EPS2 1.0e-8
float probks(float alam)
/*Kolmogorov-Smirnov probability function.*/
{
   int j;
   float a2,fac=2.0,sum=0.0,term,termbf=0.0;

   a2 = -2.0*alam*alam;
   for (j=1;j<=100;j++) {
      term=fac*exp(a2*j*j);
      sum += term;
      if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum) return sum;
      fac = -fac; /*Alternating signs in sum.*/
      termbf=fabs(term);
   }
   return 1.0; /* Get here only by failing to converge. */
}

我不知道如何在R中处理最后几行的翻译，现在我只有这些

PROBKS <- function(lambda) {

  EPS1 <- 0.001; EPS2 <- 1.0e-8;
  sum <- 0.0; fac <- 2.0; termbf <- 0.0; 
  a2 <- -2*lambda*lambda 

  for (j in 1:100) {
    term <- fac * exp(a2*j*j)
    sum <- sum + term
    if ( (abs(term) <= EPS1*termbf) || (abs(term) <= EPS2*sum) ) {
      break
    } else {
      fac <- -fac
    }
  }
  termbf <- abs(term)
  return(sum)
}

但是这会产生一个非单调的概率函数

其中应该是$Q_KS(0)=1$和$Q_KS(\infty)=0$。显然，这与如何解释/编码最后的'if'语句有关。

任何帮助将非常感激。 M

编辑1： 下面是我的会话信息

> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252 
[2] LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] reshape2_1.4.3  forcats_0.3.0   stringr_1.3.1   dplyr_0.7.7    
 [5] purrr_0.2.5     readr_1.1.1     tidyr_0.8.1     tibble_1.4.2   
 [9] ggplot2_3.1.0   tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] withr_2.1.2      rvest_0.3.2      tidyselect_0.2.5 lattice_0.20-35 
 [5] pkgconfig_2.0.2  xml2_1.2.0       compiler_3.4.4   readxl_1.1.0    
 [9] Rcpp_0.12.19     cli_1.0.1        plyr_1.8.4       cellranger_1.1.0
[13] httr_1.3.1       tools_3.4.4      nlme_3.1-131.1   broom_0.5.0     
[17] R6_2.3.0         bindrcpp_0.2.2   bindr_0.1.1      scales_1.0.0    
[21] assertthat_0.2.0 gtable_0.2.0     stringi_1.1.7    rstudioapi_0.8  
[25] backports_1.1.2  hms_0.4.2        munsell_0.5.0    grid_3.4.4      
[29] colorspace_1.3-2 glue_1.3.0       lubridate_1.7.4  rlang_0.3.0.1   
[33] magrittr_1.5     lazyeval_0.2.1   yaml_2.2.0       crayon_1.3.4    
[37] haven_1.1.2      modelr_0.1.2     pillar_1.3.0     jsonlite_1.5

编辑2 使用Konrad的函数ks_cdf和

x = seq(0, 1, by = 0.01)
plot(x, ks_cdf(x))

当输入0时仍然返回0

编辑3 升级到版本3.6.1后

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)
...

我仍然得到与上述相同的情节，即ks_cdf(0)=0，而它应该是ks_sdf(0)=1。

- mjs

3

请澄清一下——为什么不能使用内置函数 ks.test？除此之外，错误很简单，只是一个拼写错误：termfb 的更新需要放在循环内部，而不是外部（这样是没有用的），您改变了返回值，并且在这样做时没有考虑到缺乏收敛性。 - Konrad Rudolph

1

顺便提一下，我知道内置函数。我需要它用于KS检验的2D版本，虽然有一个R包支持它，但没有显著性计算。 - mjs

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Konrad Rudolph · Accepted Answer

这段代码几乎可以字面意思翻译成R语言 - 不清楚为什么你在没有理由的情况下偏离了C语言代码。以下是一个字面上的、稍微整理过的翻译：

ks_cdf = function (lambda) {
  EPS1 = 0.001
  EPS2 = 1.0e-8
  sum = 0
  fac = 2
  termbf = 0
  a2 = -2 * lambda ^ 2

  for (j in 1 : 100) {
    term = fac * exp(a2 * j ^ 2)
    sum = sum + term
    if ((abs(term) <= EPS1 * termbf) || (abs(term) <= EPS2 * sum)) {
      return(sum)
    } else {
      fac = -fac
      termbf = abs(term)
    }
  }
  1 # Failed to converge.
}

这段代码是可以工作的，但它不是矢量化的。对于真实的实现，我会对此进行改变（但是通过这样做，我们将失去早期退出的功能）。

下面是使用矢量化算术和矩阵乘法的惯用R实现：

ks_cdf = function (λ) {
  eps1 = 0.001
  eps2 = 1E-8

  range = seq(1, 100)
  terms = (-1) ^ (range - 1) * exp(-2 * range ^ 2 %*% t(λ ^ 2))
  sums = 2 * colSums(terms)
  pterms = abs(terms)
  prev_pterms = rbind(0, pterms[-nrow(pterms), , drop = FALSE])
  converged = apply(pterms <= eps1 * prev_pterms | pterms <= eps2 * sums, 2L, any)
  sums[! converged] = 1
  sums
}

为了展示它的向量化效果以及这对于 IT 技术来说是件大事：

x = seq(0, 1, by = 0.01)
plot(x, ks_cdf(x))