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

# pag. 69
f <- function(x) 2*x^2 + 3*x + 4
x <- c(0,1)
f(x)

# pag. 73
x <- c(-1,124,249,499)
f <- c(0.997,0.001,0.001,0.001)
sum(x*f)

# pag. 74
x <- c(-1,1000/7-1,2000/7-1,4000/7-1)
f <- c(0.997,0.001,0.001,0.001)
sum(x*f)

x <- c(1,-1)
f <- c(18/37,19/37)
sum(x*f)

x <- c(35,-1)
f <- c(1/37,36/37)
sum(x*f)

# pag. 75
x <- c(1,-1)
f <- c(18/37,19/37)
m <- sum(x*f)
sum((x-m)^2*f)

# pag. 76
x <- c(35,-1)
f <- c(1/37,36/37)
m <- sum(x*f)
sum((x-m)^2*f)

# pag. 77
k <- c(seq(1.5,5,0.5),6)
p <- 1 - 1/k^2
rbind(k,round(p,3))

# Uniforme discreta (pag. 79) -------------------
x <- c(1:6)
f <- rep(1/6,6)
(mu <- sum(x*f))
(sigma2 <- sum((x-mu)^2*f))

x <- c(0:9)
f <- rep(1/10,10)
(mu <- sum(x*f))
(sigma2 <- sum((x-mu)^2*f))

x <- c(1:6)
f <- rep(1/6,6)
F <- cumsum(f)
data.frame(x,f,F)

x <- c(0:9)
f <- rep(1/10,10)
F <- cumsum(f)
data.frame(x,f,F)

# Figura 2.1 (pag. 81)
# -------------------------------------------------------------------------------
pippo <- function(tmpx,tmpy,point=TRUE){
  n   <- nrow(tmpx)
  nm1 <- n - 1
  for (i in 1:n) lines(tmpx[i,],tmpy[i,])
  if (point) for (i in 1:nm1) points(tmpx[i+1,1],tmpy[i+1,2],pch=19)
  for (i in 1:nm1){
    xx <- c(tmpx[i,2],tmpx[i+1,1])
    yy <- c(tmpy[i,2],tmpy[i+1,1])
    lines(xx,yy,lty=3)
  }
}

par(mfrow=c(2,2))

# (a)
x <- c(1:6)
f <- rep(1/6,6)
F <- cumsum(f)
plot(x,f,type="h",ylab="Probabilità",main="Lancio di un dado (f)",xlab="Punteggio")
points(x,f,pch=19)
# (b)
plot(c(0,7),c(0,1),type="n",ylab="Probabilità",main="Lancio di un dado (F)",xlab="Punteggio")
tmpx <- cbind(c(0:6),c(1:7))
tmp <- c(0:6)/6
tmpy <- cbind(tmp,tmp)
pippo(tmpx,tmpy)
# (c)
x <- c(0:9)
f <- rep(1/10,10)
F <- cumsum(f)
plot(x,f,type="h",ylab="Probabilità",main="Estrazione di una pallina (f)",xlab="Punteggio")
points(x,f,pch=19)
# (d)
plot(c(-1,10),c(0,1),type="n",ylab="Probabilità",main="Estrazione di una pallina (F)",xlab="Punteggio")
tmpx <- cbind(c(-1:9),c(0:10))
tmp <- c(0:10)/10
tmpy <- cbind(tmp,tmp)
pippo(tmpx,tmpy)

par(mfrow=c(1,1))
# -------------------------------------------------------------------------------

# Triangolare discreta (pag. 83) -------------------
x <- c(1:6)
(tbl <- outer(x,x,"+"))
table(tbl)
f <- prop.table(table(tbl))
x <- c(2:12)
data.frame(x,f=as.numeric(f),F=cumsum(f))
(mu <- sum(x*f))
(sigma2 <- sum((x-mu)^2*f))

