## CARICAMENTO DATI A<-read.table(file="dati_Campi_Flegrei.txt",header=TRUE) Piro <- A[1:19,1] B<-read.table(file="test_medicina.txt",header=TRUE) Medi <- B[,2] Piro --------------------- ## ISTOGRAMMI A hist(Piro) hist(Piro,15) hist(Piro,15,freq=FALSE) hist(Medi,50,freq=FALSE) Piro.cut <- c(Piro[1:11],Piro[13:19]) hist(Piro.cut,7,freq=FALSE) ----------------------- ## CUMULATIVA plot.ecdf(Piro) plot.ecdf(Medi) ---------------------- ## LOG-NORMALE x<-1:100 y<- dlnorm(x,3,1) plot(x,y) ----------------------- ## KERNEL SMMOTHING require(KernSmooth) density <- bkde(Piro, kernel = "normal", bandwidth=20) plot(density, type="l") density.Medi <- bkde(Medi, kernel = "normal", bandwidth=20) plot(density.Medi, type="l") hist(Piro,15,freq=FALSE) lines(density, type="l") ---------------------- ## STIMA PARAMETRI require(MASS) fitdistr(Piro, "weibull") fitdistr(Piro, "weibull", list(shape = 0.5, scale = 20)) fitdistr(Piro, "weibull", list(shape = 2, scale = 100)) fitdistr(Piro, "log-normal") fitdistr(Medi, "normal") mean(Medi) sd(Medi) fitdistr(Medi, "weibull") ## esce errore Medi.sort <- sort(Medi) Medi.sort[1:10] [1] -7.00 -3.25 -1.25 1.25 1.50 2.50 3.75 4.75 5.75 [10] 6.00 Medi.plus <- Medi.sort[4:length(Medi)] fitdistr(Medi.plus, "weibull") ## risultati, per il seguito: fitdistr(Piro, "weibull") 0.85, 38.11 fitdistr(Piro, "log-normal") 3.09, 1.02 fitdistr(Medi, "normal") 34.97, 11.06 fitdistr(Medi, "weibull") 3.58, 38.84 ------------------------- ## CONFRONTO GRAFICO DENSITA'-ISTOGRAMMA a<-0.85 s<-38.11 x<-(-0:5000)/10 hist(Piro,15,freq=FALSE) y<-dweibull(x,a,s) lines(x,y) --- m<-3.09 s<-1.02 x<-(-0:5000)/10 hist(Piro,15,freq=FALSE) y<-dlnorm(x,m,s) lines(x,y) --- a<-0.8 s<-100 x<-(-0:5000)/10 hist(Piro,15,freq=FALSE) y<-dweibull(x,a,s) lines(x,y) --- --- --- m<-34.97 s<-11.06 x<-(-0:700)/10 hist(Medi,50,freq=FALSE) y<-dnorm(x,m,s) lines(x,y) --- a<-3.58 s<-38.84 x<-(-0:700)/10 hist(Medi.plus,30,freq=FALSE) y<-dweibull(x,a,s) lines(x,y) -------------------------------------- ## CONFRONTO GRAFICO CUMULATIVE a<-0.85 s<-38.11 x<-(-0:5000)/10 plot.ecdf(Piro) y<-pweibull(x,a,s) lines(x,y) --- m<-3.09 s<-1.02 x<-(-0:5000)/10 plot.ecdf(Piro) y<-plnorm(x,m,s) lines(x,y) --- a<-0.8 s<-100 x<-(-0:5000)/10 plot.ecdf(Piro) y<-pweibull(x,a,s) lines(x,y) --- --- --- m<-34.97 s<-11.06 x<-(-0:700)/10 plot.ecdf(Medi) y<-pnorm(x,m,s) lines(x,y, col="red") --- a<-3.58 s<-38.84 x<-(-0:700)/10 plot.ecdf(Medi) y<-pweibull(x,a,s) lines(x,y, col="red") -------------------------------------- ## CONFRONTO CAMPIONI CASUALI a<-0.85 s<-38.11 rweibull(19,a,s) Piro m<-3.09 s<-1.02 rlnorm(19,m,s) Piro m<-3.09 s<-1.3 rlnorm(19,m,s) m<-3.09 s<-1.3 x<-(-0:5000)/10 plot.ecdf(Piro) y<-plnorm(x,m,s) lines(x,y) -------------------------------------- ## Q-Q PLOT Dati <- Piro a<-0.85 s<-38.11 quant <- function(x) {qweibull(x,a,s)} x<- 1:500 L <- length(Dati) F.hat<- (1:L)/L - 0.5/L Dati.ord<-sort(Dati) plot(x,x, type="l") q <- quant(F.hat) lines(Dati.ord,q, type="b") ---- Dati <- Piro m<-3.09 s<-1.02 quant <- function(x) {qlnorm(x,m,s)} x<- 1:500 ---- Dati <- Piro m<-3.09 s<-1.3 quant <- function(x) {qlnorm(x,m,s)} x<- 1:500 --- --- --- Dati <- Medi m<-34.97 s<-11.06 quant <- function(x) {qnorm(x,m,s)} x<-(-0:700)/10 --- Dati <- Medi a<-3.58 s<-38.84 quant <- function(x) {qweibull(x,a,s)} x<-(-0:700)/10 ----------------------------------- ## UNA DISTANZA TRA CUMULATIVE Dati <- Piro a<-0.85 s<-38.11 F <- function(x) {pweibull(x,a,s)} x<- 1:500 plot.ecdf(Dati) y<-F(x) lines(x,y) L <- length(Dati) F.hat<- (1:L)/L - 0.5/L Dati.ord<-sort(Dati) Int<- 1:(L+1) Int[1]<-integrate(F, -Inf, Dati.ord[1])$value F.end <- function(x) {1-F(x)} Int[L+1]<-integrate(F.end, Dati.ord[L], Inf)$value for (k in 1:(L-1)) { F.k <- function(x) {abs(F(x)-(k-0.5)/L)} Int[k+1]<-integrate(F.k, Dati.ord[k], Dati.ord[k+1])$value } I <- sum(Int) Range <- max(Dati)-min(Dati) E1 <- 100*I/Range E1 --- Dati <- Piro a<-0.8 s<-100 F <- function(x) {pweibull(x,a,s)} x<- 1:500 --- Dati <- Piro a<-0.6 s<-25 F <- function(x) {pweibull(x,a,s)} x<- 1:500 --- --- Dati <- Piro m<-3.09 s<-1.02 F <- function(x) {plnorm(x,m,s)} x<- 1:500 ---- Dati <- Piro m<-3.09 s<-1.3 F <- function(x) {plnorm(x,m,s)} x<- 1:500 ------------------------------------------------ ## TEST BASATO SULLA PRECEDENTE DISTANZA Dati <- Piro m<-3.09 s<-1.02 F <- function(x) {plnorm(x,m,s)} Ran <- function(x) {rlnorm(x,m,s)} N<-1000 L <- length(Dati) F.hat<- (1:L)/L - 0.5/L Dati.ord<-sort(Dati) Int<- 1:(L+1) Int[1]<-integrate(F, -Inf, Dati.ord[1])$value F.end <- function(x) {1-F(x)} Int[L+1]<-integrate(F.end, Dati.ord[L], Inf)$value for (k in 1:(L-1)) { F.k <- function(x) {abs(F(x)-(k-0.5)/L)} Int[k+1]<-integrate(F.k, Dati.ord[k], Dati.ord[k+1])$value } I <- sum(Int) Range <- max(Dati)-min(Dati) E1 <- 100*I/Range I.rand <- 1:N Int.rand<- 1:(L+1) for (n in 1:N) { Dati.rand<-Ran(L) D.r.o<-sort(Dati.rand) Int.rand[1]<-integrate(F, -Inf, D.r.o[1])$value Int.rand[L+1]<-integrate(F.end, D.r.o[L], Inf)$value for (k in 1:(L-1)) { F.k <- function(x) {abs(F(x)-(k-0.5)/L)} Int.rand[k+1]<-integrate(F.k, D.r.o[k], D.r.o[k+1])$value } I.rand[n] <- sum(Int.rand) } hist(100*I.rand/Range) I.rand.ord <- sort(I.rand) p.value <- 1-which.min(abs(I.rand.ord-I))/N E1 p.value --- Dati <- Piro a<-0.85 s<-38.11 F <- function(x) {pweibull(x,a,s)} Ran <- function(x) {rweibull(x,a,s)} N<-1000 --------------------------------------------------- --------------------------------------------------- --------------------------------------------------- ## MOMENTI DI ORDINE SUPERIORE ## (NON LEGGERE SENZA PROVARE PER ESERCIZIO) Dati <- Piro m<-3.09 s<-1.02 F <- function(x) {plnorm(x,m,s)} Ran <- function(x) {rlnorm(x,m,s)} N<-1000 K<-2 L <- length(Dati) F.hat<- (1:L)/L - 0.5/L Dati.ord<-sort(Dati) Int<- 1:(L+1) Int.K<- 1:(L+1) F.K <- function(x) {F(x)^K} Int[1]<-integrate(F, -Inf, Dati.ord[1])$value Int.K[1]<-integrate(F.K, -Inf, Dati.ord[1])$value F.end <- function(x) {1-F(x)} F.end.K <- function(x) {(1-F(x))^K} Int[L+1]<-integrate(F.end, Dati.ord[L], Inf)$value Int.K[L+1]<-integrate(F.end.K, Dati.ord[L], Inf)$value for (k in 1:(L-1)) { F.k <- function(x) {abs(F(x)-(k-0.5)/L)} F.K.k <- function(x) {abs(F(x)-(k-0.5)/L)^K} Int[k+1]<-integrate(F.k, Dati.ord[k], Dati.ord[k+1])$value Int.K[k+1]<-integrate(F.K.k, Dati.ord[k], Dati.ord[k+1])$value } I <- sum(Int) I.K <- sum(Int.K) Range <- max(Dati)-min(Dati) E1 <- 100*I/Range Ek <- 100*(I.K/Range)^(1/K) I.rand <- 1:N I.rand.K <- 1:N Int.rand<- 1:(L+1) Int.rand.K<- 1:(L+1) for (n in 1:N) { Dati.rand<-Ran(L) D.r.o<-sort(Dati.rand) Int.rand[1]<-integrate(F, -Inf, D.r.o[1])$value Int.rand.K[1]<-integrate(F.K, -Inf, D.r.o[1])$value Int.rand[L+1]<-integrate(F.end, D.r.o[L], Inf)$value Int.rand.K[L+1]<-integrate(F.end.K, D.r.o[L], Inf)$value for (k in 1:(L-1)) { F.k <- function(x) {abs(F(x)-(k-0.5)/L)} F.K.k <- function(x) {abs(F(x)-(k-0.5)/L)^K} Int.rand[k+1]<-integrate(F.k, D.r.o[k], D.r.o[k+1])$value Int.rand.K[k+1]<-integrate(F.K.k, D.r.o[k], D.r.o[k+1])$value } I.rand[n] <- sum(Int.rand) I.rand.K[n] <- sum(Int.rand.K) } hist(100*(I.rand.K/Range)^(1/K)) I.rand.ord <- sort(I.rand) I.rand.ord.K <- sort(I.rand.K) p.value <- 1-which.min(abs(I.rand.ord-I))/N p.v.k <- 1-which.min(abs((I.rand.ord.K/Range)^(1/K)-(I.K/Range)^(1/K)))/N c(E1,Ek) c(p.value,p.v.k) ----------------------------------------------