# alcune funzioni richiamate sono nel package rmf
# che va quindi caricato come tale o come workspace di R
library(rmf)

# Uniforme continua
# Figura (pag. 125)
xx <- c(120, 120, 140, 140)
yy <- c(0, 0.5, 0.5, 0)
plot(c(110,150), c(0,1), type="n", las=1, bty="l", ylab=expression(italic(f(x))),
     xlab=expression(italic(x)), xaxs="i", yaxs="i", yaxt = "n")
axis(2, at=0.5, labels=expression(frac(1,20)), las=1)
polygon(xx, yy, col = gray(0.9))

# Figura (pag. 126)
xx2 <- c(5, 5, 7.5, 7.5, NA,
         7.5, 7.5, 12.5, 12.5, NA,
         12.5, 12.5, 15, 15)
yy2 <- c(0, 0.5, 0.5, 0, NA,
         0, 0.5, 0.5, 0, NA,
         0, 0.5, 0.5, 0)
plot(c(0,20), c(0,1), type="n", las=1, bty="l", ylab=expression(italic(f(x))),
     xlab=expression(italic(x)), xaxs="i", yaxs="i", yaxt = "n", xaxt = "n")
axis(2, at=0.5, labels=expression(frac(1,10)), las=1)
axis(1, at=c(5,7.5,12.5,15))
polygon(xx2, yy2, col=c("white", gray(0.9), "white"), lty="dashed")
segments(5, 0.5, 15, 0.5)

punif(7.5,5,15)
punif(12.5,5,15)
punif(12.5,5,15) - punif(7.5,5,15)

# Figura 3.1 (pag. 127)
set.seed(123456)
x <- runif(100000,5,15)
c(mean(x),var(x))
# -------------------------------------------------------------------------------
# (a)
tmp <- hist(x,main="",xlab="",prob=TRUE)
box()
abline(h=0.1,lty=3,lwd=2)
# (b)
p <- seq(0.001,0.999,0.001)
qc <- quantile(x,p)
qt <- (qunif(p,5,15))
plot(qt,qc,xlab="Quantili teorici",ylab="Quantili campionari")
# -------------------------------------------------------------------------------

# pag. 128
p <- seq(0.1,0.9,0.1)
qc <- quantile(x,p)
round(qc,0)
qt <- (qunif(p,5,15))
round(qt,0)

# Distribuzione triangolare (pag. 129)
# č necessario il package triangle
library(triangle)
set.seed(654321)
x1 <- runif(100000,5,15)
x2 <- runif(100000,5,15)
x <- x1 + x2
c(mean(x),var(x))
p <- seq(0.1,0.9,0.1)
qc <- quantile(x,p)
round(qc,0)
qt <- (qtriangle(p,10,30,20))
round(qt,0)

# Figura 3.2 (pag. 130)
# -------------------------------------------------------------------------------
# (a)
hist(x,main="",xlab="", prob=TRUE)
box()
lines(c(10,20),c(0,0.1),lwd=2,lty=1)
lines(c(20,30),c(0.1,0),lwd=2,lty=1)
# (b)
p <- seq(0.001,0.999,0.001)
qc <- quantile(x,p)
qt <- (qtriangle(p,10,30,20))
plot(qt,qc,xlab="Quantili teorici",ylab="Quantili campionari")
# -------------------------------------------------------------------------------

# Figura 3.3 (pag. 131)
set.seed(135246)
X <- runif(100000*18,5,15)
Y <- matrix(X,ncol=18)
x <- rowSums(Y)
c(mean(x),var(x))
# -------------------------------------------------------------------------------
# (a)
hist(x,main="",xlab="",prob=TRUE)
box()
# (b)
hist(x,breaks=100,main="",xlab="",prob=TRUE)
box()
# -------------------------------------------------------------------------------

# Esponenziale (pag. 133)
pexp(6,0.068)
pexp(18,0.068)
pexp(18,0.068) - pexp(6,0.068)