# Figura 2.2 (pag. 83)
# -------------------------------------------------------------------------------
x <- c(1:6)
(tbl <- outer(x,x,"+"))
table(tbl)
f <- prop.table(table(tbl))
x <- c(2:12)
data.frame(x,f=as.numeric(f),F=cumsum(f))
# (a)
plot(x,f,type="h",ylab="Probabilità",main="f",xlab="Punteggio totale")
points(x,f,pch=19)
# (b)
plot(c(1,13),c(0,1),type="n",ylab="Probabilità",main="F",xlab="Punteggio totale")
tmpx <- cbind(c(1:12),c(2:13))
F <- cumsum(f)
tmp <- c(0,F)
tmpy <- cbind(tmp,tmp)
pippo(tmpx,tmpy)
# -------------------------------------------------------------------------------

# Figura 2.3 (pag. 85)
# -------------------------------------------------------------------------------
dado <- c(1:6)
nrep <- 100000
ndad <- 12
set.seed(123456)
y <- sample(dado,ndad*nrep,replace=TRUE)
Y <- matrix(y,ncol=ndad)
x <- rowSums(Y)
f <- prop.table(table(x))
x <- dimnames(f); x <- as.numeric(unlist(x))
(mu <- sum(x*f))
(sigma2 <- sum((x-mu)^2*f))
data.frame(x,f=as.numeric(f),F=cumsum(f))
# (a)
plot(x,f,type="h",ylab="Probabilità",main="f",xlab="Punteggio totale")
# (b)
plot(c(15,70),c(0,1),type="n",ylab="Probabilità",main="F",xlab="Punteggio totale")
tmpx <- cbind(c(15:67),c(16:68))
F <- cumsum(f)
tmp <- c(0,F,1); tmpy <- cbind(tmp,tmp)
pippo(tmpx,tmpy,point=FALSE)
# -------------------------------------------------------------------------------

# Geometrica (pag. 87) --------------------------
x <- c(0:5)
p <- 0.6
f <- (1-p)^x * p
f
exp(diff(log(f)))

# Figura 2.4 (pag. 89)
# -------------------------------------------------------------------------------
x <- c(0:10)
p <- 0.8
f <- p * (1-p)^x
# (a)
plot(x,f,type="h",ylab="Probabilità",main="p = 0.8",xlab="Estrazione")
points(x,f,pch=19)
lines(x,f,lty=3)
# (b)
x <- c(0:10)
p <- 0.2
f <- p * (1-p)^x
plot(x,f,type="h",ylab="Probabilità",main="p = 0.2",xlab="Estrazione")
points(x,f,pch=19)
lines(x,f,lty=3)
# -------------------------------------------------------------------------------

# Binomiale (pag. 92) ---------------------------
dbinom(c(0:3),3,0.6)
Binomiale(n=3, p=0.6)

# Figura 2.5 (pag. 93)
# -------------------------------------------------------------------------------
x <- c(0:10)
par(mfrow=(c(3,3)))
f <- dbinom(x,10,0.1); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.1"); points(x,f,pch=19)
f <- dbinom(x,10,0.2); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.2"); points(x,f,pch=19)
f <- dbinom(x,10,0.3); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.3"); points(x,f,pch=19)
f <- dbinom(x,10,0.4); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.4"); points(x,f,pch=19)
f <- dbinom(x,10,0.5); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.5"); points(x,f,pch=19)
f <- dbinom(x,10,0.6); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.6"); points(x,f,pch=19)
f <- dbinom(x,10,0.7); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.7"); points(x,f,pch=19)
f <- dbinom(x,10,0.8); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.8"); points(x,f,pch=19)
f <- dbinom(x,10,0.9); plot(x,f,type="h",ylab="f(x)",main="n=10, p=0.9"); points(x,f,pch=19)
par(mfrow=(c(1,1)))
# -------------------------------------------------------------------------------

# Figura 2.6 (pag. 94)
# -------------------------------------------------------------------------------
# (a)
n <- 1000
x <- c(0:n)
x <- c(430:570)
f <- dbinom(x,n,0.5); plot(x,f,type="h",main="n=1000, p=0.5",ylab="f(x)")
# (b)
n <- 1000
x <- c(0:n)
x <- c(60:140)
f <- dbinom(x,n,0.1); plot(x,f,type="h",main="n=1000, p=0.1",ylab="f(x)")
# -------------------------------------------------------------------------------

