GaryBAYLOR / R-code

A collection of algorithms written by myself for solving statistical problems
0 stars 0 forks source link

Fligner-Policello test in R #10

Open GaryBAYLOR opened 9 years ago

GaryBAYLOR commented 9 years ago

Fligner-Policello Test in R

I try to find a package to do FP test, but found nothing. So I make a small crude function which take two vectors as input and return a p-value.

FP.test <- function(x, y) {
    n <- length(x); m <- length(y)
    Pl.x <- numeric(n); Pl.y <- numeric(m)

    for(i in 1:n) {
        Pl.x[i] <- sum(y < x[i]) + 1/2 * sum(y == x[i])
    }

    for(i in 1:m) {
        Pl.y[i] <- sum(x < y[i]) + 1/2 * sum(x == y[i])
    }

    Plxmean <- mean(Pl.x); Plymean <- mean(Pl.y)

    Vx <- (n - 1) * var(Pl.x); Vy <- (m - 1) * var(Pl.y)

    z <- (sum(Pl.y) - sum(Pl.x)) / (2 * sqrt(Vx + Vy + Plxmean * Plymean))

    p <- pnorm(abs(z), lower.tail = FALSE) * 2

    return (list(p.value = p))

}

We can try some example

> set.seed(100)
> x <- rnorm(20, mean = 0, sd = 1)
> set.seed(200)
> y <- rnorm(30, mean = 0, sd = 1)
> set.seed(201)
> y1 <- rnorm(30, mean = 1, sd = 1)
> set.seed(202)
> y2 <- rnorm(30, mean = 2, sd = 1)
> set.seed(203)
> y3 <- rnorm(30, mean = 3, sd = 1)
> FP.test(x,y)
$p.value
[1] 0.7070176

> FP.test(x,y1)
$p.value
[1] 5.794937e-05

> FP.test(x,y2)
$p.value
[1] 2.43404e-46

> FP.test(x,y3)
$p.value
[1] 1.124729e-250