Προσομοίωση Σειρά 3 - Εφαρμογές

Άσκηση 1

Έστω \(X\) τυχαία μεταβλητή που ακολουθεί εκθετική κατανομή μέ παράμετρο \(\lambda\) και μια εκτιμήτρια του \(\lambda\) που δίνεται από \[ \hat{\lambda} = \frac{1}{\bar{X}_n} \] όπου \(X_1,\ldots, X_n\) τυχαίο δείγμα από την κατανομή και \(\bar{X}_n\) ο δειγματικός μέσος. Να δημιουργηθεί μια συνάρτηση που εκτιμά μέσω προσομοίωσης τη μεροληψία της εκτιμήτριας ως συνάρτηση των \(\lambda, n\) χρησιμοποιώντας \(N\) επαναλήψεις.

Απάντηση

Η συνάρτηση που θα σχεδιάσουμε παίρνει ως ορίσματα τα \((\lambda, n, N)\), και πρσομοιώνει \(N\) τυχαία δείγματα μεγέθους \(n\) από την εκθετική κατανομή με παράμετρο \(\lambda\). Για κάθε δείγμα υπολογίζει την εκτίμηση \(\lambda\) και την απόκλιση από την πραγματική τιμή. Η εκτίμηση της μεροληψίας προκύπτει από τη μέση τιμή των αποκλίσεων πάνω στις \(N\) επαναλήψεις. Επίσης υπολογίζει το μέσο τετραγωνικό σφάλμα \(E[(\hat{\lambda}-\lambda)^2]\) από τη μέση τιμή των τετραγωνικών αποκλίσεων

expbias = function(lambda, n, N)
{
  A=rep(0,N)
  for (j in 1:N) # N simulation replications
  {
    x=rexp(n, lambda)
    lambdaest=1/mean(x)
    A[j]=lambdaest-lambda
  }
  biasest=mean(A)
  mseest=mean(A^2)
  return(c(biasest, mseest))
}

Για εφαρμογή καλούμε τη συνάρτηση για \(\lambda=1, N=1000\) και εκτιμούμε τη μεροληψία και το μέσο τετραγωνικό σφάλμα συναρτήσει του \(n\), για τιμές του \(n=10,20,\ldots,100\):

nval=seq(10,100,10)
nn= length(nval)
A=matrix(rep(0, 2*nn), nn, 2) # matrix A contains the bias in the first column and MSE in the second
for (i in 1:nn)
  A[i,]=expbias(1, nval[i],1000)
A

##             [,1]       [,2]
##  [1,] 0.12552983 0.16472979
##  [2,] 0.04691985 0.06877246
##  [3,] 0.03410658 0.03954108
##  [4,] 0.01394127 0.02676715
##  [5,] 0.02433589 0.02191314
##  [6,] 0.01928851 0.02033749
##  [7,] 0.01619736 0.01576417
##  [8,] 0.02156750 0.01347362
##  [9,] 0.01101891 0.01111367
## [10,] 0.01488097 0.01092106

Τώρα δημιουργούμε γραφήματα ως προς \(n\)

par(mfrow=c(1,2))
plot(A[,1]~nval, main="Bias", xlab="n", ylab="Bias", type="b")
plot(A[,2]~nval, main="MSE", xlab="n", ylab="MSE", type="b")

Παρατηρούμε ότι η εκτιμήτρια που εξετάζουμε συστηματικά υπερεκτιμά την πραγματική τιμή της παραμέτρου \(\lambda\), ενώ όταν αυξάνεται το μέγεθος δείγματος \(n\) τόσο η μεροληψία όσο και το μέσο τετραγωνικό σφάλμα τείνουν στο μηδέν.

Άσκηση 2

Έστω \(X\) τυχαία μεταβλητή που ακολουθεί εκθετική κατανομή μέ παράμετρο \(\lambda\) και μια εκτιμήτρια του \(\lambda\) που δίνεται από \[ \hat{\lambda} = \frac{1}{\bar{X}_n} \] όπου \(X_1,\ldots, X_n\) τυχαίο δείγμα από την κατανομή και \(\bar{X}_n\) ο δειγματικός μέσος. Να δημιουργηθεί μια συνάρτηση που δημιουργεί ένα δείγμα μεγέθους \(n\) και εκτιμά τη μεροληψία και το μέσο τετραγωνικό σφάλμα της εκτιμήτριας παίρνοντας \(N\) bootstrap δείγματα από την εμπειρική κατανομή του αρχικού δείγματος.

Απάντηση Η συνάρτηση που θα σχεδιάσουμε παίρνει ως ορίσματα τα \((\lambda, n, N)\), και πρσομοιώνει ένα τυχαίο δείγμα μεγέθους \(n\), έστω \(x=(x_1, \ldots, x_n)\) από την εκθετική κατανομή με παράμετρο \(\lambda\). Αν \(F_e\) είναι η εμπειρική κατανομή του δείγματος, τότε η προσέγγιση του \(\theta = 1/Ε(Χ)\) μέσω της \(F_e\) είναι η \(\theta(F_e) = 1/E_{F_e}(X) = 1/\bar{x}_n\).