# passeggiata casuale
# Figura 2.7 (pag. 95)
# -------------------------------------------------------------------------------
# (a)
plot(c(1,10),c(-4,3),type="n",xlab="Numero di passi",ylab="Distanza dall'origine")
abline(h=0)
n <- 10; x <- c(1:n)
set.seed(123)
y <- sample(c(-1,1),n,replace=TRUE); y
yy <- cumsum(y); lines(x,yy)
set.seed(456)
y <- sample(c(-1,1),n,replace=TRUE); y
yy <- cumsum(y); lines(x,yy,lty=2)
# (b)
plot(c(1,100),c(-10,15),type="n",xlab="Numero di passi",ylab="Distanza dall'origine")
abline(h=0)
n <- 100; x <- c(1:n)
set.seed(123)
y <- sample(c(-1,1),n,replace=TRUE); sum(y)
yy <- cumsum(y); lines(x,yy)
set.seed(456)
y <- sample(c(-1,1),n,replace=TRUE); sum(y)
yy <- cumsum(y); lines(x,yy,lty=2)
# -------------------------------------------------------------------------------

# pag. 97
set.seed(654321)
n <- c(seq(10,100,10),seq(200,1000,100))
out <- matrix(nrow=length(n),ncol=2)
k <- 0
for (i in n){
      x <- rbinom(100000,i,0.5)
      d <- 2*x - i
      k <- k + 1
      out[k,] <- c(mean(d^2),mean(abs(d)))
}
r1 <- round(out[1:10,1],0)
r2 <- round(out[1:10,2],2)
rbind(n=n[1:10],r1,r2)
fit <- lm(out[,2] ~ sqrt(n))
(b <- fit$coefficients)

# Figura 2.8 (pag. 98)
# -------------------------------------------------------------------------------
# plot(n,out[,1],xlab="Numero di passi",ylab="Distanza(^2) media dall'origine")
# (a)
plot(n,out[,2],xlab="Numero di passi",ylab="Media del valore assoluto della distanza")
# (b)
plot(sqrt(n),out[,2],xlab="Radice quadrata del numero di passi",ylab="Media del valore assoluto della distanza")
# -------------------------------------------------------------------------------

# pag. 98
n <- 10
x <- c(0:n)
d <- 2*x - n
f <- dbinom(x,n,0.5)
rbind(x,d,f=round(f,3))
sum(abs(d)*f)

# pag. 99
n <- c(1:1000)
out <- numeric(length(n))
k <- 0
for (i in n){
      x <- c(0:i)
      f <- dbinom(x,i,0.5)
      d <- 2*x - i
      k <- k + 1
      out[k] <- sum(abs(d)*f)
}
fit <- lm(out ~ sqrt(n))
(b <- fit$coefficients)

# Figura 2.9 (pag. 100)
# -------------------------------------------------------------------------------
nn <- c(seq(10,100,10),seq(200,1000,100))
# (a)
r <- fit$residuals[nn+1]
plot(nn+1,r,xlab="Numero di passi",ylab="Residui")
# (b)
r <- fit$residuals[nn]
plot(nn,r,xlab="Numero di passi",ylab="Residui")
# -------------------------------------------------------------------------------

# pag. 100
n <- 1000000
x <- c(0:n)
f <- dbinom(x,n,0.5)
d <- 2*x - n
(m <- sum(abs(d)*f))
2*n/m^2

# fluttuazioni
sum(dbinom(c(490:510),1000,0.5))
1 - sum(dbinom(c(4900:5100),10000,0.5))

# Figura 2.10 (pag. 102)
# -------------------------------------------------------------------------------
pflb <- function(n,fl){
      m <- n/2
      s <- sqrt(n/4)
      me <- n*fl
      ei <- round(m-me,0)
      es <- round(m+me,0)
      int <- c(ei:es)
      p <- sum(dbinom(int,n,0.5))
      p
}
# (a)
n <- seq(100,20000,100)
p <- numeric(length(n))
k <- 0
for (i in n){
      k <- k + 1
      p[k] <- pflb(i,0.01)
}
plot(n,1-p,ylab="probabilità")
# (b)
n <- seq(10000,2000000,10000)
p <- numeric(length(n))
k <- 0
for (i in n){
      k <- k + 1
      p[k] <- pflb(i,0.001)
}
xlb <- bquote( ~ n %*% 10^5)
plot(n,1-p,axes=FALSE,xlab=xlb,ylab="probabilità")
axis(1,at=seq(0,20,5)*100000,labels=seq(0,20,5))
axis(2,at=seq(0,1,0.2))
box()
# -------------------------------------------------------------------------------