# Figura (pag. 133)
# č necessario il package latex2exp
library(latex2exp)
lambda <- 0.068
# Grafico vuoto e assi
plot(
  NA, xlim = c(0, 38), ylim = c(0, 0.09), axes = FALSE, main = '',
  xlab = 't', ylab = latex2exp::TeX('$f(t)$')
)
axis(1, at = c(0, 38), labels = c('', ''), lwd.ticks = 0)
axis(1, at = (0:6) * 6, lwd = 0, lwd.ticks = 1)
axis(2, at = c(0, 0.1), labels = c('', ''), lwd.ticks = 0)
axis(2, lwd = 0, lwd.ticks = 1)
# Sezione 1
x <- (0:100) / 100 * 6
polygon(c(x, rev(x)), c(dexp(x, lambda), rep(0, 101)), border = NA, col = 'gray50')
# Sezione 2
x <- 6 + (0:100) / 100 * (18 - 6)
polygon(c(x, rev(x)), c(dexp(x, lambda), rep(0, 101)), border = NA, col = 'gray70')
# Sezione 3
x <- 18 + (0:100) / 100 * (38 - 18)
polygon(c(x, rev(x)), c(dexp(x, lambda), rep(0, 101)), border = NA, col = 'gray90')
# Densita'
curve(dexp(x, lambda), 0, 38, add = TRUE)
# Annotazione 1
arrows(x0 = 3, y0 = 0.07, x1 = 3, y1 = 0.6 * dexp(3, lambda), length = 0.08, lwd = 1.5)
lines(x = c(3, 6), y = c(0.07, 0.07), lwd = 1.5)
depo <- latex2exp::TeX('$P(T\\leq 6)$')
text(x = 6, y = 0.07 - 0.001, depo, pos = 4)
# Annotazione 2
arrows(x0 = 12, y0 = 0.05, x1 = 12, y1 = 0.6 * dexp(12, lambda), length = 0.08, lwd = 1.5)
lines(x = c(12, 15), y = c(0.05, 0.05), lwd = 1.5)
depo <- latex2exp::TeX('$P(6\\leq T\\leq 18)$')
text(x = 15, y = 0.05 - 0.001, depo, pos = 4)
# Linee verticali
lines(c(6, 6), c(0, dexp(6, lambda)))
lines(c(18, 18), c(0, dexp(18, lambda)))

# pag. 133
set.seed(123456)
x <- rexp(100000,0.0462)
c(median(x),mean(x))

# Figura 3.4 (pag. 134)
# -------------------------------------------------------------------------------
# (a)
hist(x,breaks=100,prob=TRUE,xlim=c(0,100),main="",xlab="Tempo (in minuti)")
box()
curve(dexp(x,0.0462),add=TRUE)
# (b)
p <- seq(0.1,0.9,0.1)
qc <- quantile(x,p)
qt <- (qexp(p,0.0462))
# round(qc,0); round(qt,0)
p <- seq(0.001,0.999,0.001)
qc <- quantile(x,p)
qt <- (qexp(p,0.0462))
plot(qt,qc,xlab="Quantili teorici",ylab="Quantili campionari")
# -------------------------------------------------------------------------------

# pag. 134
lambda=5
t <- rexp(10,lambda)
T <- cumsum(t)
round(rbind(t,T),2)

# pag. 135
lambda <- 1.5
nrep <- 10000
nexp <- 20
ris <- numeric(nrep)
tot <- matrix(nrow=nrep,ncol=nexp)
set.seed(123456)
for (i in 1:nrep){
  t <- rexp(nexp,lambda)
  T <- cumsum(t)
  if (sum(t) < 1) stop("Numero di tempi insufficiente.")
  tot[i,] <- t
  x <- length(which(T <= 1))
  ris[i] <- x
}
c(median(tot),-log(0.5)/lambda)
m <- mean(tot)
c(m,1/m)
mean(ris)

fo <- table(ris)
fa <- dpois(c(0:8),lambda)*nrep
rbind(fo,fa=round(fa,1))

