# pag. 4
factorial(52)

(x <- c(1:52))
destra <- x[1:26]
sinistra <- x[27:52]
y <- numeric(0)
for (i in 1:26) y <- c(y,destra[i],sinistra[i])
y

# la scozzata perfetta (pag. 5) -----------------
x <- c(1:52)
cat("Ordinamento iniziale del mazzo\n")
x
for (k in 1:8)
{
     a <- x[1:26]
     b <- x[27:52]
     y <- numeric(0)
     for (i in 1:26) y <- c(y,a[i],b[i])
     x <- y
     cat("Scozzata numero",k,"\n")
     print(x)
}
# -----------------------------------------------

# pag. 5
mazzo <- c(1:52)
mano <- sample(mazzo, size=52, replace=FALSE)
mano

# pag. 7
choose(5,3)
choose(21,6)
choose(52,5)

# modelli di occupazione (pag.12) ---------------
n <- seq(2,20,2)
y <- choose(n,n/2)
rbind(n,y)

n <- seq(2,20,2)
cfav <- choose(n,n/2)
cpos <- 2^n
p <- cfav/cpos
rbind(n,round(p,3))

n <- 1000
choose(n,n/2)/(2^n)

n <- 1030
choose(n,n/2)/(2^n)

n <- 1030
(lnum <- lchoose(n,n/2))
(lden <- n*log(2))
exp(lnum - lden)
# -----------------------------------------------

# estrazione di un campione da una popolazione (pag. 14)
set.seed(123456)
sample(c("1","X","2"),size=13,replace=TRUE)

set.seed(123456)
sample(c("1","X","2"),size=13,replace=TRUE,prob=c(0.5,0.4,0.1))
# -----------------------------------------------

# Il problema della divisione della posta (pag. 17)
set.seed(123456)
a <- 5; b <- 3
repeat {
  x <- sample(c("A","B"),1)
  cat(x)
  ifelse (x == "A", a <- a + 1, b <- b + 1)
  if ((a==6) | (b==6)) {
      cat("\n")
      break
}}    
c(a,b)
    
nrep <- 100000
A <- 0; B <- 0
set.seed(135246)
for (i in 1:nrep) {
 a <- 5; b <- 3
 repeat {
      x <- sample(c(1,2),1)
      ifelse (x == 1, a <- a + 1, b <- b + 1)
      if ((a==6) | (b==6)) break
}}
if(a == 6) A <- A + 100
if(b == 6) B <- B + 100
c(A,B)/100)
c(A,B)/nrep
    
for (i1 in (c("A","B"))){
    for (i2 in (c("A","B"))){
         for (i3 in (c("A","B"))){
             cat(i1,i2,i3,"\n",sep="")
}}}

nrep <- 100000    
A <- 0; B <- 0
set.seed(135246)
for (i in 1:nrep) {
 a <- 50; b <- 30
 repeat {
      x <- sample(c(1,2),1)
      ifelse (x == 1, a <- a + 1, b <- b + 1)
      if ((a==60) | (b==60)) break
}
if(a == 60) A <- A + 100
if(b == 60) B <- B + 100
}
c(A,B)/100
c(A,B)/nrep
    
sum(dbinom(c(10:39),39,0.5))
sum(dbinom(c(30:39),39,0.5))
# -----------------------------------------------
    
# Il problema di Galileo (pag. 21) --------------
for (d1 in 1:6){
	for (d2 in 1:6){
		for (d3 in 1:6){
			cat(d1,d2,d3,"\n")
}}}

for (d1 in 1:6){
	for (d2 in 1:6){
		for (d3 in 1:6){
			S <- d1+ d2 + d3
			if (S==9 | S==10) cat(d1,d2,d3,S,"\n")
}}}
# -----------------------------------------------

# Il problema del cavalier de Méré (pag. 23) ----
cfav <- 0; cpos <- 0
for (d1 in 1:6){
    for (d2 in 1:6){
        for (d3 in 1:6){
            for (d4 in 1:6){
                 cpos <- cpos + 1
                 if (d1==6 | d2==6 | d3==6 |d4==6) cfav <- cfav + 1
}}}}
c(cfav,cpos)
cfav/cpos

1 - (35/36)^24
# -----------------------------------------------

# Il problema del compleanno (pag. 26) ----------
n <- c(2:5,10,15,20,25,30,40,50)
for (i in n) cat(i,pbirthday(i),"\n")

qbirthday(0.99)

# Figura 1.1 (pag. 27)
# -------------------------------------------------------------------------------
# (a)
N <- 50
n <- c(1:N)
p <- numeric(N)
for (i in n) p[i] <- pbirthday(i)
plot(n,p,xlab="numero di persone",ylab="probabilità")
# lines(c(0,23),c(0.5,0.5),lty=3)
# lines(c(23,23),c(0,0.5),lty=3)
# (b)
p <- seq(0.001,0.999,0.001)
n <- numeric(length(p))
k <- 0
for (i in p){
      k <- k + 1
      n[k] <- qbirthday(i)
}
plot(p,n,type="l",ylab="numero di persone",xlab="probabilità")
# -------------------------------------------------------------------------------