# Ipergeometrica (pag. 104) ---------------------
dhyper(c(0:3),6,4,3)
# -----------------------------------------------

# Poisson (pag. 106) ----------------------------
round(dpois(c(0:5),0.5),6)
Poisson(mu=0.5)

# Figura 2.11 (pag. 107)
# -------------------------------------------------------------------------------
par(mfrow=(c(2,2)))
x <- c(0: 5); f <- dpois(x,0.5); plot(x,f,type="h",ylab="f(x)",main=bquote( ~ mu == 0.5)); points(x,f,pch=19,cex=0.75)
x <- c(0: 7); f <- dpois(x,  1); plot(x,f,type="h",ylab="f(x)",main=bquote( ~ mu == 1));   points(x,f,pch=19,cex=0.75)
x <- c(0:10); f <- dpois(x,  2); plot(x,f,type="h",ylab="f(x)",main=bquote( ~ mu == 2));   points(x,f,pch=19,cex=0.75)
x <- c(0:15); f <- dpois(x,  4); plot(x,f,type="h",ylab="f(x)",main=bquote( ~ mu == 4));   points(x,f,pch=19,cex=0.75)
par(mfrow=(c(1,1)))
# -------------------------------------------------------------------------------

# Figura 2.12 (pag. 108)
# -------------------------------------------------------------------------------
# (a)
m <- 10
x <- c(0:30)
f <- dpois(x,m); plot(x,f,type="h",main=bquote( ~ mu == 10),ylab="f(x)")
# (b)
m <- 100
x <- c(60:140)
f <- dpois(x,m); plot(x,f,type="h",main=bquote( ~ mu == 100),ylab="f(x)")
# -------------------------------------------------------------------------------

# Poisson troncata (pag. 109)
set.seed(654321)
x <- rpois(100000,1)
c(mean(x),var(x))

fo <- table(x)
f <- dpois(c(0:7),1)
fa <- round(f*100000,0)
rbind(fo,fa)
sum(fa)

y <- x[x > 0]
length(y)
c(mean(y),var(y))

g <- function(l,m) l/(1-exp(-l)) - m
(m <- mean(y))
l <- seq(0.8,1.6,0.1)
cbind(l,g(l,m))

l <- seq(1,1.1,0.01)
cbind(l,g(l,m))

g <- function(l,m) l/(1-exp(-l)) - m
(m <- mean(y))
uniroot(g,c(0.5,1.6),m)

dpois(0,0.5)
set.seed(654321)
x <- rpois(100000,0.5)
c(mean(x),var(x))
sum(x==0)

y <- x[x > 0]
c(mean(y),var(y))
uniroot(g,c(0.1,1.3),mean(y))

dpois(0,3)
set.seed(654321)
x <- rpois(100000,3)
c(mean(x),var(x))
sum(x==0)
y <- x[x > 0]
c(mean(y),var(y))
uniroot(g,c(1.5,3.5),mean(y))
# -----------------------------------------------

# Binomiale negativa (pag. 113) -----------------
f <- function(x,p) choose(x,2) * p^3 * (1-p)^(x-2)
x <- c(2:4)
f(x,0.6)

bn <- function(x,k,p) choose(x,(k-1)) * p^k * (1-p)^(x-k+1)
x <- c(0:4)
cbind(x,"f(x)"=bn(x,k=3,p=0.6))

x <- c(2,3,4); bn(x,k=3,p=0.6)
x <- c(2,3,4); dnbinom(x,3,0.6)

x <- c(2,3,4); bn(x,k=3,p=0.6)
x <- c(0,1,2); dnbinom(x,3,0.6)

x <- c(2:999)
f <- bn(x,k=3,p=0.6)
sum(f)
(m <- sum(x*f))
sum((x-m)^2*f)

y <- c(0:999)
f <- dnbinom(y,size=3,p=0.6)
sum(f)
(m <- sum(y*f))
sum((y-m)^2*f)
# -----------------------------------------------