# Gamma (pag. 136)
# Figura 3.5
# -------------------------------------------------------------------------------
# (a)
lambda <- 0.0462
n <- 2
set.seed(123456)
u <- rexp(100000*n,lambda)
v <- matrix(u, ncol=n)
x <- rowSums(v)
c(mean(x),var(x))
c(n/lambda,n/lambda^2)
median(x)
c(mean(x),var(x))
c(n/lambda,n/lambda^2)
median(x)
hist(x,breaks=100,prob=TRUE,xlim=c(0,200),main="",xlab="Tempo (in minuti)")
box()
curve(dgamma(x,shape=2,rate=0.0462),add=TRUE,lwd=2)
# (b)
lambda <- 0.0462
n <- 10
set.seed(123456)
u <- rexp(100000*n,lambda)
v <- matrix(u, ncol=n)
x <- rowSums(v)
c(mean(x),var(x))
c(n/lambda,n/lambda^2)
median(x)
hist(x,breaks=100,prob=TRUE,xlim=c(0,500),main="",xlab="Tempo (in minuti)")
box()
curve(dgamma(x,shape=10,rate=0.0462),add=TRUE,lwd=2)
# -------------------------------------------------------------------------------

# Weibull (pag. 141)
# Figura 3.6 (pag. 141)
# -------------------------------------------------------------------------------
H <- function(x,k,h) k*h*(h*x)^(k-1)
# (a)
plot(c(0,3),c(0,5),type="n",xlab="tempo",ylab="azzardo")
curve(H(x,0.5,1),add=TRUE,lty=2)
curve(H(x,  1,1),add=TRUE,lty=3)
curve(H(x,1.5,1),add=TRUE,lty=1)
curve(H(x,  2,1),add=TRUE,lty=4)
curve(H(x,2.5,1),add=TRUE,lty=5)
legend(0,5,c("k=0.5","k=1","k=1.5","k=2","k=2.5"),lty=c(2,3,1,4,5))
# (b)
plot(c(0,3),c(0,2),type="n",xlab="tempo",ylab="azzardo")
h <- 0.5
curve(H(x,0.5,h),add=TRUE,lty=2)
curve(H(x,  1,h),add=TRUE,lty=3)
curve(H(x,1.5,h),add=TRUE,lty=1)
curve(H(x,  2,h),add=TRUE,lty=4)
curve(H(x,2.5,h),add=TRUE,lty=5)
legend(0.2,2,c("k=0.5","k=1","k=1.5","k=2","k=2.5"),lty=c(2,3,1,4,5))
# -------------------------------------------------------------------------------

# pag. 142
set.seed(135246)
n <- 10000
x <- rweibull(n,shape=0.8,scale=2)

t <- sort(x)
F <- c(1:n)/n
S <- 1 - F

S <- S[1:(n-1)]
t <- t[1:(n-1)]
lml <- log(-log(S))

fit <- lm(lml ~ log(t))
(b <- fit$coefficients)
exp(b[1]/b[2])

# Figura 3.7 (pag. 143)
# -------------------------------------------------------------------------------
# (a)
ylb <- expression(hat(S)(t))
plot(t,S,type="s",xlab="t",ylab="")
title(ylab=ylb, line=2.5)
# (b)
ylb <- expression(log(-log(hat(S)(t))))
plot(log(t),lml,xlab="log(t)",ylab="")
title(ylab=ylb, line=2.5)
# -------------------------------------------------------------------------------

# Normale (pag. 144)
# Figura 3.8 (pag. 144)
# -------------------------------------------------------------------------------
# (a)
curve(dnorm(x,90,5),70,110,xlab="x",ylab="f(x)")
abline(h=0)
points(90,0)
points(c(85,95),c(0,0),pch=19)
y <- dnorm(85,90,5)
lines(c(85,85),c(0,y),lty=3)
lines(c(95,95),c(0,y),lty=3)
# (b)
p <- c(1:9999)/10000
plot(c(75,105),c(0,1),type="n",xlab="x",ylab="F(x)")
x <- qnorm(p,90,5); lines(x,p,lwd=lwd)
lines(c(70,90),c(0.5,0.5),lty=3,lwd=2)
lines(c(90,90),c(-1,0.5),lty=3,lwd=2)
# -------------------------------------------------------------------------------

