% Creiamo il vettore delle frequenze assolute e delle probabilità empiriche

f <- matrix(nrow=11)
p <- matrix(nrow=11)
f[1]=7
f[2]=19
f[3]=18
f[4]=25
f[5]=27
f[6]=30
f[7]=25
f[8]=21
f[9]=23
f[10]=12
f[11]=5
p<-f/sum(f)

% vediamo la distribuzione di probabilità (chi vuole può cercare
opzioni grafiche molto migliori):

plot(p)

% idealizziamola con una Poisson:

lambda<-sum((0:10)*p)
p.Poisson <- exp(-lambda)*(lambda^(0:10))/(factorial(0:10))

plot(p.Poisson)

% Creiamo i numeri N[n] di richieste casuali, relative ai giorni n=1,...,n0

n0 <- 212
N.Pois <- rpois(n0,lambda)

%%%% Esercizio: usare la distribuzione originaria p.

% Numero di unità operative:

M <- 1

% Stato X[n] = numero di richieste alla chiusura n-esima:

X <- matrix(nrow=n0)
L <- matrix(nrow=n0)
X[1] <- 0
L[1] <- 0
for (n in 2:n0) {
L[n] <- min(M,X[n-1])
X[n] <- min(X[n-1]-L[n] + N.Pois[n],2*M)
}

% per rendersi conto di come vanno le cose:

plot(X)

% Cambiamo M: ripetere i comandi impostando M<-2 ecc.

% Calcoliamo il guadagno medio giornaliero:

G<-(sum(L)*1000  - M*50000)/212
G

% eseguiamo ora tutto insieme per diversi valori di M, mediando meglio (A anni):

A<-5
n0 <- 212*A
N.Pois <- rpois(n0,lambda)
X <- matrix(nrow=n0)
L <- matrix(nrow=n0)
X[1] <- 0
L[1] <- 0
for (n in 2:n0) {
L[n] <- min(M,X[n-1])
X[n] <- min(X[n-1]-L[n] + N.Pois[n],2*M)
}
G<-(sum(L)*1000  - M*50000*A)/(212*A)
G


% conviene avere 5 unità operative. 

--------------------------------
--------------------------------

% A tempo continuo. Arrivi secondo Poisson:

A<-5
n0 <- 212*A 
T.a <- rweibull(n0,1,1)

hist (T.a)

hist (T.a,50)

% ma lambda è il numero medio di richieste giornaliere, quindi è il parametro
% dell'intertempo esponenziale (nell'unità di misura dei giorni)
% quindi la "scale" della Weibull è 1/lambda:

T.a <- rweibull(n0,1,1/lambda)

hist (T.a,50)

% Tempi di servizio (media 1 giorno):

T.s <- rweibull(n0,1,1)

% Esempio per iniziare: un servitore (M=1)
% Stato X[n] (numero ordinazioni da eseguire) in cui si arriva alla transizione n-esima

X <- matrix(nrow=n0)
X[1] <- 0

% probabilità totali:

p.up <- lambda/(lambda+1)

% p.down <- 1/(lambda+1), ma esplicitamente non serve

for (n in 2:n0) {
X[n] <- X[n-1] + 2*rbinom(1,1,p.up)-1
}

plot(X)

% è una coda esplosiva.
% C'è un errore: le transizioni sono giuste solo per X[n-1]>0.
% Inoltre, i tempi non sono stati utilizzati. A che possono servire?

% Visualizziamo un caso diverso:

lambda <-  0.9
p.up <- lambda/(lambda+1)
for (n in 2:n0) {
X[n] <- X[n-1] + 2*rbinom(1,1,p.up)-1
}

plot(X)


lambda <-  0.99
p.up <- lambda/(lambda+1)
for (n in 2:n0) {
X[n] <- X[n-1] + 2*rbinom(1,1,p.up)-1
}

plot(X)


lambda <-  1
p.up <- lambda/(lambda+1)
for (n in 2:100000) {
X[n] <- X[n-1] + 2*rbinom(1,1,p.up)-1
}

plot(X)