#
# Classe QUARTA    (uscite)
#
# (1)
#
# Approfondiremo l'uso di R, che già conoscete e potete scaricare
# sul computer collegandovi a  http://macosa.dima.unige.it/prog.htm
# Un rapido riepilogo.
# Si possono copiare le righe di comandi come le seguenti e poi
# incollarle (azionando il comando col mouse o con Crtl+V) in R.
# Nel programma i comandi appaiono preceduti da >.
# Il cancelletto precede un commento: la riga non viene eseguita dal
# computer
# Nel seguito le parti separate da # sono da eseguire una dopo l'altra:
# le parti tra un # e il successivo possono essere copiate e incollate
# in R (R le esegue automaticamente).
#
3+100*7^2
# Per separare parte intera e parte frazionaria si usa ".", non ",".
#
24; 24*60; 24*60*60; 24*60*60*100
#
2/3+1/5
#
# Come eseguire i calcoli tra frazioni in modo esatto?
library(MASS)
fractions(2/3+1/5)
#
# il simbolo "c" indica una collezione di oggetti, su cui si
# possono eseguire le stesse operazioni fatte sui singoli oggetti
#
c(24, 24*60, 24*60*60, 24*60*60*100)
3 + c(24, 24*60, 24*60*60, 24*60*60*100) * 100
#
# per indicare la notazione esponenziale viene usata "e"
2^2; 2^3; 2^20; 2^30; 2^40; 2^-2; 2^-20
#
# Schiacciando il tasto freccia ^ [v] si rivedono preced. [segu.] comandi:
# cosa COMODA per rivedere, correggere o modificare precedenti comandi.
#
# (2)
#
# Prima di fare cose più serie, un giochino. Copiate ed eseguite le righe
# seguenti. Poi le esamineremo rapidamente.
# /-----
# Preparo un "orologio" (proc.time()[3] dà il tempo in sec dall'inizio):
tic <- function(x) {sec <-proc.time()[3]; while(proc.time()[3] < sec+x) sec <- sec} # passano x sec
# Faccio un "ciclo", che parte da "for" e finisce con "}"
for (k in 0:30) { tic(0.3);
plot(c(0,100),c(0,100),type="n",xlab="", ylab="", asp=1);
symbols(50,80, circles=5, inches=FALSE, add=TRUE);
# occorre mettere inches=FALSE se no vengono usati i pollici
symbols(52,82, circles=1, inches=FALSE, add=TRUE);
lines(c(55,50,50),c(0,0,75)); tic(0.3);
plot(c(0,100),c(0,100),type="n",xlab="", ylab="", asp=1);
symbols(50,80, circles=5, inches=FALSE, add=TRUE);
symbols(52,82, circles=1, inches=FALSE, add=TRUE);
lines(c(35,30,50,70,75,70,50,50,35,50,65),c(1.5,2,40,2,2.5,2,40,75,45,75,45)) }
# Aggiungo una scritta:
lines(c(52,54.5),c(79,79)); text(71,80,"CIAO")
# \-----
# for esegue 31 volte le cose elencate tra { e }, che consitono nel
# tracciare, ogni 0.3 sec, un omino tracciato con dei segmenti (mediante
# lines) e due cerchi (mediante symbols); l'omino cambia via via
# posizione (gambe e braccia allineate al corpo, e poi allargate).
# Queste righe erano un pretesto, per introdurre alcuni comandi e passare
# a cose (forse) più serie.
#
# (3)
#
# Proponiamoci di "studiare" una funzione, ad es. questa:
q <- function(x) (2*x^2-1)/(x^2-2*x-5)
# Ne facciamo il grafico, tracciando gli assi
plot(q,-10,10); abline(h=0,v=0)
# Otteniamo un piccolo sgorbio: ci rendiamo conto che ci sono dei punti
# in cui la funzione non ha definita e ha dei salti, che il programma
# congiunge: il programma tabula la funzione per una quantità finita
# di valori e poi congiunge i punti ottenuti.  Facciamo un grafico "per
# punti piccoli" (lo facciamo con: type="p", pch=".").
plot(q,-10,10,ylim=c(-20,20),type="p", pch="."); abline(h=0,v=0)
# aumentiamo il numero dei punti, ad es. per 10000 = 10^5
plot(q,-10,10,ylim=c(-20,20),type="p",pch=".",n=10^5); abline(h=0,v=0)
# OK
# Vediamo dove si azzera il denominatore:
h <- function(x) x^2-2*x-5; plot(h,-3,5); abline(h=0,v=0)
# Il grafico attraversa l'asse delle x tra -2 e 0 e tra 2 e 4.
# Possiamo trovare le "radici" della funzione (ossia dove essa vale 0)
# col comando uniroot ("root" vuol dire "radice"), specificando la
# "tolleranza" (mettiamo un valore molto basso, in modo da avere
# tutte cifre buone):
uniroot(h,c(2,4),tol=10^-12)$root; uniroot(h,c(-2,0),tol=10^-12)$root
# Controlliamo graficamente se sono i valori in cui q non è definita:
plot(q,-10,10,ylim=c(-20,20),type="p",pch=".",n=10^5); abline(h=0,v=0)
abline(v=c(3.44949,-1.44949), col="blue")
# OK
# Vediamo come si comporta per x -> Inf (sappiamo che il limite è 2
# in quanto in (2*x^2-1)/(x^2-2*x-5) "vincono" sopra 2*x^2 e sotto x^2
# per cui tende a comportarsi come (2*x^2)/(x^2) = 2.
q(100)
# q(100) è quasi 2; vediamo meglio la cosa per 10^3,...,10^10:
q(10^c(3,4,5,6,7,8,9,10))
#
# Ora troviamo massimi e minimi. Potremmo farlo con vari metodi.
# Facciamolo calcolando la derivata di q. Anche in questo caso
# usiamo il programma:
D(expression((2*x^2-1)/(x^2-2*x-5)),"x")
# Ho specificato l'espressione della funzione e la variabile rispetto a
# cui voglio derivare. Ho trovato la derivata, che chiamo ad es. Dq
Dq <- function(x) 4*x/(x^2-2*x-5)-(2*x^2-1)*(2*x-2)/(x^2-2*x-5)^2
# Ne sovrappongo il grafico ai precedenti:
plot(Dq, -10,10, col="red", add=TRUE)
# Vediamo meglio la situazione con uno zoom (il grafico nero è della
# funzione, quello rosso della sua derivata):
plot(q,-15,10,ylim=c(-1/2,3),type="p",pch=".",n=1e5); abline(h=0,v=0)
abline(v=c(3.44949,-1.44949), col="blue")
plot(Dq,-15,10,col="red",add=TRUE)
# Vediamo meglio il grafico della derivata:
plot(Dq,-15,10,ylim=c(-1/2,1/2),col="red"); abline(h=0,v=0)
# Troviamo dove si azzera la derivata Dq, saranno i punti di massimo
# o minimo di q. Usiamo il solito comando uniroot:
uniroot(Dq, c(-10,-3),tol=10^-12)$root;  uniroot(Dq, c(-1/2,1/2),tol=10^-12)$root
# Verifichiamo la cosa tracciando le rette y=-4.386001 e y=-0.1139991:
abline(v=c(-4.386001,-0.1139991), col="blue")
# Vediamo meglio con uno zoom:
plot(Dq,-5.2,0.2,ylim=c(-0.2,0.8),col="red"); abline(h=0,v=0)
abline(v=c(-4.386001,-0.1139991), col="blue")
# OK
# Rivediamo, graficamente, quanto ottenuto:
plot(q,-6,6,ylim=c(-8,10),type="p",pch=".",n=10^5)
abline(h=0,v=0,col="red")
abline(v=c(3.44949,-1.44949), col="green")
abline(v=c(-4.386001,-0.1139991), col="blue")
#
# (4)
#
# Una applicazione alla FISICA:     la forza per ottenere
# l'allungamento di una molla (senza estenderla troppo).
F <- function(allungam) 1/20*allungam   # F in Newton, allungam in cm
plot(F,0,30); abline(h=0,v=0,col="blue")
text(16,1.4,"forza in funzione dell'allung.")
text(27,0.1,"cm"); text(1,1.45,"N")
#
# Vediamo questo esempio semplice, per illustrare tecniche che
# possono essere impiegate in ogni caso.
# Quanto lavoro occorre per allungare la molla di 15 cm?
for (x in 0:15) lines(c(x,x),c(0,F(x)),col="red")
# E' l'area di questo triangolo, che posso trovare facilmente:
F(15)*15/2
# In generale possiamo trovarla con la seguente funzione "area"
# (in matematica in genere chiamata "integrale"):
L <- function(allungam) area(F,0,allungam)
# L'area che sta tra grafico, asse x e le rette x=0 e x=allungam
# (per usarla dovrei caricare: library(MASS) ma l'ho già fatto)
L(15)
# Ottengo lo stesso valore ottenuto prima in altro modo (cosa che
# sapevo fare in questo caso particolare perché è un triangolo)
# Verifichiamo graficamente la cosa:
plot(L,0,30); abline(h=0,v=0,col="blue")
abline(v=15,h=5.625,col="red")
# OK
# In questo caso sapevamo trovare facilmente una funzione
# che abbia come derivata F; basta prendere:
L1 <- function(x) 1/20*x^2/2
# Verifica:
D(expression(1/20*x^2/2),"x")
# 1/20 * (2 * x)/2  equivale a  1/20*x
# Verifica "grafica":
plot(L1,0,30,col="green",add=TRUE)
# OK: il grafico di L1 si sovrappone a quello di L.
#
# (5)
#
# Ora facciamo qualcosina di probabilità e statistica.
#
# Avete già visto la funzione runif:
runif(1); runif(3); runif(6)
# runif(n) genera n numeri a caso tra 0 ed 1 con distribuzione uniforme.
# Infatti, generando 500*500 punti:
x <- runif(500); y <- runif(500); plot (x,y)
# Ma con questi comandi cosa ottengo? Perché?
x <- runif(500)*runif(500); y <- runif(500); plot (x,y)
# se  runif(500)  mi dà uscite che si distribuiscono uniformente in [0,1)
# runif(500)*runif(500)  mi dà uscite più vicine a 0 in quanto il
# prodotto di due numeri minori di 1 è più piccolo di essi.
# Verifichiamolo calcolando le medie:
mean(runif(10000)); mean(runif(10000)*runif(10000))
# In un caso è circa 1/2, nell'altro è circa 1/4 (le cose,
# se ci si pensa, tornano: nel secondo caso faccio il prodotto, per cui
# la media è (1/2)*(1/2) = 1/4).
#
# Il funzionamento di runif è chiaro facendo l'istogramma delle uscite
# Facciamo l'istogramma classificando le uscite nella sequenza di intervalli
# che si ottiene suddividendo [0,1] in intrevalli ampi 0.1
hist(runif(100),seq(0,1,0.1))
# aumentando le prove
hist(runif(10^6), seq(0,1,0.1))
# mettendo sull'asse y le freq. relative invece delle assolute:
hist(runif(10^6), seq(0,1,0.1), probability=TRUE)
#
# Avete visto che usando la parte intera ("floor") posso simulare
# il lancio di un dato. Ecco 30 uscite:
floor(runif(30)*6)+1
# E il lancio di due. 100 lanci:
U1 <- floor(runif(100)*6)+1; U2 <- floor(runif(100)*6)+1
hist(U1+U2, seq(1.5,12.5,1), probability=TRUE)
# e 1 milione:
U1 <- floor(runif(10^6)*6)+1; U2 <- floor(runif(10^6)*6)+1
hist(U1+U2, seq(1.5,12.5,1), probability=TRUE)
# Un modo alternativo:
U <- 0; for (i in 1:2) U <- U + floor(runif(10^6)*6)+1
hist(U, seq(1.5,12.5,1), probability=TRUE)
# Il lancio di un dado equo dava un istogramma piatto.
# Quello di due dadi equi dava un istogramma triangolare
# E quello di 3 (guarda come cambia l'insieme delle possibili
# uscite, ossia i valori che metto in "seq"):
U <- 0; for (i in 1:3) U <- U + floor(runif(10^6)*6)+1
hist(U, seq(2.5,18.5,1), probability=TRUE)
# E se ne lancio 10?
U <- 0; for (i in 1:10) U <- U + floor(runif(10^6)*6)+1
hist(U, seq(9.5,60.5,1), probability=TRUE)
# Assume una forma approsimabile con una curva particolare
# (su cui vi sofferemerete il prossimo anno), che dipende
# dalla media e da un valore particolare (scarto quadratico
# medio):
mean(U); sd(U)
# La funzione si chiama "distribuzione normale":
dn <- function(x) dnorm(x,mean=mean(U),sd=sd(U))
# Eccone il grafico e, sovrapposto, l'istogramma:
plot(dn,10,60,col="red")
hist(U, seq(9.5,60.5,1), probability=TRUE, add=TRUE)
#
# Finiamo qua.
# Buon lavoro per il prossimo anno ...!