# Figura 3.9 (pag. 145)
# -------------------------------------------------------------------------------
lwd <- 2
# (a)
plot(c(120,220),c(0,0.04),type="n",xlab="x",ylab="f(x)")
Normale(mu=170,sigma=10,add=TRUE,lwd=lwd)
Normale(mu=180,sigma=10,add=TRUE,lty=2,lwd=lwd)
Normale(mu=160,sigma=10,add=TRUE,lty=3,lwd=lwd)
points(c(160,170,180),c(0,0,0),pch=19)
# (b)
plot(c(100,240),c(0,0.04),type="n",xlab="x",ylab="f(x)")
Normale(mu=170,sigma=10,add=TRUE,lwd=lwd)
Normale(mu=170,sigma=15,add=TRUE,lty=2,lwd=lwd)
Normale(mu=170,sigma=20,add=TRUE,lty=3,lwd=lwd)
# -------------------------------------------------------------------------------

# Figura 3.10 (pag. 146)
# -------------------------------------------------------------------------------
lwd <- 2
# (a)
p <- c(1:9999)/10000
plot(c(120,220),c(0,1),type="n",xlab="x",ylab="F(x)")
x <- qnorm(p,170,10); lines(x,p,lwd=lwd)
x <- qnorm(p,180,10); lines(x,p,lty=2,lwd=lwd)
x <- qnorm(p,160,10); lines(x,p,lty=3,lwd=lwd)
# (b)
p <- c(1:9999)/10000
plot(c(120,220),c(0,1),type="n",xlab="x",ylab="F(x)")
x <- qnorm(p,170,10); lines(x,p,lwd=lwd)
x <- qnorm(p,170,15); lines(x,p,lty=2,lwd=lwd)
x <- qnorm(p,170,20); lines(x,p,lty=3,lwd=lwd)
# -------------------------------------------------------------------------------

# le tavole della normale standard (pag. 147)
z <- c(0:399)/100
p <- pnorm(z)
tbl <- matrix(p,ncol=10,byrow=TRUE)
nr <- nrow(tbl)+1
nc <- ncol(tbl)+1
out <- matrix(NA,nr,nc)
out[2:nr,2:nc] <- round(tbl,5)
out[2:nr,1] <- c(0:39)/10
out[1,2:nc] <- round(c(0:9)/100,2)
out

# Figura 3.11 (pag. 149)
# -------------------------------------------------------------------------------
par(mfrow=c(3,3))
plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(a)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=-5,xa=1.96,0,1,color="red")

plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(b)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=1.96,xa=5,0,1,color="red")

plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(c)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=-5,xa=-1.96,0,1,color="red")

plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(d)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=-1.96,xa=5,0,1,color="red")

plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(e)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=-1.96,xa=1.96,0,1,color="red")

plot(c(-3.5,3.5),c(0,0.4),type="n",xlab="",ylab="",main="(f)")
curve(dnorm(x),-3.5,3.5,add=TRUE)
abline(h=0)
FillNorm(xda=-5,xa=-1.96,0,1,color="red")
FillNorm(xda=1.96,xa=5,0,1,color="red")

m <- 170; s <- 10
plot(c(170-3.5*s,170+3.5*s),c(0,0.04),type="n",xlab="",ylab="",main="(g)")
curve(dnorm(x,m,s),170-3.5*s,170+3.5*s,add=TRUE)
abline(h=0)
FillNorm(xda=-5,xa=m+1.96*s,m,s,color="red")

m <- 170; s <- 10
plot(c(170-3.5*s,170+3.5*s),c(0,0.04),type="n",xlab="",ylab="",main="(h)")
curve(dnorm(x,m,s),170-3.5*s,170+3.5*s,add=TRUE)
abline(h=0)
FillNorm(xda=m-1.96*s,xa=m+1.96*s,m,s,color="red")

m <- 170; s <- 10
plot(c(170-3.5*s,170+3.5*s),c(0,0.04),type="n",xlab="",ylab="",main="(i)")
curve(dnorm(x,m,s),170-3.5*s,170+3.5*s,add=TRUE)
abline(h=0)
FillNorm(xda=-5,xa=m-1.96*s,m,s,color="red")
FillNorm(xda=m+1.96*s,xa=300,m,s,color="red")
par(mfrow=c(1,1))
# -------------------------------------------------------------------------------

