Attività con R - 5
# Appendice della volta scorsa. # Il grafico in coordinate polari dei millimetri di pioggia mensili # nel corso di tre anni in una data località. # Provate a capire/spiegare come è stato costruito questo diagramma: # ossia, quali commenti inserire per renderlo facilmente comprensibile? # mm <- c(44, 59, 75, 80,96,98,97,99,100,102,104,101, 103,102,108,100,92,71,49,28, 27, 45,66, 81, 92,99, 93, 85, 75,67,61,62, 64, 68,71, 72) ang <- seq(30,30*12*3,30) plot(c(-110,110),c(-110,110),type="n",xlab="",ylab="",axes=FALSE,asp=1) for (i in 1:4) symbols(0,0,circles=i*25, inches=FALSE, add=TRUE, fg="grey60") rad <- function(x) x*pi/180 for(i in 1:6) abline(0,tan(rad(30*i)),lty=3) A <- ang*pi/180 lines(cos(A)*mm,sin(A)*mm) A1 <- A[1:12]; A2 <- A[13:24]; A3 <- A[25:36] mm1 <- mm[1:12]; mm2 <- mm[13:24]; mm3 <- mm[25:36] points(cos(A1)*mm1,sin(A1)*mm1,col="red",pch=20) points(cos(A2)*mm2,sin(A2)*mm2,col="green4",pch=20) points(cos(A3)*mm3,sin(A3)*mm3,col="blue",pch=20) text(cos(A[1]+0.1)*115,sin(A[1]+0.1)*115,"Ge") text(cos(A[4]+0.1)*115,sin(A[4]+0.1)*115,"Ap") text(cos(A[7]+0.1)*115,sin(A[7]+0.1)*115,"Lu") text(cos(A[10]+0.1)*115,sin(A[10]+0.1)*115,"Ot")# (soluzione) # # Come calcolare varianza e scarto quadratico medio dei # dati 13, 15, 18, 22, 25 dati <- c(13,15,18,22,25) V <- function(dati) sum((dati-mean(dati))^2)/length(dati) sqm <- function(dati) sqrt(V(dati)) mean(dati); sqm(dati) # 18.6 4.409082 # # N dei lanci di una moneta equa da effettuare fino ad ottenere l'uscita "testa" # Al 50% esce dopo 1 lancio, al 25% dopo 2, al 12.5% dopo 3, ... # Studio della serie 1/2 + 1/2^2 + 1/2^3 +... n <- 1; s <- 0; for(i in 1:n) s <- s + 1/2^i; s # 0.5 n <- 2; s <- 0; for(i in 1:n) s <- s + 1/2^i; s # 0.75 n <- 10; s <- 0; for(i in 1:n) s <- s + 1/2^i; s # 0.9990234 n <- 100; s <- 0; for(i in 1:n) s <- s + 1/2^i; s # 1 # Potrei studiarla anche così: a <- function(n) 1/2^n # a(n) elemento n-esimo della sommatoria N <- function(n) seq(1,n,1) # N = 1 2 ... n S <- function(n) sum(a(N(n))) # somma a(1)+...a(n) S(1); S(10); S(20); S(40) # 0.5 0.9990234 0.999999 1 # Posso farne il grafico anche usando barplot con l'opzione space: a <- function(n) 1/2^n barplot(a(1:10), space=0) abline(v=(1:10), h=axTicks(2), col="brown",lty=3)
# Ovvero posso usare i comandi seguenti (type="h" sta per "tipo istogramma"): a <- function(n) 1/2^n plot(a,1,10, n=10, type="h", lwd=5, xlab="", ylab="",xlim=c(0,10)) # ovvero: plot(1:10,a(1:10), type="h", lwd=5, xlab="", ylab="",xlim=c(0,10)) abline(v=axTicks(1), h=axTicks(2), col="brown",lty=3)
# La serie 1/2 + 2*1/2^2 + 3*1/2^3 + ... # (la media del numero N dei lanci di una moneta per avere testa) n <- 1; s <- 0; for(i in 1:n) s <- s + i*1/2^i; s n <- 10; s <- 0; for(i in 1:n) s <- s + i*1/2^i; s n <- 100; s <- 0; for(i in 1:n) s <- s + i*1/2^i; s # 0.5 1.988281 2 # Alternativa: a <- function(n) n*1/2^n; S(1); S(10); S(100) # 0.5 1.988281 2 # # La somma di due uscite del generatore di numeri casuali: n <- 1e5; U1 <- runif(n); U2 <- runif(n); mean(U1+U2) # 1.002107 hist(U1+U2, probability=TRUE, col="grey90") abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3) # la "curva limite": lines( c(0,1,2), c(0,1,0) ,lty=2, col="brown", lwd=2)
# # con WolframAlpha # plot 1-abs(x-1) from x=0 to 2 # integrate 1-abs(x-1) from x=0 to 2 # # con WolframAlpha l'integrale di una funzione continua a tratti: # plot floor(x) from x=2 to 4 # integrate floor(x) from x=2 to 4 # # Gli integrali definiti # (lo sqm della distribuzione uniforme) g <- function(x) (x-1/2)^2 integrate(g,0,1) # 0.08333333 with absolute error < 9.3e-16 # Se specifichiamo di vole aver solo il valore: integrate(g,0,1)$value # 0.08333333 # Per avere il risultato sotto forma di frazione: library(MASS); fractions( integrate(g,0,1)$value ) # 1/12 # # con WolframAlpha: # integrate (x-1/2)^2 from x=0 to 1 # # Esempio di Var(X+Y) = Var(X)+Var(Y) con X ed Y indipendenti (lancio # di due dadi equi): n <- 1e7 U1 <- floor(runif(n)*6)+1; U2 <- floor(runif(n)*6)+1 V1 <- sum( (U1-mean(U1))^2/n); V2 <- sum( (U2-mean(U2))^2/n) V <- sum( ((U1+U2)-mean(U1+U2))^2/n) V; V1; V2; V1+V2 # 5.834054 2.917164 2.916623 5.833788 # # Esempio di M(X·Y) = M(X)·M(Y) con X ed Y indipendenti (esiti del lancio # di due dadi equi e del prodotto delle uscite): n <- 1e7 U1 <- floor(runif(n)*6)+1; U2 <- floor(runif(n)*6)+1 m1 <- mean(U1); m2 <- mean(U2) m <- mean(U1*U2) m; m1; m2; m1*m2 # 12.25129 3.500237 3.500425 12.25232 #