Discussione:Distribuzione di Kolmogorov-Smirnov

Script R per calcolare il valori delle distribuzioni

Al fine di non dover copiare le tavole dei valori critici da qualche libro o qualche software, si riporta lo script in R con il quale è possibile calcolare le distribuzioni esatte di D_n,m e di conseguenza le tavole dei valori critici.

In questa versione calcola solo la distribuzione da usare nel test di Kolmogorov-Smirnov a due code.

##################################################################
library(combinat)

ks <- function(n=2,m=2) {
    c <-nCm(n+m,m)
    comb <- t(combn(n+m, m, tabulate, nbins = n+m)) 

    D <- rep(0,c)
    Dnm <- D

    for (i in 1:c) {
        a <- (1:(n+m))[comb[i,]==1]
        b <- (1:(n+m))[comb[i,]==0]
        ks <- ks.test(a,b)
        D[i] <- ks$statistic
        Dnm[i] <- round(D[i]*n*m)
    }
    tableDnm <- table(Dnm)
    fD <- data.frame(tableDnm)
    fD$Dnm <- sort(unique(Dnm))
    k <- length(tableDnm)
    FD <- fD
    alfaD <- fD
    for (i in 2:k) {FD$Freq[i]<- FD$Freq[i-1]+fD$Freq[i]}
    for (i in (k-1):1) {alfaD$Freq[i]<- alfaD$Freq[i+1]+alfaD$Freq[i]}
    cbind(n=rep(n,k),m=rep(m,k),D=fD$Dnm/(n*m),fD,nCm=c,f=fD$Freq/c,F=FD$Freq/c,alfa=alfaD$Freq/c)
}

# valori limite per rifiutare l'ipotesi nulla con una confidenza pari a 1-alfa
N <- 10
Dalfa_01 <- matrix(0,N,N)
Dalfa_05 <- matrix(0,N,N)
Dalfa_10 <- matrix(0,N,N)
print(c('n','m','0.90','0.95','0.99'))
for (n in 2:N) {for (m in 2:n) {
    butta<-ks(n,m); 
    print(c(n,m, n*m*c(min((butta$D)[butta$alfa<0.1]),
                       min((butta$D)[butta$alfa<0.05]),
                       min((butta$D)[butta$alfa<0.01])
                       )
         ))
    Dalfa_01[n,m] <- n*m*min((butta$D)[butta$alfa<0.01]); Dalfa_01[m,n] <- Dalfa_01[n,m];
    Dalfa_05[n,m] <- n*m*min((butta$D)[butta$alfa<0.05]); Dalfa_05[m,n] <- Dalfa_05[n,m];
    Dalfa_10[n,m] <- n*m*min((butta$D)[butta$alfa<0.10]); Dalfa_10[m,n] <- Dalfa_10[n,m];
    } 
} 

# stampa a video delle tabelle con i valori critici
print(Dalfa_01)
print(Dalfa_05)
print(Dalfa_10)
##################################################################