# pag. 150
pnorm(1.96)
pnorm(1.96,lower.tail=FALSE)
pnorm(-1.96)
pnorm(-1.96,lower.tail=FALSE)
pnorm(1.96)-pnorm(-1.96)
pnorm(-1.96)+pnorm(1.96,lower.tail=FALSE)

ProbNorm(a=1.96,color="red")
ProbNorm(da=1.96,color="red")
ProbNorm(a=-1.96,color="red")
ProbNorm(da=-1.96,color="red")
ProbNorm(da=-1.96,a=1.96,color="red")

# pag. 153
qnorm(c(0.25,0.5,0.75),mean=0,sd=1)
qnorm(c(0.25,0.5,0.75),mean=170,sd=10)

# approssimazione alla normale e correzione per la continuitā (pag. 154)
dbinom(20,size=50,prob=0.5)
pnorm(20.5,25,sqrt(12.5))-pnorm(19.5,25,sqrt(12.5))

(p <- sum(dbinom(c(0:14),size=50,prob=0.5)))
(pa <- pnorm(14,25,sqrt(12.5))); (pa-p)/p
(pa <- pnorm(15,25,sqrt(12.5))); (pa-p)/p
(pa <- pnorm(14.5,25,sqrt(12.5))); (pa-p)/p

(p <- sum(dbinom(c(0:15),size=50,prob=0.5)))
(pa <- pnorm(15,25,sqrt(12.5))); (pa-p)/p
(pa <- pnorm((15.5),25,sqrt(12.5))); (pa-p)/p

(p <- sum(dbinom(c(31:50),size=50,prob=0.5)))
(pa <- pnorm(30.5,25,sqrt(12.5),lower.tail=FALSE)); (pa-p)/p

(p <- sum(dbinom(c(30:50),size=50,prob=0.5)))
(pa <- pnorm(29.5,25,sqrt(12.5),lower.tail=FALSE)); (pa-p)/p

# ancora sulle fluttuazioni (pag. 156)
# Figura 3.12 (pag. 156)
# -------------------------------------------------------------------------------
# (a)
n <- 100
p <- 1/2
x <- c(0:n)
y <- dbinom(x, n, p)
plot(x, y, type = "h", ylab = "p", xlab="x")
# (b)
n <- 1000
p <- 1/2
x <- c(0:n)
y <- dbinom(x, n, p)
plot(x, y, type = "h", ylab = "p", xlab="x")
# -------------------------------------------------------------------------------

# pag. 156
n <- 2; x <- qbinom(c(0.005,0.995),10^n,0.5); x
(x[2] - x[1])/10^n

n <- 3; x <- qbinom(c(0.005,0.995),10^n,0.5); x
(x[2] - x[1])/10^n

n <- 4; x <- qbinom(c(0.005,0.995),10^n,0.5); x
(x[2] - x[1])/10^n

n <- 6; x <- qbinom(c(0.005,0.995),10^n,0.5); x
(x[2] - x[1])/10^n

n <- 10^9; m <- n/2; s <- sqrt(n)/2
p <- 10^(-6)
q <- qnorm(c(p/2,1-p/2),m,s)
(q[2]-q[1])/n

# Figura 3.13 (pag. 158)
# -------------------------------------------------------------------------------
qNor <- function(n,alfa=0.001){
      m <- 10^n/2
      v <- 10^n/4
      s <- sqrt(v)
      x <- qnorm(c(alfa/2,1-alfa/2),m,s)
      dlt <- x[2] - x[1]
      dp <- (x[2] - x[1])/10^n
      c(dlt,dp)
}

n <- c(10:28)
y <- numeric(length(n))
k <- 0
for (i in n){
      k <- k + 1
      y[k] <- qNor(i,10^(-6))[2]
}

# (a)
plot(n,y,xlab="k",ylab="ampiezza percentuale")
# (b)
plot(n,log10(y),xlab="k",ylab="ampiezza percentuale (logaritmo in base 10)")
# -------------------------------------------------------------------------------