Για την εκτίμηση μέσω bootstrap, παίρνουμε \(N\) προσομοιωμένα δείγματα μεγέθους \(n\) από την \(F_e\). Όπως έχουμε δει, η \(F_e\) είναι η ομοιόμορφη διακριτή στο σύνολο \(\{x_1, \ldots, x_n \}\). Για το \(j\) δείγμα bootstrap, έστω \(X_j=(X_{1j}, \ldots,X_{nj}\), υπολογίζουμε την ποσότητα \(1/\bar{X_j}\), την απόκλιση και την τετραγωνική απόκλιση από το \(\theta(F_e)\): \[ B_j=1/\bar{X_j} - \theta(F_e) \] και \[ S_j=(1/\bar{X_j} - \theta(F_e))^2. \] Η εκτίμηση της μεροληψίας είναι \[ \hat{B}=\frac{1}{N}\sum_{j=1}^N B_j \] και του μέσου τετραγωνικού σφάλματος \[ \hat{MSE}=\frac{1}{N}\sum_{j=1}^N S_j. \] Η επόμενη συνάρτηση υλοποιεί τα παραπάνω:

expbootstrap = function(lambda, n, N)
{
  x=rexp(n, lambda)
  thetafe=1/mean(x)
  bootstrapestimate=rep(0,N)
  for (j in 1:N)
  {
    X=sample(x,n,replace=T)
    bootstrapestimate[j]=1/mean(X)
  }
  biasest=mean(bootstrapestimate-thetafe)
  mseest=mean((bootstrapestimate-thetafe)^2)
  return(c(biasest, mseest))
}

Για εφαρμογή, καλούμε τη συνάρτηση με \(\lambda=1, n=40, N=1000\):

bootstrapest=expbootstrap(1,40,1000)
bootstrapest

## [1] 0.02465174 0.02112729

Άσκηση 3

θεωρούμε μια ουρά αναμονής \(GI/G/m\) με \(m\) σταθμούς εξυπηρέτησης και τα εξής στατιστικά χαρακτηριστικά: Οι αφίξεις των πελατών συμβαίνουν σύμφωνα με μια ανανεωτική διαδικασία, της οποίας οι ενδιάμεσοι χρόνοι ακολουθούν κατανομή Gamma\((a,l)\). Οι χρόνοι εξυπηρέτησης ακολουθούν ομοιόμορφη κατανομή στο διάστημα \((b1, b2)\). Η σειρά εξυπηρέτησης είναι First-Come-First-Served και το σύστημα έχει χωρητικότητα \(K\) που μπορεί να είναι και άπειρη. Αν \(K<\infty\), τότε όσοι πελάτες φτάνουν ενώ το σύστημα είναι πλήρες απορρίπτονται. Κάθε πελάτης που εξυπηρετείται αφήνει στο σύστημα ένα έσοδο ίσο με \(R\). Το κόστος ανά μονάδα χρόνου για τη λειτουργία του συστήματος είναι ίσο με \(Q(\mu) = c m \mu^2\) (ντετερμινιστικό).

Υποθέτουμε ότι το σύστημα λειτουργεί στο χρονικό διάστημα \([0,T]\), με αρχική κατάσταση μηδενικού αριθμού πελατών. Ζητείται να εκτιμηθούν μέσω προσομοίωσης τα παρακάτω μεγέθη:

1. Ο μέσος αριθμός πελατών που βρίσκονται στο σύστημα τη χρονική στιγμή \(t\), για \(n\) ισοκατανεμημένες τιμές του \(t \in [0,T]\), δηλαδή για \(t= jT/n, j=1, 2,\ldots,n\).

2. Ο μέσος χρόνος παραμονής ενός πελάτη στο σύστημα για όλους τους πελάτες που έχουν εξυπηρετηθεί μέχρι και τη χρονική στιγμή \(T\).

3. Ο μέσος αριθμός πελατών που έχουν απορριφθεί στο διάστημα \([0,T]\).

4. Το μέσο καθαρό κέρδος του συστήματος στο διάστημα \([0,T]\).

Απάντηση Για την προσομοίωση ενός σεναρίου από την παραπάνω ουρά αναμονής δημιουργούμε τη συνάρτηση ggmsimul που φαίνεται παρακάτω. Για το σχεδιασμό της προσομοίωσης σημειώνουμε τα εξής. Επειδή η ουρά έχει γενικά περισσότερους από έναν σταθμούς εξυπηρέτησης, η σειρά αναχώρησης των πελατών δεν είναι απαραίτητα η ίδια με τη σειρά άφιξης. Επομένως για να υπολογίσουμε το χρόνο παραμονής ενός πελάτ στο σύστημα, πρέπει να παρακολουθούμε τον πελάτη σε όλη τη διάρκεια του χρόνου μεταξύ άφιξης και αναχώρησης. Για το λόγο αυτό κάθε πελάτης που φτάνει στο σύστημα παίρνει ένα αριθμό ταυτοποίησης που είναι ο αριθμός άφιξής του (ο πρώτος 1 ο δεύτερος 2 κλπ). Η κατάσταση του συστήματος περιλαμβάνει τον αριθμό πελατών που βρίσκονται στο σύστημα (n), ένα διάνυσμα μεγέθους \(m\) (servers) που περιέχει τον αριθμό του πελάτη που εξυπηρετείται σε κάθε σταθμό (0 αν ο σταθμός είναι άδειος) και ένα διάνυσμα μεταβλητού μεγέθους που περιέχει τον αριθμό κάθε πελάτη που βρίσκεται στην ουρά αναμονής σύμφωνα με τη σειρά άφιξης. Όταν συμβεί μια άφιξη αυτή απορρίπτεται αν \(n=K\), διαφορετικά ο πελάτης είτε μπαίνει τελευταίος στην ουρά αν όλοι οι σταθμοί είναι κατειλημμένοι είτε αρχίζει εξυπηρέτηση στον πρώτο ελεύθερο σταθμό. Όταν συμβεί μια αναχώρηση ο σταθμός είτε μπαίνει σε κατάσταση αδράνειας (0) αν ουρά είναι άδεια, είτε αρχίζει να εξυπηρετεί τον πρώτο πελάτη που περιμένει στην ουρά. Σε κάθε περίπτωση μετά από μια άφιξη προσομοιώνεται ο χρόνος επόμενης άφιξης και κατά την έναρξη της εξυπηρέτησης σε ένα σταθμό προσομοιώνεται ο χρόνος αναχώρησης του αντίστοιχου πελάτη. Oι χρόνοι των επόμενων γεγονότων παριστάνονται αντίστοιχα από τη μεταβλητή tnexarr και το διάνυσμα tnextdep μεγέθους \(m\).

Για τη συλλογή των στατιστικών κατά τη διάρκεια της προσομοίωσης ορίζουμε τις εξής μεταβλητές. Οι μεταβλητές narrivals, nserved και nrejected περιέχουν τον αριθμό των πελατών που έχουν φτάσει στο σύστημα, έχουν αναχωρήσει αφού εξυπηρετηθούν και έχουν απορριφθεί λόγω πληρότητας του συστήματος. Ο πίνακας custdata περιέχει το χρόνο άφιξης και το χρόνο αναχώρησης κάθε πελάτη όπως επίσης και μια 0/1 μεταβλητή rej που δείχνει αν ο πελάτης απορρίφθηκε ή όχι. Αν ο πελάτης απορριφθεί οι δύο χρόνοι ταυτίζονται. Στο τέλος της προσομοίωσης οι χρόνοι αναχώρησης των πελατών που βρίσκονται ακόμη στο σύστημα είναι μη ορισμένοι (NA). Τέλος το διάνυσμα systemlength περιέχει τον αριθμό πελατών που βρίσκονται στο σύστημα τις χρονικές στιγμές \(t_j, j=1,\ldots, n\), όπως έχουν οριστεί παραπάνω. Το διάνυσμα αυτό ενημερώνεται ως εξής. Σε κάθε βήμα της προσομοίωσης καθορίζονται οι στιγμές \(t_j\) που περιέχονται στο διάστημα μεταξύ του προηγούμενου και του τρέχοντος γεγονότος. Επειδή κατά τη διάρκεια αυτού του διαστήματος ο αριθμός πελατών στο σύστημα έχει παραμείνει σταθερός, όλες οι συνιστώσες του διανύσματος systemlength που αντιστοιχούν σε αυτές τις χρονικές στιγμές παίρνουν τιμή ίση με τον αριθμό πελατών \(n\) σε αυτό το διάστημα. Τέλος για τον υπολογισμό του συνολικού κέρδους χρησιμοοποείται ο τύπος \(TotalProfit = R\cdot nserved - T c m \left ( \frac{2}{b_1+b_2} \right )^2\). Ο μέσος ρυθμός εξυπηρέτησης είναι ίσος με \(\mu=1/E(S)\), όπου \(E(S)=(b_1+b_2)/2\) είναι ο μέσος χρόνος εξυπηρέτησης.

Η συνάρτηση επιστρέφει ως αποτέλεσμα μια λίστα που περιέχει τον πίνακα custdata, το διάνυσμα systemlength και τις μεταβλητές narrivals, nserved, nrejectes, totalprofit.

#ggm simulation
ggmsimul=function(m, a, l, b1, b2, K, T, c, R, n)
{
  # GI/G/m queue simulation
  # m=number of servers
  # Interarrival times Ta~Gamma(a,l)
  # Service times S~U(b1, b2)
  # System capacity = K
  # Service cost /unit time = c m^2, m=1/E(S)
  # Simulation length = T
  # Record queue length at times tj = jT/n, j=1,2,...,n
  #
  #
  # Event clocks (next arrival, next departure from each server, next event, system clock)
  tnextarr=rgamma(1, a, l)
  tnextdep=rep(Inf,m)
  tnextevent=tnextarr
  t=0
  # Type of next event : 0=arrival, j=departure from server j
  nextevent=0
  custdeparting=NULL
  #
  # System state 
  # List : n=number of customers in the system
  # q : vector of customer id's in the queue in FCFS order
  # servers : vector of customer id's being served at each server (0=server idle)
  systemstate=list(n=0, queue=NULL, servers=rep(0,m))
  #
  # Simulation Statistics
  #
  # Customer Data Matrix
  # For each customer records : Customer id, arrival time, departure time, rejected (0/1) 
  # If rejected, arrival time=departure time
  custdata=NULL
  # Vector of system lengths at times tj=jT/n, j=1,...,n
  systemlength=rep(NA,n)
  # Total number of arrivals
  narrivals=0
  # Total number of served customers
  nserved=0
  # Total number of rejected customers
  nrejected=0
  #
  while (tnextevent<=T)
  {
    # Advance clock
    t0=t
    t=tnextevent
    # 
    # Update system length vector for tj in the interval [t0, t)
    # j = floor(n*t0/T)+1,..., floor(n*t/T)
    j0=floor(n*t0/T)+1
    j1=floor(n*t/T)
    if (j0 <= j1)
      for (j in j0:j1) 
        systemlength[j]=systemstate$n
    # 
    if (nextevent==0) #simulate arrival
    {
      # Record keeping for arrival 
      narrivals=narrivals+1
      tnextarr=t+rgamma(1,a,l)
      custid=narrivals
      if (systemstate$n==K) #arrival rejected
      {
        custdata=rbind(custdata, c(custid, t,t, 1 ))
        nrejected=nrejected+1
      }
        else #arrival joins system
        {
          custdata=rbind(custdata, c(custid, t,NA, 0 ))
          systemstate$n=systemstate$n+1
          if (systemstate$n<=m) # at least one server is idle (previous n was <=m-1)
          {
            firstidleserver=which.max(systemstate$servers==0)
            systemstate$servers[firstidleserver]=custid
            tnextdep[firstidleserver]=t+runif(1,b1,b2)
          }
          else # arrival joins queue
          {
            systemstate$queue=c(systemstate$queue, custid)
          }
        }
    } # simulate arrival
    else #simulate departure
    {
      systemstate$n=systemstate$n-1
      custdeparting=systemstate$servers[nextevent]
      custdata[custdeparting,3]=t
      nserved=nserved+1
      if (length(systemstate$queue)==0) #queue is empty server becomes idle
      {
        systemstate$servers[nextevent]=0
        tnextdep[nextevent]=Inf
      }
      else #first customer in queue enters the server
      {
        systemstate$servers[nextevent]=systemstate$queue[1]
        tnextdep[nextevent]=t+runif(1,b1,b2) 
        if (length(systemstate$queue)==1) #queue becomes empty
          systemstate$queue=NULL
        else # first customer is removed from queue
          systemstate$queue=systemstate$queue[2:length(systemstate$queue)]
      }
    } # simulate departure
    #
    #print(c(t, nextevent,narrivals, nserved, custdeparting))
    #print(c(t, systemstate$queue))
    #print(c(t, systemstate$servers))
    # Determine if next event is an arrival (0) or departure from some server j
    # and time of next event
    if (tnextarr<min(tnextdep)) 
    {
      nextevent=0
      tnextevent=tnextarr
    } else 
    {
      nextevent=which.min(tnextdep)
      tnextevent=min(tnextdep)
    }
  }  #while
  # Complete statistics collection for last interval [t,T]
  t0=t
  t=T
  j0=floor(n*t0/T)+1
  j1=n
  if (j0 <= j1)
    for (j in j0:j1) 
      systemlength[j]=systemstate$n
  totalprofit=R*nserved - T*c*m*(2/(b1+b2))^2
  simoutput=list(custdata=custdata, narrivals=narrivals, nserved=nserved, nrejected=nrejected, systemlength=systemlength, totalprofit=totalprofit)
  return(simoutput)
}
#

Η παραπάνω συνάρτηση προσομοιώνει μια επανάληψη. Για την εκτίμηση των στατιστικών μεγεθών του συστήματος, δημιουργούμε την παρακάτω συνάρτηση που προσομοιώνει \(N\) επαναλήψεις και υπολογίζει τα αντίστοιχα μέσα μεγέθη για το μέσο χρόνο παραμονής ενός πελάτη στο σύστημα, το μέσο αριθμό πελατών που απορρίπτονται, το μέσο κέρδος στο διάστημα \([0,T]\) και το μέσο αριθμό πελατών στο σύστημα στις χρονικές στιγμές \(t_j, j=1,\ldots,n\). Για τον υπολογισμό της μέσης καθυστέρησης, μετά από κάθε επανάληψη αφαιρούμε από τον πίνακα custdata τις γραμμές που αντιστοιχούν σε πελάτες που έχουν απορριφθεί ή δεν έχουν αναχωρήσει από το σύστημα στο τέλος του χρόνου \(T\) και υπολογίζουμε το μέσο χρόνο παραμονής των πελατών κατά τη διάρκεια της συγκεκριμένης επανάληψης από τη μέση τιμή της διαφοράς των χρόνων άφιξης και αναχώρησης στο μικρότερο πίνακα custdata. Ο εκτιμώμενος μέσος χρόνος είναι η μέση τιμή των \(N\) μέσων χρόνων παραμονής των επιμέρους επαναλήψεων.

# Simulate ggm for N iterations and compute estimates of performance measures
ggmqueueest=function(m, a, l, b1, b2, K, T, c, R, n,N)
{
  syslength=matrix(rep(0,N*n), N,n)
  profit=rep(0,N)
  delay=rep(0,N)
  rejected=rep(0,N)
  for (j in 1:N)
  {
    sout=ggmsimul(m,a,l,b1,b2,K,T,c,R,n)
    syslength[j,]=sout$systemlength
    cdata=sout$custdata[sout$custdata[,4]==0,]
    cdata=cdata[!(is.na(cdata[,3])),]
    delay[j]=mean(cdata[,3]-cdata[,2])
    profit[j]=sout$totalprofit
    rejected[j]=sout$nrejected
  }
  msyslength=apply(syslength,2,mean)
  mdelay=mean(delay)
  mprofit=mean(profit)
  mrejected=mean(rejected)
  estout=list(syslength=msyslength, delay=mdelay, profit=mprofit, rejected=mrejected)
  return(estout)
}

Για εφαρμογή των παραπάνω θεωρούμε μια ουρά με \(m=4\) servers και χρόνους άφιξης Erlang(2,1) δηλαδή \(a=2, l=1\) επομένως ο μέσος χρόνος μεταξύ αφίξεων είναι ίσος με \(a/l=2\) και ο μέσος ρυθμός άφιξης \(\lambda=1/2\). Οι χρόνοι εξυπηρέτησης ακολουθούν \(U(0,14)\), επομένως \(E(S)=7\) και \(\mu=1/7\). Υποθέτουμε απεριόριστη χωρητικότητα \(K=\infty\). Ο συντελεστής ευστάθειας του συστήματος είναι ίσος με \(\rho= \frac{\lambda}{m \mu} = \frac{7}{8}<1\), επομένως το σύστημα είναι ευσταθές και περιμένουμε ο μέσος αριθμός πελατών να σταθεροποιείται έπειτα από μεγάλο χρονικό διάστημα μετά την έναρξη της λειτουργίας του. Το έσοδο από την εξυπηρέτηση είναι \(R=4\) και το κόστος λειτουργίας ανά μονάδα χρόνου ίσο με \(m \mu^2\), δηλαδή \(c=1\).

Για το παραπάνω σύστημα προσομοιώνουμε \(N=1000\) επαναλήψεις διάρκειας \(T=100\). Ο αριθμός ατόμων στο σύστημα καταγράφεται σε \(n=100\) ισαπέχουσες χρονικές στιγμές στο διάστημα \([0,T]\).

#
a=2
l=1
b1=0
b2=14
K=Inf
n=100
T=100
m=4
c=1
R=4
N=1000
#
estout=ggmqueueest(m,a,l,b1,b2,K,T,c,R,n,N)
estout

## $syslength
##   [1] 0.268 0.699 1.150 1.566 1.912 2.270 2.567 2.869 3.059 3.239 3.395 3.520
##  [13] 3.620 3.681 3.713 3.734 3.779 3.836 3.881 3.942 3.987 4.000 4.102 4.150
##  [25] 4.155 4.175 4.233 4.250 4.280 4.283 4.305 4.357 4.389 4.378 4.414 4.403
##  [37] 4.442 4.488 4.524 4.517 4.518 4.577 4.569 4.585 4.608 4.619 4.591 4.601
##  [49] 4.620 4.635 4.656 4.700 4.724 4.727 4.762 4.751 4.783 4.784 4.766 4.789
##  [61] 4.801 4.809 4.793 4.831 4.827 4.833 4.868 4.850 4.862 4.845 4.882 4.875
##  [73] 4.883 4.865 4.885 4.886 4.911 4.908 4.912 4.955 4.951 4.967 4.925 4.992
##  [85] 5.013 4.994 4.982 5.006 4.999 5.025 5.003 4.997 4.997 5.003 5.006 5.039
##  [97] 5.018 5.023 5.008 5.000
## 
## $delay
## [1] 8.824389
## 
## $profit
## [1] 171.0527
## 
## $rejected
## [1] 0

Μπορούμε να δημιουργήσουμε ένα γράφημα του μέσου αριθμού πελατών ως προς \(t\), με βάση το διάνυσμα estout$syslength:

cl=estout$syslength
t=(1:n)*T/n
plot(cl~t, type="l")

Παρατηρούμε ότι ο μέσος αριθμός πελατών όντως συγκλίνει σε μια οριακή τιμή περίπου ίση με 5.5. Αν επαναλάβουμε το πείραμα με \(S~U(0,30)\) οπότε \(\rho>1\), βλέπουμε ότι ο μέσος αριθμός πελατών αυξάνει χωρίς να σταθεροποιείται σε κάποια οριακή τιμή.

b2=30
estout=ggmqueueest(m,a,l,b1,b2,K,T,c,R,n,1000)
estout

## $syslength
##   [1]  0.272  0.749  1.194  1.646  2.090  2.474  2.873  3.253  3.624  3.966
##  [11]  4.325  4.665  5.025  5.421  5.757  6.063  6.346  6.709  7.013  7.327
##  [21]  7.611  7.854  8.137  8.386  8.585  8.834  9.024  9.270  9.489  9.713
##  [31]  9.916 10.111 10.332 10.521 10.731 11.006 11.243 11.451 11.703 11.946
##  [41] 12.174 12.412 12.689 12.944 13.135 13.413 13.643 13.880 14.082 14.316
##  [51] 14.513 14.757 14.979 15.196 15.441 15.708 15.918 16.122 16.361 16.621
##  [61] 16.868 17.111 17.267 17.482 17.702 17.921 18.157 18.367 18.594 18.857
##  [71] 19.098 19.321 19.582 19.809 20.047 20.255 20.430 20.666 20.915 21.151
##  [81] 21.342 21.567 21.793 22.029 22.212 22.466 22.742 22.931 23.132 23.362
##  [91] 23.612 23.866 24.084 24.302 24.555 24.721 24.972 25.230 25.429 25.664
## 
## $delay
## [1] 30.34536
## 
## $profit
## [1] 93.50222
## 
## $rejected
## [1] 0

cl=estout$syslength
t=(1:n)*T/n
plot(cl~t, type="l")

Τέλος, στο παραπάνω σύστημα θέτουμε πεπερασμένη χωρητικότητα π.χ. \(Κ=20\) :

b2=30
K=20
estout=ggmqueueest(m,a,l,b1,b2,K,T,c,R,n,N)
estout

## $syslength
##   [1]  0.294  0.764  1.255  1.693  2.125  2.558  2.985  3.344  3.707  4.061
##  [11]  4.420  4.723  5.059  5.422  5.737  6.050  6.370  6.704  7.009  7.288
##  [21]  7.541  7.802  8.079  8.309  8.562  8.799  9.002  9.257  9.451  9.641
##  [31]  9.864 10.070 10.274 10.489 10.701 10.937 11.148 11.381 11.600 11.812
##  [41] 12.096 12.315 12.517 12.766 13.011 13.211 13.442 13.627 13.862 14.022
##  [51] 14.232 14.442 14.644 14.852 15.031 15.221 15.391 15.556 15.719 15.882
##  [61] 16.034 16.206 16.375 16.505 16.639 16.770 16.865 16.965 17.083 17.193
##  [71] 17.318 17.442 17.528 17.628 17.724 17.815 17.909 17.956 18.035 18.101
##  [81] 18.185 18.258 18.330 18.411 18.475 18.504 18.552 18.599 18.638 18.691
##  [91] 18.715 18.750 18.749 18.778 18.801 18.835 18.825 18.871 18.887 18.891
## 
## $delay
## [1] 30.08456
## 
## $profit
## [1] 93.97022
## 
## $rejected
## [1] 6.732

cl=estout$syslength
t=(1:n)*T/n
plot(cl~t, type="l")

Τώρα ο μέσος αριθμός ατόμων στο σύστημα σταθεροποιείται επειδή ο χώρος καταστάσεων είναι πεπερασμένος και το σύστημα είναι ευσταθές για κάθε τιμή του \(\rho < \infty\).

Άσκηση 4

Θεωρούμε μια Μαρκοβιανή αλυσίδα διακριτού χρόνου \(\{X_n, n=0,1,\ldots\}\) με χώρο καταστάσεων \(S=\{1,\ldots,M\}\), πίνακα πιθανοτήτων μετάβασης ενός βήματος \(P=(p_{ij})_{i,j \in S}\) και κατανομή αρχικής κατάστασης \(p(0) =(p_i(0))_{i \in S}\).

Να δημιουργήσετε τις κατάλληλες συναρτήσεις προσομοίωσης για την εκτίμηση των παρακάτω ποσοτήτων:

1. Κατανομή πιθανότητας της κατάστασης στην περίοδο \(n\): \(p(n)=(p_{j}(n))_{j \in S}\), όπου \(p_j(n) = P(X_n=j)\).

2. Μέσος χρόνος πρώτης μετάβασης από την κατάσταση \(j\) στην κατάσταση 1 \(\tau_j = E(T |X_0=j)\), όπου \(T = \min\{n>0 : X_n=1\}\), ως συνάρτηση του \(j\). Υποθέτουμε ότι η κατάσταση 1 είναι θετικά επαναληπτική. Για \(j=1\) η \(\tau_1\) είναι ο μέσος χρόνος επανόδου στην κατάσταση 1. Η εκτίμηση να γίνει με ένα διάστημα εμπιστοσύνης 95%.

3. Έστω ότι όποτε η διαδικασία βρεθεί στην κατάσταση \(j\) υπάρχει ένα κόστος περιόδου ίσο με \(c(j)\) όπου \(c=(c(j))_{j \in S}\) είναι δοσμένο διάνυσμα. Να εκτιμηθεί το αναμενόμενο μέσο κόστος ανά περίοδο \(C_n\) για τις \(n\) πρώτες περιόδους. Η εκτίμηση να γίνει με ένα διάστημα εμπιστοσύνης 95%.

\[ C_n = E \left ( \frac{1}{n} \sum_{t=0}^{n-1} c(X_t) \right). \] Σε όλα τα παραπάνω η εκτίμηση να βασιστεί σε \(Ν\) επαναλήψεις, και ο αριθμός \(N\) περιλαμβάνεται στα ορίσματα της κάθε συνάρτησης. Για κάθε ερώτημα να δημιουργηθεί μια διαφορετική συνάρτηση που εκτελεί τις προσομοιώσεις και εκτιμά την αντίστοιχη ποσότητα.

Απάντηση

Σε όλα τα ερωτήματα θα χρειαστεί μια γεννήτρια τυχαίων αριθμών από μια διακριτή κατανομή στους ακέραιους \(\{1, \ldots, M\}\) με διάνυσμα πιθανότητας \(p\). Η παρακάτω συνάρτηση υλοποιεί τη γεννήτρια αντίστροφου μετασχηματισμού και θα χρησιμοποιηθεί για την προσομοίωση των μεταβάσεων της Μαρκοβιανής αλυσίδας. \[ X=\max\{n : p_1+\cdots+p_n \leq U \} \] όπου \(U \sim U(0,1)\).

gendiscrdistr=function(p)
{
  M=length(p)
  j=1
  s=p[1]
  u=runif(1,0,1)
  while(s<=u)
  {
    j=j+1
    s=s+p[j]
  }
  return(j)
}

Ερώτημα 1 Η παρακάτω συνάρτηση δημιουργεί \(N\) υλοποιήσεις της διαδικασίας μήκους \(n\) περιόδων η κάθε μια, με αρχική κατάσταση από την κατανομή \(p(0)\). Στο τέλος κάθε υλοποίησης καταγράφεται η κατάσταση \(X_n\) και τελικά η κατανομή της \(X_n\) εκτιμάται μέσω της εμπειρικής κατανομής των δεδομένων \((X_n)_k, k=1,\ldots,N\).

markovchainnstepdistr=function(p0, p, n, N)
{
  # Simulate a discrete time MC with transition matrix p and initial distribution p0. Compute empirical distribution of X_n
  M=length(p0)
  # Generate initial state from p0
  xn=rep(0,N)
  for (k in 1:N)
  {
      # Generate initial state from p0
      x=gendiscrdistr(p0)
      for (t in 1:n)
      {
        # Transition probabilities from row x of p
        px=p[x,]
        x=gendiscrdistr(px)
      }
      xn[k]=x
  }
  # Compute empirical frequencies of X_n
  pnest=table(xn)/N
  return(pnest)
}

Για έλεγχο της συνάρτησης, θεωρούμε μια Μαρκοβιανή αλυσίδα με \(M=3\), πίνακα ενός βήματος \[ P=\left ( \begin{array}{ccc} 1/2&1/4&1/4\\ 1/3&1/3&1/3\\ 0&1/4&3/4\\ \end{array} \right ) \] και διάνυσμα αρχικής κατανομής \(p(0)=(1/4, 1/2, 1/4)\).

P=matrix(c(1/2, 1/4, 1/4, 1/3, 1/3, 1/3, 0, 1/4, 3/4), 3, 3, byrow=T)
p0=c(1/4, 1/2, 1/4)

Από τη θεωρία Μαρκοβιανών Αλυσίδων γνωρίζουμε ότι η κατανομή της κατάστασης την περίοδο \(n\) δίνεται από τον τύπο \(p(n)^T = p(0)^T P^n\). Επομένως, π.χ. για \(n=5\),

n=5
pn=p0
for (t in 1:n) pn=pn%*%P
pn

##         [,1]      [,2]      [,3]
## [1,] 0.18963 0.2727282 0.5376419

Τώρα προσομοιώνουμε \(N=1000\) υλοποιήσεις με αρχική κατανομή \(p(0)\) χρησιμοποιώντας τη συνάρτηση markovchainnstepdistr

n=5
N=1000
pnest=markovchainnstepdistr(p0,P,n,N)
pnest

## xn
##     1     2     3 
## 0.176 0.272 0.552

Η εκτίμηση είναι αρκετά κοντά στην ακριβή κατανομή.

Ερώτημα 2 Η παρακάτω συνάρτηση δημιουργεί \(N\) υλοποιήσεις της διαδικασίας με αρχική κατάσταση \(j\) και διάρκεια μέχρι την περίοδο της πρώτης επίσκεψης στην κατάσταση 1. Στο τέλος κάθε υλοποίησης καταγράφεται ο αριθμός βημάτων μέχρι το σταμάτημα.

markovchainfirstpassagetime=function(p, j, N)
{
  # Simulate a discrete time MC with transition matrix p and initial state j. Estimate mean first passage time from j to 1.
  t1=rep(0,N)
  for (k in 1:N)
  {
    # Generate initial state from p0
    n=1
    px=p[j,]
    x=gendiscrdistr(px)
    while (!(x==1))
    {
      n=n+1
      px=p[x,]
      x=gendiscrdistr(px)
    }
    t1[k]=n
  }
  # Compute empirical frequencies of X_n
  t1est=mean(t1)
  sd1=sd(t1)
  hl=qt(0.975, N-1)*sd1/sqrt(N)
  cint=c(t1est-hl, t1est, t1est+hl)
  return(cint)
}

Εφαρμόζουμε τη συνάρτηση για να υπολογίσουμε τα διαστήματα εμπιστοσύνης για τα \(\tau_1, \ldots, \tau_3\) από το προηγούμενο παράδειγμα. Οι ακριβείς τιμές μπορούν να υπολογιστούν λύνοντας κατάλληλο σύστημα γραμμικών εξισώσεων και είναι \(\tau_1=11/2, \tau_2=7, \tau_3=11\).

M=length(p0)
N=1000
tau=matrix(rep(0,3*M), M,3)
for (j in 1:M)
{
  tau[j,]=markovchainfirstpassagetime(P, j, N)
}
tau

##           [,1]   [,2]      [,3]
## [1,]  4.973151  5.451  5.928849
## [2,]  6.432756  6.931  7.429244
## [3,] 10.601734 11.162 11.722266

Ερώτημα 3 Για την εκτίμηση του αναμενόμενου μέσου κόστους ανά περίοδο για τις πρώτες \(n\) περιόδους \(C_n\) δημιουργούμε την παρακάτω συνάρτηση που προσομοιώνει \(N\) υλοποιήσεις της διαδικασίας με αρχική κατανομή \(p(0)\) και υπολογίζει το διάστημα εμπιστοσύνης για το \(C_n\):

markovchainavcost=function(p0, p,c, n,N)
{
  # Simulate a discrete time MC with transition matrix p, initial distribution p0 and one period cost function c(j), j in S. Estimate average cost per period for periods 0,1,...n-1 .
  # Generate initial state from p0
  cn=rep(0,N)
  for (k in 1:N)
  {
      # Generate initial state from p0
      x=gendiscrdistr(p0)
      totalcost=c[x]
      for (t in 1:n-1)
      {
        # Transition probabilities from row x of p
        px=p[x,]
        x=gendiscrdistr(px)
        totalcost=totalcost+c[x]
      }
      cn[k]=totalcost/n
  }
  # Compute confidence interval for cn
  cnest=mean(cn)
  cnsd=sd(cn)
  hl=qt(0.975, N-1)*cnest/sqrt(N)
  cint=c(cnest-hl, cnest, cnest+hl)
  return(cint)
}

Εφαρμόζουμε τη συνάρτηση στο προηγούμενο παράδειγμα με τρεις καταστάσεις και συνάρτηση κόστους μιας περιόδου \(c(1)=15, c(2)=30, c(3)=50\) και εκτιμούμε το μέσο κόστος 20 περιόδων:

c=c(15,30,50)
n=20
N=1000
costint=markovchainavcost(p0,P,c,n, N)
costint

## [1] 36.71446 39.14350 41.57254

Άσκηση 5

Θεωρούμε μια αδιαχώριστη θετικά επαναληπτική Μαρκοβιανή αλυσίδα διακριτού χρόνου \(\{X_n, n=0,1,\ldots\}\) με χώρο καταστάσεων \(S=\{1,\ldots,M\}\) και πίνακα πιθανοτήτων μετάβασης ενός βήματος \(P=(p_{ij})_{i,j \in S}\).

1 Να εκτιμηθεί μέσω προσομοίωσης η οριακή κατανομή \(\pi=(\pi_1, \ldots, \pi_M)\).

2 Να εκτιμηθεί το μέσο κόστος ανά μονάδα χρόνου σε άπειρο ορίζοντα.

Απάντηση Η ακριβής τιμή της στάσιμης κατανομής μπορεί επίσης να υπολογιστεί ως εξής: Το σύστημα εξισώσεων \(\pi = \pi P, \sum p_j =1\) μπορεί να γραφτεί πινακικά ως εξής: \[ (I-P)^T \pi=0\\ 1^T \pi = 1 \] ή ισοδύναμα \(A \pi = b\), όπου \(A\) ο πίνακας που προκύπτει από τον \((I-P)^T\) αν αντικατασταθεί η τελευταία γραμμή με το διάνυσμα \(1^Τ\) και \(b = (0,0,\ldots,1)\). Αφού δημιουργήσουμε αυτά τα στοιχεία, μπορούμε να λύσουμε το γραμμικό σύστημα με την εντολή solve:

A=t(diag(3) - P)
A=rbind(A[1:2,], c(1,1,1))
A

##       [,1]       [,2]  [,3]
## [1,]  0.50 -0.3333333  0.00
## [2,] -0.25  0.6666667 -0.25
## [3,]  1.00  1.0000000  1.00

b=c(0,0,1)
b

## [1] 0 0 1

pi=solve(A,b)
pi

## [1] 0.1818182 0.2727273 0.5454545

Επίσης η ακριβής τιμή του μέσου κόστους σε άπειρο ορίζοντα είναι ίση με \(C=\sum_j \pi_j c_j\), επομένως

avcost= sum(pi*c)
avcost

## [1] 38.18182

Τώρα επιστρέφουμε στην προσομοίωση. Ένας πρώτος άμεσος τρόπος να υπολογιστούν προσεγγιστικά τα παραπάνω μεγέθη μέσω προσομοιίωσης είναι να χρησιμοποιήσουμε την προηγούμενη άσκηση θέτοντας ένα μεγάλο αριθμό στο μήκος κάθε υλοποίησης. Αυτό βασίζεται στην ιδότητα ότι σε μια αδιαχώριστη θετικά επαναληπτική αλυσίδα η οριακή κατανομή είναι \(\pi_j = \lim_{n \to \infty} p_{j}(n)\) για οποιαδήποτε αρχική κατανομή \(p(0)\). Όμοια το μέσο κόστος ανά μονάδα χρόνου σε άπειρο ορίζοντα είναι ίσο με \(C = \lim_{n \to \infty C_n}\). Για παράδειγμα για την Μαρκοβιανή αλυσίδα της προηγούμενης άσκησης θέτοντας \(n=1000\) και αριθμό υλοποιήσεων \(N=1000\):

n=1000
N=1000
piest=markovchainnstepdistr(p0,P,n,N)
piest

## xn
##     1     2     3 
## 0.198 0.265 0.537

avcostint=markovchainavcost(p0,P,c,n,N)
avcostint

## [1] 35.85682 38.22912 40.60142

Βλέπουμε ότι έχει επιτευχθεί σχετικά καλή προσέγγιση, αυτό όμως έγινε με μεγάλο υπολογιστικό κόστος γιατί χρειάστηκε να υλοποιηθούν 1000 υλοποιήσεις μήκους 1000 περιόδων η κάθε μια.

Ένας εναλλακτικός τρόπος προσομοίωσης είναι να εκμεταλλευτούμε την εργοδικότητα, δηλαδή την ερμηνεία της στάσιμης κατανομής ως ποσοστών χρόνου σε κάθε κατάσταση και το μέσο κόστος ως αντίστοιχο όριο μιας υλοποίησης με πιθανότητα 1: \[ \pi_j = \lim_{n \to \infty} \frac{1}{n} \sum_{t=0}^{n-1} 1(X_t=j), j \in S\\ C= \lim_{n \to \infty} \frac{1}{n} \sum_{t=0}^{n-1} c(X_t) \] με πιθανότητα 1. Επομένως μπορούμε να εκτιμήσουμε τα παραπάνω μεγέθη προσομοιώνοντας μόνο μια υλοποίηση της διαδικασίας με μεγάλο μήκος ορίζοντα. Η παρακάτω συνάρτηση κάνει αυτή την προσομοίωση και χρησιμοποιεί και ένα αρχικό burnin period \(r n\) περιόδων όπου \(r \in [0,1)\) είναι μια παράμετρος επιλογής μας. Ο σκοπός του burnin είναι το μέτρημα των στατιστικών στοιχείων να αρχίσει μετά από τις πρώτες \(\lfloor r n \rfloor\) περιόδους (\(t=0,\ldtos,\lfloor r n \rfloor)-1\)), έτσι ώστε να δοθεί στο σύστημα χρόνος να συγκλίνει στη στάσιμη κατανομή έτσι ώστε οι εκτιμήσεις να μην επηρεάζονται από τη μεταβατική κατανομή των αρχικών περιόδων. Με βάση αυτό οι εκτιμήσεις των μεγεθών θα είναι : \[ \hat{\pi}_j = \frac{1}{n - \lfloor r n \rfloor} \sum_{t=\lfloor r n \rfloor }^{n-1} 1(X_t=j), j \in S\\ \hat{C}= \frac{1}{n - \lfloor r n \rfloor} \sum_{t=\lfloor r n \rfloor}^{n-1} c(X_t) \]

