## EXAMPLE OF TEST (FOR THE AVERAGE)

delay <- c(5,7,4,10,6,5,8,2,8,6)
mu.declared <- 5
average.delay <- mean(delay)
average.delay
sd.delay <- sd(delay)
sd.delay
x.bar <- 1:10000
for (i in 1:10000) {
x.bar[i] <- mean(rnorm(10,5,sd.delay))
}
hist(x.bar,100)


1-which.min(abs(sort(x.bar)-average.delay))/10000


## DENSITY FIT TESTS FOR PIRO DATA 
----------------------------
## PREPARATION

A<-read.table(file="dati_Campi_Flegrei.txt",header=TRUE)
Piro <- A[1:19,1]

require(MASS)

Dati.18 <- sort(Piro)[1:18]
fitdistr(Dati.18, "weibull")

     shape        scale   
   1.4404892   28.2534758 
 ( 0.2669708) ( 4.8848966)

Dati.19 <- Piro
fitdistr(Dati.19, "weibull")

     shape        scale   
   0.8538612   38.0241762 
 ( 0.1325941) (10.8614965)



fitdistr(Dati.19, "log-normal")

   meanlog      sdlog  
  3.0994303   1.0235558 
 (0.2348198) (0.1660427)


---------------------------------------
## CHI SQUARE TEST WITHOUT OUTLIER, ML WEIBULL

## search for a reasonable range

min <- min(Dati.18)
max <- max(Dati.18)
c(min,max)

## split the range in n.class intervals

n.class <- 5
points <- (0:n.class)/n.class*80

## compute empirical frequences

hist <- hist(Dati.18, br = points, freq = FALSE)
hist
f <- hist$counts

## compute theoretical frequences

a <- 1.44
s <- 28.25

ecdf <- pweibull(points[2:(n.class+1)],a,s)
ecdf.shift <- pweibull(points[1:n.class],a,s)
p <- (ecdf-ecdf.shift)/sum((ecdf-ecdf.shift))
plot(p)

## finally, the test

chisq.test(f, p=p, rescale.p = TRUE)


---------------------------------------
## CHI SQUARE TEST WITHOUT OUTLIER, OTHER EXAMPLES

## the result depends on n.class

n.class <- 2
points <- (0:n.class)/n.class*80
hist <- hist(Dati.18, br = points, freq = FALSE)
hist
f <- hist$counts
a <- 1.44
s <- 28.25
ecdf <- pweibull(points[2:(n.class+1)],a,s)
ecdf.shift <- pweibull(points[1:n.class],a,s)
p <- (ecdf-ecdf.shift)/sum((ecdf-ecdf.shift))
chisq.test(f, p=p, rescale.p = TRUE)

n.class <- 10
points <- (0:n.class)/n.class*80
hist <- hist(Dati.18, br = points, freq = FALSE)
hist
f <- hist$counts
a <- 1.44
s <- 28.25
ecdf <- pweibull(points[2:(n.class+1)],a,s)
ecdf.shift <- pweibull(points[1:n.class],a,s)
p <- (ecdf-ecdf.shift)/sum((ecdf-ecdf.shift))
chisq.test(f, p=p, rescale.p = TRUE)

## the result is still good with little change of parameters

n.class <- 5
points <- (0:n.class)/n.class*80
hist <- hist(Dati.18, br = points, freq = FALSE)
hist
f <- hist$counts
a <- 1
s <- 30
ecdf <- pweibull(points[2:(n.class+1)],a,s)
ecdf.shift <- pweibull(points[1:n.class],a,s)
p <- (ecdf-ecdf.shift)/sum((ecdf-ecdf.shift))
chisq.test(f, p=p, rescale.p = TRUE)

## the result deteriorates with large change of parameters

n.class <- 5
points <- (0:n.class)/n.class*80
hist <- hist(Dati.18, br = points, freq = FALSE)
hist
f <- hist$counts
a <- 3
s <- 40
ecdf <- pweibull(points[2:(n.class+1)],a,s)
ecdf.shift <- pweibull(points[1:n.class],a,s)
p <- (ecdf-ecdf.shift)/sum((ecdf-ecdf.shift))
chisq.test(f, p=p, rescale.p = TRUE)




---------------------------------------
## KS TEST WITHOUT OUTLIER

a <- 1.44
s <- 28.25

ks.test(Dati.18, "pweibull", a, s)

Piro.mod <- Piro
Piro.mod[18]<- Piro.mod[18]+0.01
Piro.mod[16]<- Piro.mod[16]+0.01
Piro.mod[15]<- Piro.mod[15]+0.01
Dati.mod.18 <- sort(Piro.mod)[1:18]

ks.test(Dati.mod.18, "pweibull", a, s)

---------------------------------------
## KS TEST WITH OUTLIER, ML Weibull

a <- 0.85
s <- 38

ks.test(Piro.mod, "pweibull", a, s)

---------------------------------------
## KS TEST WITH OUTLIER, ML log-normal

m <- 3.1
s <- 1.02

ks.test(Piro.mod, "plnorm", m, s)

---------------------------------------
## KS TEST WITH OTHER HYPOTESIS

a <- 1
s <- 30

ks.test(Piro.mod, "pweibull", a, s)

a <- 3
s <- 40

ks.test(Piro.mod, "pweibull", a, s)

-----------------------------------------
-----------------------------------------