# pag. 158
fit <- lm(log10(y) ~ n)
fit$coefficients
qnorm(10^(-6)/2,lower.tail=FALSE)

# Log-normale (pag. 159)
mu <- 0; sigma <- 1
set.seed(123456)
y <- rnorm(100000,mu,sigma)
c(mean(y),sd(y))
x <- exp(y)
c(mean(x),median(x))
c(mean(x),exp(mu+sigma^2/2))
c(median(x),exp(mu))
v <- exp(2*mu+sigma^2) * (exp(sigma^2)-1)
c(var(x),v)
# Figura 3.14a (pag. 160)
# -------------------------------------------------------------------------------
hist(x,breaks=1000,xlim=c(0,6),prob=TRUE,xlab="",main="")
box()
curve(dlnorm(x,0,1),0,6,add=TRUE)
# -------------------------------------------------------------------------------

# pag. 160
mu <- 1; sigma <- 0.5
set.seed(123456)
y <- rnorm(100000,mu,sigma)
x <- exp(y)
c(mean(x),exp(mu+sigma^2/2))
c(median(x),exp(mu))
v <- exp(2*mu+sigma^2) * (exp(sigma^2)-1)
c(var(x),v)
# Figura 3.14b (pag. 160)
# -------------------------------------------------------------------------------
hist(x,breaks=100,xlim=c(0,10),prob=TRUE,xlab="",main="")
box()
curve(dlnorm(x,1,0.5),0,10,add=TRUE)
# -------------------------------------------------------------------------------

# Figura 3.15 (pag. 161)
# -------------------------------------------------------------------------------
# (a)
set.seed(654321)
n <- 10
y <- runif(100000*n,0.5,1.5)
Y <- matrix(y,ncol=n)
x <- rowSums(Y)
c(mean(x),var(x))
hist(x,breaks=100,xlim=c(6,14),prob=TRUE,xlab="",ylab="Densitā",main="")
box()
curve(dnorm(x,10,sqrt(10/12)),add=TRUE)
# (b)
x <- apply(Y,1,prod)
hist(x,breaks=200,xlim=c(0,4),prob=TRUE)
box()
m <- mean(log(x))
s <- sd(log(x))
curve(dlnorm(x,m,s),0,4,add=TRUE,xlab="",main="")
# -------------------------------------------------------------------------------

# Figura 3.16 (pag. 163)
# -------------------------------------------------------------------------------
# occorre leggere i dati indicando la cartella che li contiene
# setwd("e:/!zero/r/3/")
df <- read.table("comuni.txt")
x <- df$x
(m <- mean(log(x)))
(s <- sd(log(x)))
p <- c(1:9)/10
(qo <- round(quantile(x,p),0))
(qa <- round(qlnorm(p,m,s),0))
# (a)
hist(x,breaks=5000,prob=TRUE,xlim=c(0,20000),xlab="Popolazione residente (F+M)",main="")
box()
curve(dlnorm(x,m,s),0,20000,add=TRUE)
# (b)
p <- c(1:95)/100
qo <- round(quantile(x,p),0)
qa <- round(qlnorm(p,m,s),0)
plot(qa,qo,xlab="Quantili teorici",ylab="Quantili campionari")
abline(0,1)
# -------------------------------------------------------------------------------

# Figura 3.17 (pag. 165)
# -------------------------------------------------------------------------------
par(mfrow=c(2,2))
curve(dchisq(x,1),0,5,xlab="",ylab="",main="n = 1")
curve(dchisq(x,2),0,10,xlab="",ylab="",main="n = 2")
curve(dchisq(x,3),0,10,xlab="",ylab="",main="n = 3")
curve(dchisq(x,100),50,150,xlab="",ylab="",main="n = 100")
curve(dnorm(x,100,sqrt(200)),add=TRUE,lty=3,lwd=2)
par(mfrow=c(1,1))
# -------------------------------------------------------------------------------

# t di Student (pag. 165)
set.seed(123456)
x <- rt(1000000,1); mean(x)
x <- rt(1000000,1); mean(x)
x <- rt(1000000,1); mean(x)
x <- rt(1000000,1); mean(x)
x <- rt(1000000,1); mean(x)

