Skip to content

Commit b1a8872

Browse files
author
hornik
committed
Add dpq functions for the exact conditional/permutation distribution of the Wilcoxon signed rank statistic.
git-svn-id: https://svn.r-project.org/R/trunk@88752 00db46b3-68df-0310-9c12-caf00c1e9a41
1 parent 30dfb1f commit b1a8872

File tree

1 file changed

+66
-0
lines changed

1 file changed

+66
-0
lines changed

src/library/stats/R/wilcox.test.R

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -684,3 +684,69 @@ function(p, m, n, z = NULL, lower.tail = TRUE)
684684
y[i] <- vapply(p[i], function(e) s[v >= e][1L], 0)
685685
y
686686
}
687+
688+
.dsignrank <-
689+
function(x, n, z = NULL)
690+
{
691+
if(is.null(z))
692+
return(dsignrank(x, n))
693+
694+
stopifnot(length(z) == n)
695+
if (!all(2 * z == floor(2 * z)) || any(z < 1))
696+
stop("'z' is not a rank vector")
697+
y <- rep.int(NA_real_, length(x))
698+
i <- which(!is.na(x))
699+
if (!any(i))
700+
return(y)
701+
f <- 2 - (max(z - floor(z)) == 0)
702+
d <- .Call(C_cpermdist1,
703+
as.integer(sort(floor(f * z))))
704+
w <- seq.int(0, length(d) - 1L)
705+
x <- f * x[i]
706+
w <- w[match(x, w)] + 1L
707+
y[i] <- ifelse(is.na(w), 0, d[w])
708+
y
709+
}
710+
711+
.psignrank <-
712+
function(q, n, z = NULL, lower.tail = TRUE)
713+
{
714+
if(is.null(z))
715+
return(psignrank(q, n, lower.tail = lower.tail))
716+
717+
y <- rep.int(NA_real_, length(q))
718+
i <- which(!is.na(q))
719+
if(!any(i))
720+
return(y)
721+
722+
## Support of V
723+
s <- seq.int(0, n * (n + 1)) / 2
724+
## Density
725+
d <- .dsignrank(s, n, z)
726+
y[i] <- vapply(q[i],
727+
function(e)
728+
sum(d[s < e + sqrt(.Machine$double.eps)]),
729+
0)
730+
if(lower.tail) y else 1 - y
731+
}
732+
733+
.qsignrank <-
734+
function(p, n, z = NULL, lower.tail = TRUE)
735+
{
736+
if(is.null(z))
737+
return(qsignrank(p, n, lower.tail = lower.tail))
738+
739+
if (any(i <- (p < 0) | (p > 1)))
740+
y[i] <- NaN
741+
i <- !is.na(p) & !i
742+
if (!any(i))
743+
return(y)
744+
745+
s <- seq.int(0, n * (n + 1)) / 2
746+
v <- .psignrank(s, n, z)
747+
if (!lower.tail)
748+
p <- 1 - p
749+
p <- p - 10 * .Machine$double.eps
750+
y[i] <- vapply(p[i], function(e) s[v >= e][1L], 0)
751+
y
752+
}

0 commit comments

Comments
 (0)