# pag. 28
anno <- c(1:365)
set.seed(654321)
nrep <- 100000
n <- 23
x <- sample(anno,nrep*n,replace=TRUE)
X <- matrix(x,nrow=nrep,byrow=TRUE)
dum <- function(y) ifelse(sum(duplicated(y))>0,1,0)
ris <- apply(X,1,dum)
sum(ris)/nrep

pbirthday(50,classes = 365, coincident = 3)
qbirthday(0.9,classes = 365, coincident = 3)
# -----------------------------------------------

# Il modello di Ehrenfest (pag. 31) -------------
Ehr <- function(n,nmax=500000) {
	x <- c(1:n)
	y <- rep("A",n)
	out <- numeric(nmax)
	nest <- 0
	repeat {
		k <- sample(x,size=1)
		nest <- nest + 1
		if (nest > nmax) stop("Numero massimo di passi raggiunto.")
		ifelse((y[k] == "A"), y[k] <- "B", y[k] <- "A")
		chk <- sum(y == "A")
		out[nest] <- chk
		if (chk == n) break
	}
	cat("Numero di estrazioni eseguite:",nest,"\n")
	out <- out[1:nest]
	return(out)
}

set.seed(123); out <- Ehr(4); out 
plot(out,type="l",lty=3)

set.seed(456); out <- Ehr(5); out
set.seed(567); out <- Ehr(10)

# Figura 1.2 (pag. 33)
# -------------------------------------------------------------------------------
# (a)
set.seed(678); out <- Ehr(14)
y <- table(out)/sum(table(out))*100
x <- c(0:14)
plot(x,y, type="h", xlab="Numero di palline nella scatola A", ylab="(%)")
# (b)
set.seed(789); out <- Ehr(20,10000000)
y <- table(out)/sum(table(out))*100
x <- c(0:20)
plot(x,y, type="h", xlab="Numero di palline nella scatola A", ylab="(%)")
# -----------------------------------------------

# pag. 36
p <- 3/5
n <- 25551
x <- rbinom(1,n,p)
x/n

# Figura (pag. 50)
n <- c(1:25)
den <- 0.001+0.999*0.5^n
p <- 0.001/den
plot(n,p)
cbind(n[13:14],p[13:14])

# probabilità, odds e scommesse (pag. 56)
p <- 0.4
odds <- p/(1-p)
quota <- 1/odds
x <- c(quota,-1)
nrep <- 100000
set.seed(123456)
a <- sample(x,nrep,prob=c(p,(1-p)),replace=TRUE)
(tbl <- table(a))
1.5*tbl[2] - tbl[1]
a <- sample(x,nrep,prob=c(p,(1-p)),replace=TRUE)
(tbl <- table(a))
1.5*tbl[2] - tbl[1]

# Figura 1.7 (pag. 57)
# -------------------------------------------------------------------------------
# (a)
curve(x/(1-x),0,0.9,xlab="Probabilità",ylab="Odds",
cex.lab=1.25,cex.axis=1.25)
abline(v=0.5,h=1,lty=3)
# (b)
curve(log(x/(1-x)),0.001,0.999,xlab="Probabilità",
ylab="Logit",cex.lab=1.25,cex.axis=1.25)
abline(v=0.5,h=0,lty=3)
# -------------------------------------------------------------------------------

# pag. 58
(p <- seq(0,1,by=0.1))
round(p/(1-p),3)

# Figura 1.8 (pag. 60)
# -------------------------------------------------------------------------------
# (a)
plot(c(0,100),c(-10,10),type="n",xlab="livello di stress psico-fisico",ylab="logit")
abline(-7,0.14)
# (b)
curve(exp(-7+0.14*x)/(1+exp(-7+0.14*x)),0,100,xlab="livello di stress psico-fisico",ylab="probabilità")
# -------------------------------------------------------------------------------

# pag. 62
x <- c(11.5,7,1.22)
(quote <- x - 1)
odds <- 1/quote
round(odds,3)
prob <- odds/(1 + odds)
round(prob,3)
sum(prob)
(K <- 1/sum(prob))
round(prob*K,3) 
sum(prob*K)

# pag. 64
x <- c(11.5,7.00,1.22)
quote <- x - 1
odds <- 1/quote
prob <- odds/(1 + odds)
K <- 1/sum(prob)
p <- prob*K
(X <- cbind(quote*100,p,quote*100*p,-100,(1-p),-100*(1-p)))
Y <- X[,c(3,6)]
rowSums(Y)

# pag. 65
x <- c(11.5,7.00,1.22)
quote <- x - 1
odds <- 1/quote
prob <- odds/(1 + odds)
(K <- 1/sum(prob))
K - 1

# pag. 66
x <- c(1.75,3,12,13,17,21,51,251,501,501,501,751,751,751,751,751,1001,1001)
quote <- x - 1
quote
odds <- 1/quote
round(odds,3)
prob <- odds/(1 + odds)
round(prob,3)
sum(prob)
round(prob/sum(prob),3)