R-1/2/3 R-4
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
#