# 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)