set.seed(123456)
x <- rt(1000000,2); mean(x)
x <- rt(1000000,2); mean(x)
x <- rt(1000000,2); mean(x)
x <- rt(1000000,2); mean(x)
x <- rt(1000000,2); mean(x)

set.seed(123456)
x <- rt(1000000,2); var(x)
x <- rt(1000000,2); var(x)
x <- rt(1000000,2); var(x)
x <- rt(1000000,2); var(x)
x <- rt(1000000,2); var(x)

set.seed(123456)
x <- rt(1000000,3); var(x)
x <- rt(1000000,12); var(x)
x <- rt(1000000,102); var(x)
x <- rt(1000000,1002); var(x)

# Figura 3.18 (pag. 167)
# -------------------------------------------------------------------------------
par(mfrow=c(2,2))
curve(dnorm(x),-4,4,lty=3,lwd=2,xlab="",ylab="",main="n = 1")
curve(dt(x,1),-4,4,add=TRUE)
curve(dnorm(x),-4,4,lty=3,lwd=2,xlab="",ylab="",main="n = 2")
curve(dt(x,2),-4,4,add=TRUE)
curve(dnorm(x),-4,4,lty=3,lwd=2,xlab="",ylab="",main="n = 5")
curve(dt(x,5),-4,4,add=TRUE,main="n = 5")
curve(dnorm(x),-4,4,lty=3,lwd=2,xlab="",ylab="",main="n = 10")
curve(dt(x,10),-4,4,add=TRUE,main="n = 2")
par(mfrow=c(1,1))
# -------------------------------------------------------------------------------

# La famiglia esponenziale e i modelli lineari generalizzati
#
# Esempio 3.1 (pag. 170)
set.seed(654321)
a <- -85
b <- 0.88
x <- seq(150,190,5)
n <- length(x)
m <- a + b*x
e <- rnorm(5*n,0,10)
e <- matrix(e,nrow=n)
y <- m + e
Y <- as.vector(y)
X <- rep(x,5)
fit <- lm(Y ~ X)
fit$coefficients
confint(fit)[2,]

# Figura 3.19a (pag. 171)
# -------------------------------------------------------------------------------
plot(c(150,190),c(20,110),type="n",xlab="altezza (cm)",ylab="peso (kg)")
title("(a)")
k <- 0
for (i in x){
      k <- k + 1
      xx <- rep(i,ncol(y))
      points(xx,y[k,])
}
abline(a,b,lwd=1)
abline(fit, lwd=2, lty=3)
# -------------------------------------------------------------------------------

# Esempio 3.2 (pag. 174)
x <- seq(1.7,2.3,0.05)
pl <- -40 + 20*x
p <- exp(pl)/(1+exp(pl))
k <- length(p)
y <- numeric(k)
n <- rep(100,k)
set.seed(135246)
for (i in 1:k) y [i] <- rbinom(1,n[i],p[i])
rbind(x,y)
fit <- glm(y/n ~ x, family=binomial("logit"), weights=n)
fit$coefficients

# Figura 3.19b (pag. 171)
# -------------------------------------------------------------------------------
plot(x,y/n, xlab="dose", ylab="probabilitā", main="(b)")
# curva vera
xx <- seq(1.7,2.3,0.005)
pl <- -40 + 20*xx
ph <- exp(pl)/(1 + exp(pl))
lines(xx,ph)
# osservata
b <- fit$coefficients
pl <- b[1] + b[2]*x
ph <- exp(pl)/(1 + exp(pl))
lines(x,ph,lty=3,lwd=2)
# -------------------------------------------------------------------------------

# Esempio 3.3 (pag. 175)
y <- c(12196,58535)
n <- c(1788800,402516)
x <- c(0,1)
fit <- glm(y/n ~ x, family=binomial("logit"), weights=n)
fit$coefficients
exp(fit$coefficients)

p <- y/n
p[1]/(1 - p[1])
b <- fit$coefficients
exp(b[1])

p <- y/n
p[2]/(1 - p[2])
b <- fit$coefficients
exp(b[1] + b[2])

odds <- p/(1-p)
odds[2]/odds[1]
exp(b[2])
exp(confint(fit))

