#--------------------------------------------------------------
# Propagazione delle incertezze sui  dati del foglio A4
# - valore atteso e varianza di combinazioni lineari 
#   (eventualmente previa linearizzazione)
# - propagazione mediante Monte Carlo
#
# - Variante rispetto a foglio_A4.R: 
#   la propagazione tramite Monte Carlo usa una distr. uniforme   
#
#  NOTA: Versione senza i conti, a parte per il perimetro
#        -> Problema mr 38.1
# 
#  GdA, giugn0 2023
#--------------------------------------------------------------

pausa <- function() { cat (">> press Enter to continue\n"); scan() }

# sovrascriviamo rnorm() con una uniforme (DA RIMUOVERE alla fine!)
rnorm <- function(N, mu=0, sigma=1) {
    mu + runif(N, -1/2, 1/2) * sigma / (1/sqrt(12))
}

Ea <- 29.73  # cm
sa <- 0.03   # cm
Eb <- 21.45  # cm
sb <- 0.04   # cm

N <- 10000
# sampling
a <- rnorm(N, Ea, sa)
b <- rnorm(N, Eb, sb)
cat(sprintf(" a = %.2f +- %.2f cm  (MC:  %.2f +- %.2f cm)) \n",
            Ea, sa, mean(a), sd(a) ))
cat(sprintf(" b = %.2f +- %.2f cm  (MC:  %.2f +- %.2f cm))\n",
            Eb, sb, mean(b), sd(b) ))


# perimetro
Ep <- 2*Ea + 2*Eb
sp <- sqrt(4*sa^2 + 4*sb^2)
p <- 2*a + 2*b
cat(sprintf(" p = %.2f +- %.2f cm  (MC:  %.2f +- %.2f cm)\n",
            Ep, sp, mean(p), sd(p)))

# Area
EA <- NA # lasciato come esercizio
sA <- NA # lasciato come esercizio
A <- a * b
cat(sprintf(" A = %.1f +- %.1f cm^2   (MC:  %.1f +- %.1f cm^2)\n",
            EA, sA, mean(A), sd(A)))

# diagonale
Ed <- NA # lasciato come esercizio
sd <- NA # lasciato come esercizio
d <- sqrt(a^2+b^2)
cat(sprintf(" d = %.2f +- %.2f cm   (MC:  %.2f +- %.2f cm)\n",
            Ed, sd, mean(d), sd(d)))
# nota: per R 'sd' e 'sd()' sono due cose diverse!

# rapporto
Er <-  NA # lasciato come esercizio
sr <-  NA # lasciato come esercizio
r <- a/b
cat(sprintf(" r = %.3f +- %.3f   (MC:  %.3f +- %.3f)\n",
            Er, sr, mean(r), sd(r)))

pausa()
plot(a, b, cex=0.5, col='blue')

pausa()
# Correlazioni (alcuni esempi)
plot(a,p, pch=19, cex=0.5, col='blue')
cat(sprintf(" rho(p, a) = %.3f \n", cor(p,a) )) 
abline(v=Ea, lty=2, col='blue')
abline(h=Ep, lty=2, col='blue')

pausa()
plot(a,A, pch=19, cex=0.5, col='blue')
cat(sprintf(" rho(A, a) = %.3f \n", cor(A,a) )) 
abline(v=Ea, lty=2, col='blue')
abline(h=EA, lty=2, col='blue')

pausa()
plot(a,r, pch=19, cex=0.5, col='blue')
cat(sprintf(" rho(r, a) = %.3f \n", cor(r,a) )) 
abline(v=Ea, lty=2, col='blue')
abline(h=Er, lty=2, col='blue')

pausa()
plot(b,r, pch=19, cex=0.5, col='blue')
cat(sprintf(" rho(r, b) = %.3f \n", cor(r,b) )) 
abline(v=Eb, lty=2, col='blue')
abline(h=Er, lty=2, col='blue')

pausa()
plot(r,A, pch=19, cex=0.5, col='blue')
cat(sprintf(" rho(A, r) = %.3f \n", cor(r,A) )) 
abline(v=Er, lty=2, col='blue')
abline(h=EA, lty=2, col='blue')

# RIMUOVIAMO rnorm !
rm(rnorm)