lorenzo — Oct 24, 2013, 10:15 PM
# definisco la directory di lavoro
setwd("/Users/lorenzo/Documents/Lavori guidati e corsi matematica/ST410 a.a. 2013-2014/Esercitazione3")
# Esempi di utilizzo di altri dispositivi grafici:
# utilizzare il dispositivo grafico di RStudio:
curve(dnorm(x),-3,3)
# conoscere quale dispositivo grafico e' in uso
dev.cur()
pdf
2
# aprire un altro dispositivo grafico
x11()
# chiudere il dispositivo grafico
dev.off()
pdf
2
# aprire il dispositivo grafico "pdf"
pdf(file="prova.pdf", paper="a4")
# disegnare il grafico nel file
curve(dnorm(x),-3,3)
# chiudere il dispositivo grafico "pdf"
dev.off()
pdf
2
# installare pacchetti aggiuntivi
# ad esempio installiamo il pacchetto actuar che contiene la funzione "mgfnorm" (funzione generatrice dei momenti della normale)
#install.packages("actuar")
library("actuar")
Attaching package: 'actuar'
The following object is masked from 'package:grDevices':
cm
?mgfnorm
mgfnorm(0)
[1] 1
#pdf(file="chi.pdf", paper="a4r")
curve(dchisq(x,df=1),0,10, ylim=c(0,0.5), lwd=3, col="black", main="densita' di probabilita'", xlab="x", ylab="f(x)")
for(i in 1:9) {curve(dchisq(x,df=i),0,10, lwd=3, col=rainbow(9)[i], add=T)}
legend(1,0.5,col=rainbow(9),lty=1, legend=paste(rep("df=",9), 1:9), ncol=3)
#dev.off()
#pdf(file="chi2.pdf", paper="a4r")
curve(dchisq(x,df=100), 50,150, main="densita' di probabilita'", xlab="x", ylab="f(x)", lwd=3)
curve(dnorm(x,mean=100, sd=sqrt(200)), add=T, lty=2, col="grey", lwd=3)
abline(v=100, col="red")
abline(v=98, col="blue")
legend(x=110,y=0.028,col=c("black", "grey"), lty=c(1,2), lwd=3, cex=0.8,legend=c("chi quadrato - df=100.", "normale media=100, var=200"))
#dev.off()
#Esercizio 1
n<-5 # numero variabili normali indipendenti, di media 0 e varianza 1
somma<-function(x) sum((x-mean(x))^2)
y<-rep(0,1000)
for(i in 1:1000) y[i]<-somma(rnorm(n))
hist(y, freq=F, lty=1, col="lightgrey")
lines(density(y), col="black")
curve(dchisq(x, df=n-1), add=T, col="red")
#quartz()
to<-10
step<-1
curve(dt(x,df=1), main="", xlab="x", ylab="f(x)", lwd=3, ylim=c(0,0.4), xlim=c(-3,3))
curve(dnorm(x), add=T, col="black", lty=2, lwd=1)
Sys.sleep(2)
for(i in seq(1,to,step)) {{curve(dt(x,df=i), lty=1, lwd=3, col=rainbow(to)[i], add=T);Sys.sleep(1)}}
# Esercizio 2
# Z ha distribuzione t di student: (X_1 + X_2) e (X_1-X_2) sono normali, indipendenti
# di media 0 e varianza 2. Dividendo numeratore e denominatore per sqrt(2)
# Z=(X_1+X_2)/sqrt((X_1-X_2)^2) si puo' vedere come
# rapporto tra una normale di media 0 e varianza 1 e il quadrato di una normale
# di media 0 e varianza 1, i.e. una chi-quadro con un grado di liberta'.
# Percio' Z e' una t di student con un grado di liberta'.
x1<-rnorm(300, mean=0, sd=1)
x2<-rnorm(300, mean=0, sd=1)
z<-(x1+x2)/(sqrt((x1-x2)^2))
plot(density(z), xlim=c(-10,10))
curve(dt(x, df=1), add=T, col="red")
hist(z,lty=2, border="lightgrey", breaks=100,add=T, freq=F)
# Esercizio 3
statura<-read.csv("statura.csv")
str(statura)
'data.frame': 1802 obs. of 2 variables:
$ sesso : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
$ altezze: num 181 174 164 165 183 ...
table(statura$sesso)
F M
801 1001
tapply(statura$altezze, statura$sesso, mean)
F M
160.0 175.3
tapply(statura$altezze, statura$sesso, sd)
F M
7.231 11.079
hist(statura$altezze, freq=F, ylim=c(0,0.052), border="black", main="")
hist(statura$altezze[statura$sesso=="F"], freq=F, add=T, lty=2, border="pink")
hist(statura$altezze[statura$sesso=="M"], freq=F, add=T, lty=2, border="blue")
lines(density(statura$altezze), col="black", lwd=3)
lines(density(statura$altezze[statura$sesso=="F"]), col="pink", lwd=3)
lines(density(statura$altezze[statura$sesso=="M"]), col="blue", lwd=3)
maschi<-statura$altezze[statura$sesso=="M"]
str(maschi)
num [1:1001] 181 174 164 165 183 ...
(sqrt(1001)*sqrt(1000)*(mean(maschi)-174.5))/(sqrt(sum((maschi-mean(maschi))^2)))
[1] 2.362
#pdf(file="test.pdf", paper="a4r")
curve(dt(x,df=1000), -3,3, yaxt="n", xlab="", bty="n", ylab="" )
curve(dt(x,df=1000), -3,qt(0.05, df=1000), add=T, type="h", col="lightgrey", n=1000)
curve(dt(x,df=1000), -3,qt(0.025, df=1000), add=T, type="h", col="lightblue", n=1000)
curve(dt(x,df=1000), -3,qt(0.005, df=1000), add=T, type="h", col="blue", n=1000)
curve(dt(x,df=1000), qt(0.95, df=1000),3, add=T, type="h", col="lightgrey", n=1000)
curve(dt(x,df=1000), qt(0.975, df=1000),3, add=T, type="h", col="lightblue", n=1000)
curve(dt(x,df=1000), qt(0.995, df=1000),3, add=T, type="h", col="blue", n=1000)
points(2.361658,dt(2.361658,df=1000), col="red", pch=16)
legend(x=1.5,y=0.4,legend=c("1%", "5%", "10%"),col=c("blue", "lightblue", "lightgrey"), lwd=c(3,3,3))
lines(x=c(-2.361658,2.361658), y=dt(c(2.361658,2.361658),df=1000), type="h",lwd=2,lty=2, col="red")
(1-pt(2.361658,df=1000))*2
[1] 0.01838
#dev.off()
t.test(maschi, mu=174.5)
One Sample t-test
data: maschi
t = 2.362, df = 1000, p-value = 0.01838
alternative hypothesis: true mean is not equal to 174.5
95 percent confidence interval:
174.6 176.0
sample estimates:
mean of x
175.3