# Esempio 3.4 (pag. 178)
x <- seq(1.8,2.06,0.02)
pl <- -40 + 20*x
p <- 1 - exp(-exp(pl))
k <- length(p)
y <- numeric(k)
n <- rep(100,k)
set.seed(246135)
for (i in 1:k) y [i] <- rbinom(1,n[i],p[i])
rbind(x,y)
fit <- glm(y/n ~ x, family=binomial("cloglog"), weights=n)
fit$coefficients

# Figura 3.19c (pag. 171)
# -------------------------------------------------------------------------------
plot(x,y/n, xlab="dose", ylab="probabilitā", main="(c)")
# curva vera
xx <- seq(1.7,2.3,0.005)
pl <- -40 + 20*xx
ph <- 1 - exp(-exp(pl))
lines(xx,ph)
# stimata
b <- fit$coefficients
pl <- b[1] + b[2]*x
ph <- 1 - exp(-exp(pl))
lines(x,ph,lty=3,lwd=2)
# con un modello logistico
fit <- glm(y/n ~ x, family=binomial("logit"), weights=n)
fit$coefficients
# stimata
b <- fit$coefficients
pl <- b[1] + b[2]*x
ph <- exp(pl)/(1 + exp(pl))
lines(x,ph,lty=2,lwd=2)
# -------------------------------------------------------------------------------

# le due analisi precedenti sulla probabilitā di insuccesso
x <- seq(1.7,2.3,0.05)
pl <- -40 + 20*x
p <- exp(pl)/(1+exp(pl))
k <- length(p)
y <- numeric(k)
n <- rep(100,k)
set.seed(135246)
for (i in 1:k) y[i] <- rbinom(1,n[i],p[i])
fit <- glm(y/n ~ x, family=binomial("logit"), weights=n)
fit$coefficients
fit <- glm((n-y)/n ~ x, family=binomial("logit"), weights=n)

x <- seq(1.8,2.06,0.02)
pl <- -40 + 20*x
p <- 1 - exp(-exp(pl))
k <- length(p)
y <- numeric(k)
n <- rep(100,k)
set.seed(246135)
for (i in 1:k) y [i] <- rbinom(1,n[i],p[i])
fit <- glm(y/n ~ x, family=binomial("cloglog"), weights=n)
fit$coefficients
fit <- glm((n-y)/n ~ x, family=binomial("cloglog"), weights=n)
fit$coefficients

# Figura a pag. 180
x <- seq(1.7,2.3,0.001)
pl <- -40 + 20*x
p <- exp(pl)/(1+exp(pl))
plot(x,p,type="l",xlab="probabilitā")
pp <- 1 - exp(-exp(pl))
lines(x,pp,lty=3,lwd=2)
legend(1.7,1,c("Link logit","Link cloglog"),lty=c(1,3))

# Esempio 3.5 (pag. 181)
x <- c(5:35)
pl <- 5.83 + 0.077*x
mu <- exp(pl)
k <- length(mu)
y <- numeric(k)
set.seed(246135)
for (i in 1:k) y[i] <- rpois(1,mu[i])
rbind(x,y)
fit <- glm(y ~ x, family=poisson("log"))
fit$coefficients

# Figura 3.19d (pag. 171)
# -------------------------------------------------------------------------------
plot(x,y, xlab="temperatura (gradi)", ylab="numero di ciclisti", main="(d)")
# curva vera
lines(x,mu)
# stimata
b <- fit$coefficients
pl <- b[1] + b[2]*x
ph <- exp(pl)
lines(x,ph,lty=3,lwd=2)
# -------------------------------------------------------------------------------

# Esempio 3.6 (pag. 183)
y <- c(1510463,1058684)
n <- c(27616213,32025275)/1000
x <- c(1,0)
fit <- glm(y ~ x + offset(log(n)), family=poisson("log"))
fit$coefficients
exp(fit$coefficients)
inc <- y/n
inc[2]
inc <- y/n
inc[1]
b <- fit$coefficients
exp(b[1] + b[2])
inc[1]/inc[2]
exp(b[2])
exp(confint(fit)[2,])