Nothing
#-*- R -*-
##########################################################
### ###
### Script tratti da `Laboratorio di statistica con R' ###
### ###
### Stefano M. Iacus & Guido Masaratto ###
### ###
### CAPITOLO 4 ###
##########################################################
require(labstatR)
### Sez 4.1 CALCOLO DELLE PROBABILITA' E SPAZIO CAMPIONARIO
Omega <- c(1,2,3,4,5,6)
sample(Omega, 15, replace=TRUE)
sample(Omega, 6)
sample(Omega, 6)
sample(Omega, 6)
sample(Omega, 7)
# Scommessa di De Mere
ptm <- proc.time()
gioco1a()
proc.time() - ptm
ptm <- proc.time()
gioco2a()
proc.time() - ptm
ptm <- proc.time()
gioco1()
proc.time() - ptm
ptm <- proc.time()
gioco2()
proc.time() - ptm
system.time( gioco1a() )
system.time( gioco2a() )
system.time( gioco1() )
system.time( gioco2() )
# Probabilita' di compleanni coincidenti
n <- c(5,10,15,20,21,22,23,24,25,30,50,60,
70,80,90,100,200,300,365)
for(i in n)
cat("\n n=",i,"P(A)=",birthday(i))
pbirthday(23, classes = 365, coincident = 4)
qbirthday(prob= 0.5, classes = 365, coincident = 4)
pbirthday(168, classes = 365, coincident = 4)
pbirthday(169, classes = 365, coincident = 4)
qbirthday(prob= 0.5, classes = 365, coincident = 2)
pbirthday(23, classes = 365, coincident = 2)
# calcolo combinatorio
prod(1:4)
choose(4,3)
word <- c("R","O", "M", "A")
sample(word,3)
sample(word,3)
sample(word,3)
prod(1:4)
gamma(5)
prod(1:150)
gamma(151)
exp(lgamma(151))
prod(1:170)
prod(1:171)
gamma(172)
log(gamma(172))
lgamma(200)
lgamma(2000)
lgamma(2000000000)
### Sez 4.2.2 MODELLI MEDIA-VARIANZA
p0 <- 25
pr <- c(0.1,0.2,0.4,0.2,0.1)
p <- c(20,22.5,25,30,40)
pm <- sum(p*pr)
pm
vp <- sum((p-pm)^2*pr)
vp
(pm-p0)/p0
vp/p0^2
x <- c(11,9,25,7,-2)/100
y <- c(-3,15,2,20,6)/100
pxy <- matrix(rep(0,25),5,5)
pxy[1,1] <- 0.2
pxy[2,2] <- 0.2
pxy[3,3] <- 0.2
pxy[4,4] <- 0.2
pxy[5,5] <- 0.2
Rpa(0.1,x,y,pxy)
Rpa(0.5,x,y,pxy)
Rpa(0.7,x,y,pxy)
a <- seq(0,1,0.1)
rr <- Rpa(a,x,y,pxy)
plot(a,rr$Rm,main="rendimento medio", ylab="Rm",type="b")
plot(a,rr$VR,main="varianza del redimento", ylab="Vr",type="b")
plot(sqrt(rr$VR),rr$Rm,main="Trade-off media Varianza",
ylab=expression(E(R[p])),xlab=expression(sigma(R[p])), type="b")
Rp(x,y,pxy)
VRpa <- function(a,vx,vy,mx,my,rxy){
vv <- sqrt(a^2*vx+(1-a)^2*vy+2*a*(1-a)*rxy*sqrt(vx)*sqrt(vy))
mm <- a*mx+(1-a)*my
return(list(vv=vv,mm=mm))
}
a <- seq(0,1,0.1)
mx <- 0.10
my <- 0.08
vx <- 0.0872^2
vy <- 0.0841^2
rr1 <- VRpa(a,vx,vy,mx,my,0.2)
rr2 <- VRpa(a,vx,vy,mx,my,-0.3)
rr3 <- VRpa(a,vx,vy,mx,my,-1)
rr4 <- VRpa(a,vx,vy,mx,my,1)
rr5 <- VRpa(a,vx,vy,mx,my,0)
plot(c(rr4$vv,rr3$vv),c(rr4$mm,rr3$mm),type="l",
xlim=c(0,0.10),lwd=2,xlab=expression(sigma(R[p])),
ylab=expression(E(R[p])))
lines(rr1$vv,rr1$mm)
lines(rr2$vv,rr2$mm)
lines(rr3$vv,rr3$mm)
lines(rr5$vv,rr5$mm,lty=3,lwd=2)
text(0.06,0.09,expression(rho==0))
text(0.04,0.09,expression(rho==+0.2))
text(0.075,0.09,expression(rho==-0.3))
text(0.02,0.095,expression(rho==-1))
text(0.02,0.085,expression(rho==-1))
text(0.095,0.09,expression(rho==+1))
### Sez 4.2.3 ESPERIMENTO DI BERNOULLI E VARIABILI CASUALI DERIVATE
pbinom(3,10,0.3)
pbinom(3,10,0.3, lower.tail=FALSE)
pbinom(3,10,0.3) + pbinom(3,10,0.3, lower.tail=FALSE)
dbinom(3,10,0.3)
k <- 0:10
p <- dgeom(k,1/8)
plot(k,p,type="h",lwd=10)
pnbinom(3,5,0.3)
dnbinom(3,1,0.3)
dgeom(3,0.3)
### Sez 4.2.5 VARIABILE CASUALE DI POISSON
k <- 0:20
p <- dpois(k,lambda=5)
plot(k,p,type="h",lwd=10)
dpois(0,20/60*5)
ppois(10,20/60*10)
### Sez 4.2.10 VARIABILE CASUALE NORMALE
pnorm(3,mean=5,sd=sqrt(3))
pnorm((3-5)/sqrt(3))
curve(dnorm(x,mean=-4),-10,12,ylab="",axes=FALSE)
curve(dnorm(x,mean=7),-10,12,ylab="",add=TRUE)
mu <- 5
sigma <- 2
pnorm(mu+sigma,mean=mu,sd=sigma) - pnorm(mu-sigma,mean=mu,sd=sigma)
pnorm(1) - pnorm(-1)
pnorm(mu+2*sigma,mean=mu,sd=sigma) - pnorm(mu-2*sigma,mean=mu,sd=sigma)
pnorm(2) - pnorm(-2)
pnorm(mu+3*sigma,mean=mu,sd=sigma) - pnorm(mu-3*sigma,mean=mu,sd=sigma)
pnorm(3) - pnorm(-3)
### Sez 4.2.12 VARIABILE CHI-QUADRATO
curve(dchisq(x,df=3),0.,20,ylab="densita'")
curve(dchisq(x,df=5),0.,20,add=TRUE,lty=2)
curve(dchisq(x,df=7),0.,20,add=TRUE,lty=3)
legend(10,0.2,c("gdl = 3","gdl = 5","gdl = 7"), lty=c(1,2,3))
### Sez 4.2.13 VARIABILE t DI STUDENT
curve(dnorm(x),-5,5,ylab="densita'")
curve(dt(x,df=1),-6,6,lty=3,add=TRUE)
legend(2,0.3,c("Z","t"), lty=c(1,3))
### Sez 4.2.14 VARIABILE F DI FISHER
curve(df(x,df1=3,df2=1),0,2,ylab="densita'")
### Sez 4.2.15 VARIABILI CASUALI GAMMA E BETA
gamma(1/2)
sqrt(pi)
gamma(7)
prod(1:6)
beta(1,1)
beta(pi,2*pi)
curve(dbeta(x,1,1),ylim=c(0,3),ylab="densita'")
curve(dbeta(x,0.1,1), add=TRUE,lty=3)
curve(dbeta(x,1,.1), add=TRUE,lty=3)
curve(dbeta(x,.1,.1), add=TRUE,lty=2,lwd=2)
curve(dbeta(x,4,4), add=TRUE,lty=2,lwd=2)
curve(dbeta(x,2,6), add=TRUE,lty=2,lwd=3)
curve(dbeta(x,6,2), add=TRUE,lty=2,lwd=3)
curve(dbeta(x,2,2), add=TRUE,lty=3,lwd=3)
dgamma(1,2,3)
dgamma(1,2,scale=3)
dgamma(1,2,1/3)
# Sez 4.3.1 IL METODO DELL'INVERSIONE
1*(runif(5)<1/3)
a <- runif(5)
a
a<1/3
1*(a<1/3)
as.integer(a<1/3)
sum((runif(10)<1/3))
gen.vc2 <- function(x,p)
x[min(which(cumsum(p)>runif(1)))]
x <- c(-2,3,7,10,12)
p <- c(0.2, 0.1, 0.4, 0.2, 0.1)
y <- numeric(1000)
for(i in 1:1000) y[i] <- gen.vc(x,p)
table(y)/1000
# secondo generatore
for(i in 1:1000) y[i] <- gen.vc2(x,p)
table(y)/1000
# meglio usare la funzione sample di R
y <- sample( c(-2,3,7,10,12), 1000, c(0.2, 0.1, 0.4, 0.2, 0.1), replace = TRUE)
table(y)/1000
lambda <- 0.25
y <- -log(runif(1000))/lambda
hist(y,freq=FALSE)
curve(dexp(x,lambda),0,25,add=TRUE)
### Sez 4.3.2 IL METODO DEL RIFIUTO
u <- runif(3000)
y <- runif(3000)
w <- 256/27*y*(1-y)^3
z <- which(u <= w)
x <- y[z]
length(x)
plot(density(x),main="")
curve(20*x*(1-x)^3,0,1,add=TRUE,lty=2,lwd=2)
x <- rnorm(1000)
plot(density(x))
plot(density(x), main="Dati simulati e vero modello")
curve(dnorm(x), add=TRUE, lty=2)
### Sez 4.4 I PROCESSI STOCASTICI
x <- rbinom(15,1,0.3)
x
plot(x,type="s",main="Processo di Bernoulli", ylab="Spazio degli stati",xlab="tempo")
points(1:15,x)
### Sez 4.4.1 LA PASSEGGIATA ALEATORIA
n <- 50
x <- rbinom(n,1,0.5)
x
x[which(x==0)] <- -1
x
y <- cumsum(x)
plot(1:n,y,type="l", main="passeggiata aleatoria", xlab="tempo",ylab="posizione")
abline(h=0,lty=3)
n <- 500
x <- rbinom(n,1,0.45)
x[which(x==0)] <- -1
y <- cumsum(x)
plot(1:n,y,type="l", main="passeggiata aleatoria",
xlab="tempo",ylab="posizione")
abline(h=0,lty=3)
n <- 500
x <- rbinom(n,1,0.51)
x[which(x==0)] <- -1
y <- cumsum(x)
plot(1:n,y,type="l", main="passeggiata aleatoria", xlab="tempo",ylab="posizione")
abline(h=0,lty=3)
# una barriera
n <- 50000
L <- 40
t <- numeric(100)
t.na <- numeric(100)
for(i in 1:100){
x <- rbinom(n,1,0.5)
x <- 2*x -1
y <- cumsum(x)
t1 <- min(which(y==L))
t2 <- t1
if(t1>n) { t1 <- n; t2 <- NA;}
t[i] <- t1
t.na[i] <- t2
}
mean(t)
mean(t.na,na.rm=TRUE)
# una barriera rifettente
n <- 500
L <- 10
continua <- TRUE
x <- 2*rbinom(n,1,0.5) - 1
x1 <- x
while(continua){
y <- cumsum(x)
bar <- which(y==L+1)
if(length(bar) == 0)
continua = FALSE
else{
h <- min(bar)
x[h] <- -1
}
}
plot(1:n,cumsum(x1),type="l",lty=3,ylab="", xlab="n",ylim=c(-30,30))
lines(1:n,y)
abline(h=L,lty=2)
text(5,L+2,"L")
legend(60,-15,c("libera","riflessa"),lty=c(3,1))
# Passeggiata con due barriere
n <- 1000
L1 <- 10
L2 <- -5
continua <- TRUE
x <- 2*rbinom(n,1,0.5) - 1
x1 <- x
while(continua){
y <- cumsum(x)
bar1 <- which(y==L1+1)
bar2 <- which(y==L2-1)
if( (length(bar1) == 0) & (length(bar2) == 0))
continua = FALSE
else{
h1 <- min(bar1)
h2 <- min(bar2)
if(min(h1,h2) == h1)
x[h1] <- -1
else
x[h2] <- +1
}
}
plot(1:n,cumsum(x1),type="l",lty=3,ylab="",
xlab="n",ylim=c(-30,30))
lines(1:n,y)
abline(h=L1,lty=2)
abline(h=L2,lty=2)
text(5,L1+2,"L1")
text(5,L2-4,"L2")
legend(600,-15,c("libera","riflessa"),lty=c(3,1))
### Sez 4.4.2 CATENTE DI MARKOV
p0 <- c(0,1,0)
P <- matrix(c(0.5,0.5,0.25,0.25,0,0.25,0.25,0.5,0.5),3,3)
p0 %*% P
p0 %*% (P %*% P)
p0 %*% (P %*% P %*% P)
x <- c("P","S","N")
P
Markov("S",15,x,P) -> traj
traj
plot(traj$t,codes(factor(traj$X)),type="s",axes=FALSE, xlab="t",ylab="Che tempo fa")
axis(1)
axis(2,c(1,2,3),levels(factor(traj$X)))
box()
P
P %*% P
P <- matrix( c(1,0.5,0,0.5), 2,2)
P
P %*% P
P %*% P %*% P
# tempo medio di ritorno
P <- matrix( c(0.5,0.5,0.25,0.25,0,0.25,0.25,0.5,0.5), 3,3)
x <- c("P","S","N")
mm <- Markov("S", 10000, x, P)
r1 <- which(mm$X=="P")
n1 <- length(r1)
p1 <- 1/(mean(r1[-1]-r1[-n1]))
p1
r2 <- which(mm$X=="S")
p2 <- 1/(mean( diff(r2) )
r3 <- which(mm$X=="N")
p3 <- 1/(mean( diff(r3) )
p2
p3
### Sez 4.3.3 PROCESSI AUTOREGRESSIVI
n <- 100
lambda <- 0.3
x <- rnorm(n)
y <- numeric(n)
y[1] <- 0
for(i in 2:n) y[i] <- y[i-1] * lambda + x[i]
plot(1:n,y,type="l",xlab="n", ylab=expression(X[n]), main="Modello AR(1)")
curve(0.3^x, 0,4, main=expression(rho(h)==lambda^h), ylab=expression(rho(h)),xlab="h")
library(ts)
acf(y)
### Sez 4.4.4 PROCESSO DI POISSON
n <- 10
x <- rexp(n,rate=1/10)
y <- c(0,cumsum(x))
plot(y,0:n,type="s",xlim=c(0,max(y)),ylim=c(0,n),
xlab="t",ylab="N(t)",main="Processo di Poisson")
n <- 100
t <- 50
x <- rexp(n,rate=1/10)
y <- c(0,cumsum(x))
evt <- max(which(y<t)) - 1
evt
plot(y[0:(evt+2)],0:(evt+1),type="s",xlim=c(0,y[evt+2]),
ylim=c(0,evt+2),xlab="t",ylab="N(t)", main="Processo di Poisson")
abline(v=t,lty=3)
# processo di Poisson non omogeneo
lambda <- 1.1
T <- 20
E <- 0
t <- 0
while(t<T){
t <- t - 1/lambda * log(runif(1))
if( runif(1) < sin(t)/lambda )
E <- c(E, t)
}
length(E)
E
plot(E,0:(length(E)-1),type="s",ylim=c(-4,length(E)))
curve(-3+sin(x),0,20,add=TRUE,lty=2,lwd=2)
lewis(T,sin)
lewis(T,sin,FALSE)
### Sez 4.4.5 PROCESSI DI DIFFUSIONE
n <- 100
T <- 1
dt <- T/n
y <- numeric(n+1)
for(i in 2:(n+1))
y[i] <- y[i-1] + rnorm(1) * sqrt(dt)
plot(seq(0,T,dt),y,type="l",main="Moto browniano", xlab="t",ylab="W(t)")
n <- 100
T <- 1
dt <- T/n
x <- c(0,rnorm(n,sd=sqrt(dt)))
y <- cumsum(x)
n <- 100
T <- 1
dt <- T/n
x0 <- 1
mu <- function(x,t) { -x*t }
sigma <- function(x,t) { x*t }
y <- numeric(n+1)
y[1] <- x0
for(i in 2:(n+1)){
t <- dt*(i-1)
y[i] <- y[i-1] + mu(y[i-1], t) *dt + sigma(y[i-1], t) * rnorm(1,sd=sqrt(dt))
}
plot(seq(0,T,dt),y,type="l",main="Processo di diffusione", xlab="t",ylab="X(t)")
diffusione <- trajectory(1,0,1,mu,sigma,100)
plot(diffusione$t,diffusione$y,type="l")
acf(diffusione$y, main="Processo di diffusione")
# EOF Cap4.R
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.