Merge branch 'development' (version 1.3.1)

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: brunnermunzel
 Type: Package
 Title: (Permuted) Brunner-Munzel Test
-Version: 1.3.0
-Date: 2019-03-10
+Version: 1.3.1
+Date: 2019-03-11
 Authors@R: person("Toshiaki", "Ara", email = "toshiaki.ara@gmail.com",
   role = c("aut", "cre"))
 License: GPL-2 | GPL-3

NEWS.md

@@ -1,5 +1,9 @@
 # ChangeLog
-## version 1.3.0
+## version 1.3.1 (2019-03-11)
+* rewrote completely with FORTRAN (PR#12)
+  * calculation of P value
+
+## version 1.3.0 (2019-03-10)
 * bugfix (PR#9)
   * return appropriate values using non-overlapping data
     * estimate, confidence interval and P value

R/brunnermunzel.perutation.test.R

@@ -161,25 +161,24 @@ brunnermunzel.permutation.test.default <-
             return(res)
         }
 
-        res <- .Fortran("bm_permutation_stat",
+        alter <- switch(alternative,
+                        "two.sided" = 1L,
+                        "greater" = 2L,
+                        "less" = 3L)
+
+        res <- .Fortran("bm_permutation_test",
                         n = as.integer(nx + ny),
                         r = as.integer(nx),
                         n_nCr = as.integer(n_nCr),
                         dat = as.double(c(x, y)),
-                        statistics = numeric(n_nCr),
+                        alter = as.integer(alter),
+                        pval = numeric(1),
                         PACKAGE = "brunnermunzel")
-        z1 <- res$statistics
-        z0 <- z1[1L]
-
-        p.value <- switch(alternative,
-                          "two.sided" = mean(abs(z1) >= abs(z0)),
-                          "greater" = mean(z1 <= z0),
-                          "less" = mean(z1 >= z0))
 
         structure(
             list(
                 method = "permuted Brunner-Munzel Test",
-                p.value = p.value,
+                p.value = res$pval,
                 data.name = DNAME),
             class = "htest")
     }

src/bm_permutation_test.f

@@ -0,0 +1,51 @@
+*
+*     bm_permutation_test:
+*     perform permuted Brunner-Munzel test
+*
+*     (input)
+*     n    : length of data
+*     r    : length of provided data
+*     n_nCr: choose(n, r) in R function
+*     dat  : provided data (length is n)
+*     alter : alterative (1: "two.sided", 2: "greater", 3: "less")
+*     (output)
+*     pval : P value
+*
+      subroutine bm_permutation_test(n, r, n_nCr, dat, alter, pval)
+      implicit none
+      integer,intent(in)::n,r,n_nCr,alter
+      double precision,intent(in)::dat(n)
+      double precision,intent(out)::pval
+
+      integer::i,count
+      double precision statistics(n_nCr)
+
+      call bm_permutation_stat(n, r, n_nCr, dat, statistics)
+
+!     description by R
+!     - z0: statistic by observed data, statistic(1)
+!     - z1: statistic by permuted data, statistic(i)
+!     - "two.sided": mean(abs(z1) >= abs(z0))
+!     - "greater"  : mean(z1 <= z0)  i.e. mean(-z1 >= -z0)
+!     - "less"     : mean(z1 >= z0)
+      if (alter.eq.1) then      ! "two sided"
+         statistics(1:n_nCr) = abs(statistics(1:n_nCr))
+      else if (alter.eq.2) then ! "greater"
+         statistics(1:n_nCr) = -statistics(1:n_nCr)
+      endif
+
+      count = 0
+      do i = 1, n_nCr
+         if (statistics(i).ge.statistics(1)) then
+            count = count + 1
+         endif
+      enddo
+
+      pval = dble(count) / n_nCr
+      return
+      end
+
+
+
+
+

src/calc_statistics.f

@@ -51,11 +51,6 @@ subroutine calc_statistics(nx, ny, dat, const, idx, stat)
       enddo
 
       v = const(3) * vx + const(4) * vy
-      ! to avoid division by zero
-      if (v.lt.0.000001) then
-         v = 0.00001
-      endif
-
       stat = (my - mx) / sqrt(v)
       return
       end