picostergodicestimation=function(p0,p,c,n,r)
{
  # Simulate a discrete time MC with transition matrix p, initial distribution p0 and one period cost function c(j), j in S. Estimate stationary probabilities and average cost per period as corresponding time averages from a single realization of length n using an initial burnin period r n
  # Generate initial state from p0
  M=length(p0)
  # Generate initial state from p0
  x=gendiscrdistr(p0)
  n0=floor(r*n)    
  # Simulate first n0 transitions 0,...,n0-1 without keeping count 
  for (t in 1:(n0-1))
      {
        # Transition probabilities from row x of p
        px=p[x,]
        x=gendiscrdistr(px)
      }
  # Continue simulation keeping statistics
  statefrequency=rep(0,M)
  totalcost=0
  for (t in n0:(n-1))
      {
        # Transition probabilities from row x of p
        px=p[x,]
        x=gendiscrdistr(px)
        statefrequency[x]=statefrequency[x]+1
        totalcost=totalcost+c[x]
      }  # Compute confidence interval for cn
piestimate=statefrequency/(n-n0)
costestimate=totalcost/(n-n0)
return(list(piest=piestimate, avcostest=costestimate))
}

Για εφαρμογή χρησιμοποιούμε τις προηγούμενες παραμέτρους και καλούμε τη συνάρτηση με \(n=10000\) και \(r=0.1\):

n=10000
r=0.1
res=picostergodicestimation(p0,P,c,n,r)
res$piest

## [1] 0.1830000 0.2707778 0.5462222

res$avcostest

## [1] 38.17944

Τώρα έχουμε επιτύχει παρόμοια προσέγγιση με την προηγούμενη μέθοδο, αλλά με το 1/100 των επαναλήψεων.