# distribuzione triangolare - GdA, gennaio 05 # attivare con source("triang.R"), quindi eseguire per es. le prove # commentate in fondo ############################################################# dtriang <- function(x, x0=0, dm=1, dp=1) { # per evitare errori prendiamo i moduli di dm e dp dp <- abs(dp); dm <- abs(dm) h <- 2/(dp+dm) xmin <- x0 - dm; xmax <- x0 + dp m1 <- h/dm m2 <- -h/dp # prima valutazione incurante dei limiti (trucco per # poter agire su un vettore) f <- m1 * (x-x0) + h f [x > x0 ] <- m2 * (x[ x > x0]-x0) + h #f [x < xmin | x > max ] <- 0 f [x < xmin ] <- 0 f [x > xmax ] <- 0 return(f) } ptriang <- function(x, x0=0, dm=1, dp=1) { dp <- abs(dp); dm <- abs(dm) h <- 2/(dp+dm) pm <- dm * h /2 xmin <- x0 - dm; xmin2 <- xmin^2 xmax <- x0 + dp; xmax2 <- xmax^2 m1 <- h/dm; m1half <- m1/2; hmm1x0 <- h - m1 * x0 ; k1 <- - hmm1x0 * xmin - m1half * xmin2 m2 <- -h/dp; m2half <- m2/2; hmm2x0 <- h - m2 * x0 ; k2 <- - hmm2x0 * x0 - m2half * x0 f <- k1 + hmm1x0 * x + m1half * x^2 f [x > x0 ] <- k2 + hmm2x0 * x[ x > x0] + m2half * x[ x > x0]^2 + pm f [x < xmin] <- 0 f [x > xmax] <- 1 return(f) } rtriang <- function(n, x0=0, dm=1, dp=1) { dp <- abs(dp); dm <- abs(dm) h <- 2/(dp+dm) pm <- dm * h /2 xmin <- x0 - dm; xmax <- x0 + dp sw <- runif(n) r <- sqrt(runif(n)) * dm + xmin r[sw > pm] <- xmax - sqrt( runif(length(r[sw > pm]) ) ) * dp return(r) } # manca la 'qtriang': chi si vuole cimentare nell'esercizio? #prove (i comandi che seguono sono commentati per permettere di eseguire # lo script per definire dtraing, ptriang e rtriang senza # che esso esegua anche questi ultimi) # # curve(ptriang(x,0,2,1),-2.5, 1.5, ylab='f(x) e F(x)', main='Triangolare', col='green') # curve(dtriang(x,0,2,1),-2.5, 1.5, col='blue', add=T) # # r <- rtriang(1000000,1,1,2) # hist(r,ncl=100,prob=T) # curve(dtriang(x,1,1,2),add=T) # # r1 <- rtriang(1000000, 0.5, 1.5, 0.5) # mean(r1) # sd(r1) # summary(r1) # r2 <- rtriang(1000000, 0.5, 1.5, 0.5) # mean(r2) # sd(r2) # summary(r2) # s <- r1+r2 # mean(s) # sd(s) # summary(s) # # r3 <- rtriang(1000000, 0.5, 1.5, 0.5) # r4 <- rtriang(1000000, 0.5, 1.5, 0.5) # r5 <- rtriang(1000000, 0.5, 1.5, 0.5) # r6 <- rtriang(1000000, 0.5, 1.5, 0.5) # r7 <- rtriang(1000000, 0.5, 1.5, 0.5) # r8 <- rtriang(1000000, 0.5, 1.5, 0.5) # r9 <- rtriang(1000000, 0.5, 1.5, 0.5) # r10 <- rtriang(1000000, 0.5, 1.5, 0.5) # s10 <- r1+r2+r3+r4+r5+r6+r7+r8+r9+r10 # layout(matrix(1:3, 3, 1)) # hist(r1, ncl=25,xlim=c(-3,6), main="r_i triangolare", prob=T) # hist(s, ncl=40,xlim=c(-3,6), main="Somma 2 triangolari", xlab="s2 = r1 + r2") # hist(s10, ncl=100,xlim=c(-3,6), main="Somma 10 triangolari", xlab="s2 = r1 + ..+ r10",prob=T) # curve(dnorm(x,mean(s10),sd(s10)),add=T,col='red') # layout(1)