Άσκηση 1 (Ross 4.1)

Έστω \(p_j = P(X=j), j=1,2,3,4\). Έχουμε \(p_3 > p_1 > p_2 > p_4\). Δημιουργούμε τη γεννήτρια αντίστροφου μετασχηματισμού, με \(U \sim U(0,1)\): \[ X=\left \{ \begin{array}{ll} 3, & 0 \leq U \leq 0.35 \\ 1, & 0.35 <U \leq 0.65 \\ 2, & 0.65 <U \leq 0.85 \\ 4, & 0.85 <U \leq 1 \end{array} \right . \] Η παρακάτω συνάρτηση δημιουργεί \(n\) παρατηρήσεις από αυτή τη γεννήτρια:

genx=function(n)
{
  x=rep(0,n)
  for (j in 1:n)
  {
    u=runif(1)
    if (u <= 0.35) x[j]=3 else 
      if (u <= 0.65) x[j]=1 else
        if (u <= 0.85) x[j]=2 else 
          x[j]=4
  }
  return(x)
}

Για επαλήθευση καλούμε αυτή τη συνάρτηση με \(n=1000\) και δημιουργούμε ένα πίνακα συχνοτήτων των προσομοιωμένων τιμών:

x=genx(1000)
freq=table(x)/1000
freq
## x
##     1     2     3     4 
## 0.301 0.218 0.330 0.151

Άσκηση 2 (Ross 4.4)

Έστω \(n\) ο αριθμός των παικτών και \(v=(v_1, \ldots, v_n)\) το διάνυσμα των παραμέτρων των παικτών. Για να προσομοιώσουμε τη διαδικασία του διαγωνισμού, θεωρούμε αρχικά το διάνυσμα \(x = 1:n\) το οποίο προσαρμόζουμε μετά από κάθε παιχνίδι ως εξής. Yποθέτουμε ότι στο βήμα \(k\) της διαδικασίας (πριν παιχτεί το παιχνίδι \(k\)) οι θέσεις \(1:k-1\) του διανύσματος \(x\) περιέχουν τους παίκτες που έχουν ήδη αποκλειστεί. Επιπλέον, αν στο παιχνίδι \(k\) αποκλειστεί ο παίκτης ο οποίος βρίσκεται στη θέση \(j\) του διανύσματος \(x\), τότε αυτός αλλάζει θέση με τον παίκτη που βρίσκεται στη θέση \(k\). Επομένως μετά το παιχνίδι \(k\) οι θέσεις \(1:k\) περιέχουν τους παίκτες που έχουν αποκλειστεί. Για την προσομοίωση του παιχνιδιού \(k\) επιλέγονται τυχαία δύο θέσεις \(i,j \in \{k, \ldots, n\}\) και κερδίζει ο \(x_i\) με πιθανότητα \(\frac{v_{x_i}}{v_{x_i}+v_{x_j}}\). (Για την επιλογή δύο τυχαίων αριθμών από το σύνολο \(k:n\) χρησιμοποιούμε τη συνάρτηση sample.)

Για παράδειγμα αν \(n=5\), στο πρώτο παιχνίδι κληρώθηκαν οι 2,5 και αποκλείστηκε ο 5 και στο δεύτερο κληρώθηκαν οι 4,2 και αποκλείστηκε ο 2, μετά από το πρώτο παιχνίδι το διάνυσμα \(x\) θα είναι \(x=(5,2,3,4,1)\) και μετά από το δεύτερο παιχνίδι θα είναι \(x=(5,4,3,2,1)\).

Επομένως στο τέλος του διαγωνισμού ο νικητής θα είναι ο παίκτης που βρίσκεται στη θέση \(n\) του διανύσματος \(x\). Για την προσομοίωση μιας επανάληψης του διαγωνισμού δημιουργούμε την παρακάτω συνάρτηση

simcompetition=function(v)
{
  n=length(v)
  x=1:n
  for (k in 1:(n-1))
  {
    pair=sample(k:n, 2)
    i=pair[1]
    j=pair[2]
    winprobi=v[x[i]]/(v[x[i]]+v[x[j]])
    if (runif(1,0,1)<=winprobi) # j loses
    {
      a=x[k]
      x[k]=x[j] 
      x[j]=a
    } else # i loses
    {
      a=x[k]
      x[k]=x[i]
      x[i]=a
    }
  }
  return(x)
}

Η συνάρτηση επιστρέφει το διάνυσμα \(x\) με τη σειρά των παικτών που αποκλείστηκαν. Για παράδειγμα έστω \(n=5\) και διάνυσμα παραμέτρων \(v=(5,6,1,2,3)\).

v=c(5,6,1,2,3)
x=simcompetition(v)
x
## [1] 1 3 5 4 2

Εδώ αποκλείστηκε πρώτος ο παίκτης

x[1]
## [1] 1

ενώ νικητής του διαγωνισμού ήταν ο

x[length(x)]
## [1] 2

Έχοντας τη συνάρτηση που προσομοιώνει ένα διαγωνισμό, μπορούμε να εκτιμήσουμε την πιθανότητα να κερδίσει τον διαγωνισμό π.χ. ο παίκτης 1, εκτελώντας \(Ν\) επαναλήψεις και εκτιμώντας την πιθανότητα μέσω του αντίστοιχου ποσοστού. Για παράδειγμα έστω ότι ο παίκτης 1 έχει πολύ μεγαλύτερη τιμή της παραμέτρου του από τους άλλους, επομένως έχει το πλεονέκτημα. Αν υπάρχουν λίγοι παίκτες βλέπουμε

N=10000
v=c(10,1,1)
n=length(v)
results=rep(0,N)
for (j in 1:N)
{
  x=simcompetition(v)
  results[j]=as.numeric(x[n]==1)
}
pwin1=mean(results)
pwin1
## [1] 0.8558

επομένως ο 1 έχει πολύ μεγάλη πιθανότητα να κερδίσει το διαγωνισμό. Όμως αν ο αριθμός των άλλων (αδύνατων) παικτών αυξηθεί βρίσκουμε ότι η πιθανότητα μειώνεται. Για παράδειγμα με 20 παίκτες συνολικά:

N=10000
v=c(10,rep(1,19))
n=length(v)
results=rep(0,N)
for (j in 1:N)
{
  x=simcompetition(v)
  results[j]=as.numeric(x[n]==1)
}
pwin1=mean(results)
pwin1
## [1] 0.6106

Άσκηση 3 (Ross 4.11)

Θα δούμε πρώτα την προσομοίωση επιλογής ενός υποσυνόλου \(r\) στοιχείων από ένα σύνολο \(n\) στοιχείων χωρίς επανάθεση (επιλογή ενός τυχαίου συνδυασμού).

Έστω το σύνολο \(\{1, 2, \ldots, n\}\). Ακολουθώντας το παράδειγμα (4b), η επιλογή συνδυασμού μπορεί να υλοποιηθεί ακολουθώντας τον αλγόριθμο δημιουργίας μιας τυχαίας μετάθεσης \(n\) στοιχείων και σταματώντας όταν έχουν επιλεγεί τα στοιχεία που θα μπουν στις \(r\) πρώτες θέσεις. Επομένως ξεκινάμε με το αρχικό διάνυσμα \(x=(1,2,\ldots,n)\) και στο βήμα \(k\) επιλέγουμε τυχαία ένα στοιχείο από τα \(x_{k}, x_{k+2}, \ldots, x_n\) και το εναλλάσσουμε με το στοιχείο \(x_k\). Προφανώς αν επιλεγεί το \(x_k\) παραμένει στη θέση του και προχωρούμε στο επόμενο βήμα. Επαναλαμβάνουμε αυτή τη διαδικασία για \(r\) βήματα και παίρνουμε τα πρώτα \(r\) στοιχεία του διανύσματος \(x\):

randomcombination=function(n,r)
{
  x=1:n
  for (k in 1:r)
  {
    j=sample(k:n, 1)
    a=x[k]
    x[k]=x[j]
    x[j]=a
  }
  y=x[1:r]
  return(y)
}
randomcombination(10,4)
## [1] 3 7 2 1

Για τη δημιουργία ενός τυχαίου συνδυασμού \(r\) στοιχείων από \(n\), δεδομένου ότι περιέχει τουλάχιστον ένα στοιχείο από τα πρώτα \(k\), χρησιμοποιούμε τη μέθοδο αποδοχής-απόρριψης. Συγκεκριμένα επιλέγουμε ένα συνδυασμό χρησιμοποιώντας την παραπάνω συνάρτηση και αν αυτός περιέχει τουλάχιστον ένα από τα στοιχεία \(1,\ldots,k\) τότε γίνεται δεκτός, διαφορετικά απορρίπτεται και επαναλαμβάνουμε τη διαδικασία. Παρατηρούμε ότι ένα διάνυσμα ακεραίων από τους \(1,\ldots,n\) περιέχει τουλάχιστον έναν ακέραιο από τους \(1,\ldots,k\) αν και μόνο αν \(\min(x) \leq k\).

randomcombinationconditional=function(n,r,k)
{
  x=randomcombination(n,r)
  while (min(x)>k)
    x=randomcombination(n,r)
  return(x)
}
randomcombinationconditional(10, 4,3)
## [1] 8 3 2 1

Άσκηση 4 (Ross 4.13)

Παρατηρούμε ότι \[ P(X=i) = P(Y=i | Y \leq k), \] όπου \(Y\sim Poisson(\lambda)\). Η πρώτη γεννήτρια που θα εφαρμόσουμε είναι μέσω αντίστροφου μετασχηματισμού. Για δεδομένα \(\lambda, k\) μπορούμε να υπολογίσουμε το διάνυσμα πιθανότητας \(p_i = P(X=i)\) της \(X\) και να εφαρμόσουμε τη μέθοδο του αντίστροφου μετασχηματισμού. Η επόμενη συνάρτηση εφαρμόζει αυτή τη μέθοδο για τη δημιουργία \(n\) παρατηρήσεων από αυτή την κατανομή. Ο λόγος που ενσωματώνεται το \(n\) στη συνάρτηση είναι έτσι ώστε αν απαιτούνται περισσότερες παρατηρήσεις αυτές να δημιουργούνται απευθείας με ένα μόνο υπολογισμό του διανύσματος \(p\) και όχι να υπολογίζεται αυτό κάθε φορά που δημιουργείται μια παρατήρηση:

gen_cond_Poisson_1=function(n,l,k)
{
  p=rep(0,k+1)
  P=ppois(k,l)
  p[1]=exp(-l)
  for (j in 2:(k+1))
    p[j]=p[j-1]*l/(j-1)
  p=p/P
  x=rep(0,n)
  for (j in 1:n)
  {
    i=1
    s=p[1]
    u=runif(1,0,1)
    while(s<u)
    {
      i=i+1
      s=s+p[i]
    }
    x[j]=i-1
  }
  return(x)
}
gen_cond_Poisson_1(10,4,8)
##  [1] 6 5 1 3 3 1 4 4 6 5
gen_cond_Poisson_1(10,20,8)
##  [1] 7 8 8 8 8 8 8 8 7 8

Η δεύτερη γεννήτρια βασίζεται στη μέθοδο αποδοχής-απόρριψης και περιγράφεται ως εξής: Δημιουργούμε μια παρατήρηση \(Y\) από την Poisson($). Αν \(Y\leq k\) τη δεχόμαστε και τότε \(X=Y\), διαφορετικά την απορρίπτουμε και επαναλαμβάνουμε τη διαδικασία.

gen_cond_Poisson_2=function(n,l,k)
{
  x=rep(0,n)
  for (j in 1:n)
  {
    y=rpois(1,l)
    while (y > k)
    {
      y=rpois(1,l)
    }
    x[j]=y
  }
  return(x)
}
gen_cond_Poisson_2(10,4,8)
##  [1] 3 2 3 2 5 3 5 6 6 3
gen_cond_Poisson_2(10,20,8)
##  [1] 8 8 8 7 7 7 8 7 8 8

Για να ελέγξουμε την αξιοπιστία των δύο γεννητριών, έστω \(\lambda=10, k=7\). Τότε η συνάρτηση πιθανότητας της \(X\) είναι

l=10
k=7
p=rep(0,k+1)
P=ppois(k,l)
p[1]=exp(-l)
for (j in 2:(k+1))
  p[j]=p[j-1]*l/(j-1)
p=p/P
p
## [1] 0.0002061566 0.0020615655 0.0103078277 0.0343594258 0.0858985644
## [6] 0.1717971289 0.2863285481 0.4090407830

Τώρα δημιουργούμε 10000 παρατηρήσεις από κάθε γεννήτρια και υπολογίζουμε τους πίνακες συχνοτήτων, τους οποίους μπορούμε να συγκρίνουμε με τη θεωρητικά υπολογισμένη συνάρτηση πιθανότητας:

x1=gen_cond_Poisson_1(1000,l,k)
x2=gen_cond_Poisson_2(1000,l,k)
freq1=table(x1)/1000
freq2=table(x2)/1000
freq1
## x1
##     2     3     4     5     6     7 
## 0.012 0.038 0.091 0.170 0.298 0.391
freq2
## x2
##     0     2     3     4     5     6     7 
## 0.001 0.013 0.039 0.091 0.178 0.275 0.403

Άσκηση 5 (Ross 4.17)

Έστω η συνάρτηση πιθανότητας της τυχαίας μεταβλητής \(X\), η οποία μπορεί να γραφτεί ως εξής: \[\begin{align} p(j) &=& \left ( \frac{1}{2} \right)^{j+1} + \frac{1}{2} \frac{1}{3} \left ( \frac{2}{3} \right)^{j-1} \\ &=& \frac{1}{2} \left[ \frac{1}{2} \left ( \frac{1}{2} \right)^{j-1} \right ]+ \frac{1}{2} \left [ \frac{1}{3} \left ( \frac{2}{3} \right)^{j-1} \right ] \\ &=& \frac{1}{2} P(X_1=j) + \frac{1}{2} P(X_2=j), j=1,2,\ldots \end{align}\] όπου \(X_1, X_2\) ανεξάρτητες τυχαίες μεταβλητές που ακολουθούν γεωμετρική κατανομή (αριθμό δοκιμών), με πιθανότητα επιτυχίας 1/2 και 1/3 αντίστοιχα. Eπομένως η \(p\) είναι μια ισοπίθανη μείξη των δύο κατανομών. Για τη δημιουργία γεννήτριας μπορούμε να δημιουργήσουμε \(N\) Bernoulli με πιθανότητα 1/2, \(U\sim U(0,1)\) και θέτουμε \[ X= \left \{ \begin{array}{ll} \left \lfloor \frac{\log(U)}{\log(1/2)} \right \rfloor, & N=0 \\ \left \lfloor \frac{\log(U)}{\log(2/3)} \right \rfloor, & N=1 \end{array} \right . \] Η παρακάτω συνάρτηση γενικεύει την κατανομή ως μείξη \(n\) γεωμετρικών με πιθανότητες επιτυχίας \(p_1,\ldots, p_n\) και πιθανότητες μείξης \(a_j = P(X = X_j), j=1, \ldots, n\). Δημιουργεί την \(N \in \{1,2,\ldots, n\}\) από την κατανομή \(a\) μέσω αντίστροφου μετασχηματισμού και αν \(N=j\) δημιουργεί παρατήρηση από την \(Geom(p_j)\).

geommix=function(a, p)
{
  u=runif(1,0,1)
  j=0
  s=0
  while (s<u)
  {
    j=j+1
    s=s+a[j]
  }
  v=runif(1,0,1)
  x=floor(log(v)/log(1-p[j]))
  return(x)
}

Εφαρμόζουμε τη σνάρτηση για \(a=(1/2, 1/2), p=(1/2, 1/3)\) και δημιουργούμε \(n=10\) παρατηρήσεις

n=10
a=c(1/2, 1/2)
p=c(1/2, 1/3)
x=rep(0,n)
for (j in 1:n)
  x[j]=geommix(a,p)
x
##  [1] 0 1 0 0 4 1 2 0 8 0

Άσκηση 6 (Ross 4.18)

(a) Η πιθανότητα \(P(X=j)\) για \(j\geq1\) μπορεί να εκφραστεί ως \[ P(X=j) = P(X > j-1) P(X=j | X>j-1) = \lambda_j P(X>j-1) \] ενώ η πιθανότητα \(P(X>j)\) \[ P(X>j) = P(X>j-1) P(X>j | X>j-1). \] Όμως δεδομένου \(X>j-1\) ισχύει ένα από τα \(X=j\) και \(X>j\), δηλαδή \[ P(X>j | X>j-1) = 1- P(X=j | X>j-1) = 1-\lambda_j, \] επομένως \[ P(X>j) = (1-\lambda_j) P(X>j-1), j=1,2,\ldots \] Επειδή \(P(X>0)=1\), ισχύει επίσης ότι \(P(X=1) = P(X=1 | X>0) = \lambda_1\) και \(P(X>1) = 1- \lambda_1\). Επομένως από την παραπάνω αναδρομική σχέση προκύπτει εύκολα ότι \[ P(X>n) = \prod_{j=1}^{n} (1-\lambda_j) \] και \[ P(X=n) = P(X>n-1) - P(X>n) = \prod_{j=1}^{n-1} (1-\lambda_j) - \prod_{j=1}^{n} (1-\lambda_j) = \prod_{j=1}^{n-1} (1-\lambda_j) \lambda_n \]

(b) Με βάση το (a), μπορούμε να δημιουργήσουμε μια παρατήρηση από την κατανομή της \(X\) προσομοιώνοντας διαδοχικά τα γεγονότα \(X=1, X=2|X>1, X=3|X>2, \ldots,\) μέχρι το πρώτο γεγονός από αυτά που θα συμβεί. Κάθε ένα από τα γεγονότα αυτά έχει πιθανότητα να συμβεί \(P(X=j|X>j-1) = \lambda_j\). Η παρακάτω συνάρτηση υλοποιεί την διαδικασία αυτή. Η λογική μεταβλητή greater χρησιμεύει για να δώσει τη συνθήκη σταματήματος στο while loop, την πρώτη φορά που θα συμβεί το γεγονός \(X=j\). Η μεταβλητή \(F\) υπολογίζει διαδοχικά τις \(F(j) = P(X\leq j)\) και τις χρησιμοποιεί στον υπολογισμό των \(\lambda_j\). Παρατηρούμε ότι αν η \(Χ\) είναι φραγμένη με \(P(X<=M) =1\) για κάποιο \(M < \infty\), τότε \(\lambda_M=1\) επομένως τουλάχιστον το γεγονός \(X=M | X>M-1\) έχει πιθανότητα 1, επομένως ο αλγόριθμος οπωσδήποτε θα σταματήσει σε πεπερασμένο αριθμό βημάτων.

gen_disc_hazard_rate= function(p)
{
  j=0
  F=0
  greater=TRUE
  while (greater)
  {
    j=j+1
    l=p[j]/(1-F)
    F=F+p[j]
    if (runif(1,0,1) <= l) greater=FALSE
  }
  return(j)
}

Έστω η \(X\) με συνάρτηση πιθανότητας \(p_1=1/2, p_2=1/4, p_3=p_4=1/8\). Η συνάρτηση του διακριτού ρυθμού κινδύνου είναι \(\lambda_1 = p_1 = 1/2, \lambda_2 = p_2/(1-p_1) = 1/2, \lambda_3 = p_3/(1-p_1-p_2) = 1/2, \lambda_4 = p_4 /(1-p_1 - p_2 -p_3 = 1)\). Για επαλήθευση της γεννήτριας δημιουργούμε \(n=1000\) παρατηρήσεις από αυτή την κατανομή μέσω της συνάρτησης gen_disc_hazard_rate και υπολογίζουμε τον πίνακα συχνοτήτων :

p=c(1/2, 1/4, 1/8, 1/8)
n=1000
x=rep(0,n)
for (j in 1:n) x[j]=gen_disc_hazard_rate(p)
freq=table(x)/n
freq
## x
##     1     2     3     4 
## 0.508 0.244 0.125 0.123

Άσκηση 7 (Ross 5.7)

(a) Η \(F(x)\) γράφεται ως \[ F(x) = \frac{1}{3} F_1(x) + \frac{1}{3} F_2(x) + \frac{1}{3} F_3(x) \] όπου \(F_1(x) = x, F_2(x) = x^3, F_3(x)=x^5, x \in [0,1]\). Επομένως πρόκειται για μείξη των τριών κατανομών με ισες πιθανότητες μείξης.

genmix1=function()
{
  u=runif(1,0,1)
  N=floor(3*u+1)
  if (N==1) X=runif(1,0,1) else if (N==2) X=runif(1,0,1)^(1/3) else X=runif(1,0,1)^(1/5)
  return(X)
}

Για επαλήθευση της γεννήτριας δημιουργούμε ένα διάνυσμα \(Χ\) με \(n=1000\) παρατηρήσεις και σχηματίζουμε τα plots της θεωρητικής και της εμπειρικής συνάρτησης κατανομής. Ορίζουμε ένα σύνολο τιμών του \(x\) στο διάστημα \([0,1]\), έστω \(x=(0, 0.01, 0.02, \ldots, 1)\) δηλαδή \(x_j = (j-1)/100, j=1,\ldots,101\). Για την εμπειρική συνάρτηση κατανομής, υπολογίζουμε το ποσοστό των στοιχείων του \(X\) που είναι μικρότερα από \(x_j, j=1,\ldots, 101\): \[ \hat{F}(x_j) = \frac{1}{n} \sum_{k=1}^n 1(X_k \leq x_j), j=1,\ldots,101. \]

n=1000
X=rep(0,n)
for (j in 1:n) X[j]=genmix1()
x=seq(0,1,0.01)
nx=length(x)
trueF=(x+x^3+x^5)/3
hatF=rep(0,nx)
for (j in 1:nx) hatF[j] = mean(X<x[j])
plot(trueF~x, col="blue", type="l")
lines(hatF~x, col="red")
legend(0, 0.8, legend=c("true", "genmix1"), col=c("blue","red"), lty=c(1,1), cex=0.8)

Παρατηρούμε ότι οι δύο κατανομές σχεδόν ταυτίζονται, επομένως η γεννήτρια είναι αξιόπιστη.

(b) Η \(F(x), x\geq 0\) γράφεται ως \[ F(x) = \frac{1}{3} F_1(x) + \frac{2}{3} F_2(x), \] όπου \[ F_1(x) = 1-e^{-2x}, \ \ F_2(x) = \left \{ \begin{array}{ll} x,&0 \leq x \leq 1\\ 1,& x>1 \end{array} \right . \]

Επομένως η κατανομή της \(X\) είναι μείξη της \(X_1\sim Exp(2)\) και της \(X_2\sim U(0,1)\) με πιθανότητες 1/3 και 2/3.

genmix2=function()
{
  if (runif(1,0,1) <= 1/3) X=-log(runif(1,0,1))/2 else X=runif(1,0,1)
}

Κάνουμε και εδώ γραφικές παραστάσεις της θεωρητικής και της εμπειρικής κατανομής:

n=1000
X=rep(0,n)
for (j in 1:n) X[j]=genmix2()
x=seq(0,5,0.05)
nx=length(x)
trueF=(1/3)*(1-exp(-2*x))+(2/3)*(x*(x<=1)+(x>1))
hatF=rep(0,nx)
for (j in 1:nx) hatF[j] = mean(X<x[j])
plot(trueF~x, col="blue", type="l")
lines(hatF~x, col="red")
legend(2, 0.8, legend=c("true", "genmix2"), col=c("blue","red"), lty=c(1,1), cex=0.8)

Άσκηση 8 (Ross 5.16)

Έστω \(X\) συνεχής τυχαία μεταβλητή με συνάρτηση κατανομής \(F(x) = 1-e^{-x} - e^{-2x} + e^{-3x}, x>0\). Δημιουργούμε δύο γεννήτριες της \(X\) ως εξής:

(α) Η \(F(x)\) μπορεί να γραφτεί ως εξής: \[ F(x) = (1-e^{-x})(1-e^{-2x}) = F_1(x) F_2(x), \] όπου \(F_1(x), F_2(x)\) είναι οι συναρτήσεις κατανομής δύο ανεξάρτητων τ.μ. \(X_1 \sim Exp(1), X_2 \sim Exp(2)\). Επομένως \(F(x) = P(X_1 \leq x, X_2 \leq x)\), δηλαδή η \(X\) ακολουθεί την κατανομή της \(\max(X_1, X_2)\). Η αντίστοιχη γεννήτρια:

genmaxexp1=function()
{
  x1=-log(runif(1,0,1))
  x2=-log(runif(1,0,1))/2
  return(max(x1,x2))
}

(β) Η συνάρτηση πυκνότητας πιθανότητας της \(X\) είναι \[ f(x) = e^{-x} + 2 e^{-2x} - 3 e^{-3x}, x>0 \] Θα δημιουργήσουμε γεννήτρια με βάση τη μέθοδο αποδοχής-απόρριψης. Παίρνουμε ως βοηθητική κατανομή την εκθετική με παράμετρο 1, δηλαδή \(g(x) = e^{-x}\). Έστω \[ h(x) = \frac{f(x)}{g(x)} = 1 + 2 e^{-x} - 3 e^{-2x} \] και \(c=\sup_{x>0} h(x)\). Η \(h(x)\) έχει μορφή

x=seq(0,10,0.2)
h=function(x) 1+2*exp(-2*x)-3*exp(-3*x)
y=h(x)
plot(y~x, type="l")

Για να προσδιορίσουμε την τιμή του \(c\) αναλυτικά έχουμε: \[ h'(x) = -2 e^{-x} + 6 e^{-2x} = e^{-2x} (6 - 2 e^x). \] Επομένως \(h'(x)=0\) για \(x=\log3\). Επίσης \(h''(x) = 2 e^{-x} - 12 e^{-2x}\) και \(h''(log3) = -2/3\), επομένως η \(h(x)\) μεγιστοποιείται στο \((0,\infty)\) για \(x=\log3\) και \[ c=h(\log3) = 4/3. \] Επομένως η γεννήτρια αποδοχής απόρριψης δημιουργεί μια παρατήρηση \(Y=y\) από την \(Exp(1)\) και την αποδέχεται με πιθανότητα \(h(y)/c = \frac{3}{4}(1 + 2 e^{-y} - 3 e^{-2y})\).

genmaxexp2=function()
{
  y=-log(runif(1,0,1))
  while (runif(1,0,1) > 3*h(y)/4 ) y=-log(runif(1,0,1))
  return(y)
}

Για να ελέγξουμε τις δύο γεννήτριες δημιουργούμε \(n=1000\) παρατηρήσεις από καθεμιά και σχεδιάζουμε τα γραφήματα της θεωρητικής και των δύο εμπειρικών συναρτήσεων κατανομής:

n=1000
X1=rep(0,n)
for (j in 1:n) X1[j]=genmaxexp1()
X2=rep(0,n)
for (j in 1:n) X2[j]=genmaxexp2()
x=seq(0,5,0.01)
nx=length(x)
trueF=1-exp(-x) - exp(-2*x) + exp(-3*x)
hatF1=rep(0,nx)
for (j in 1:nx) hatF1[j] = mean(X1<=x[j])
hatF2=rep(0,nx)
for (j in 1:nx) hatF2[j] = mean(X2<=x[j])
plot(trueF~x, type="l", col="black", ylab=expression(F(x)))
lines(hatF1~x, col="blue")
lines(hatF2~x, col="red")
legend(3, 0.3, legend=c("true", "genmax1", "genmax2"), col=c("black","blue","red"), lty=c(1,1,1), cex=0.8)

Άσκηση 9 (Ross 5.28)

Η κατανομή της \(X\) προσδιορίζεται από ένα θετικό ακέραιο \(k\), ένα διάνυσμα πιθανoτήτων σταματήματος \(a_1, \ldots, a_{k-1}\) και ένα διάνυσμα ρυθμών \(\lambda_1, \ldots, \lambda_k\). Η γεννήτρια της κατανομής της \(X\) δίνεται από την παρακάτω συνάρτηση:

genCox=function(a,l)
{
  k=length(l)
  j=0
  x=0
  stop=FALSE
  while(!(stop))
  {
    j=j+1
    x=x + rexp(1,l[j])
    if (j==k) stop=TRUE else stop=(runif(1,0,1)<=a[j])
  }
  return(x)
}

Έστω \(k=5, a_j=1/2, j=1,\ldots,4, \lambda_j=1, j=1, \ldots,5\). Δημιουργούμε 100 παρατηρήσεις από την κατανομή Cox με τις παραπάνω παραμέτρους και σχηματίζουμε ένα ιστόγραμμα και την εμπειρική συνάρτηση πυκνότητας που υπολογίζεται μέσω της density().

n=10000
a=rep(1/2,4)
l=rep(1,5)
x=rep(0,n)
for (j in 1:n) x[j]=genCox(a,l)
hist(x, probability = T)
lines(density(x), col="red")

Άσκηση 10 (Ross 5.29)

Θα δημιουργήσουμε μια γενικότερη γεννήτρια ως εξής: Έστω \(\lambda\) ο ρυθμός αφίξεων λεωφορείων. Επίσης έστω ότι ο αριθμός επιβατών σε κάθε λεωφορείο ακολουθεί διακριτή ομοιόμορφη κατανομή μεταξύ \(n_1\) και \(n_2\). Η συνάρτηση προσομοιώνει τις αφίξεις των λεωφορείων σε ένα διάστημα \([0,T]\) και επιστρέφει ένα πίνακα δύο στηλών, του οποίου η πρώτη στήλη περιέχει τις στιγμές αφίξεων και η δεύτερη τον αριθμό επιβατών σε κάθε άφιξη. Για την προσομοίωση των χρόνων αφίξεων ακολουθούμε τη μέθοδο δημιουργίας εκθετικών χρόνων μεταξύ αφίξεων. Για την προσομοίωση της διακριτής ομοιόμορφης μεταξύ \(n_1\) και \(n_2\) χρησιμοποιούμε τον τύπο \(N=\lfloor Y \rfloor\), όπου \(Y \sim U(n_1, n_2+1)\) (αποδείξτε ότι ισοδυναμεί με τη γεννήτρια αντιστροφής).

Compound_Poisson_arrivals = function(l, n1, n2, T)
{
  X=matrix(c(0,0), 1,2)
  t=0
  j=0
  stop=F
  while(!(stop))
  {
    j=j+1
    t=t+rexp(1,l)
    if (t <=T)
    {
      x=floor(runif(1,n1,n2+1))
      X=rbind(X,c(t,x))
    } else
      stop=T
  }
  X=X[-1,] #delete the first line of 0,0
  return(X[-1,])
}

Δημιουργούμε μια προσομοίωση των αφίξεων για \(\lambda=5, n_1=20, n_2=40, T=1\):

Compound_Poisson_arrivals(5,20,40,1)
##           [,1] [,2]
## [1,] 0.4752788   24
## [2,] 0.5889955   40
## [3,] 0.8527665   28
## [4,] 0.9986753   28

Άσκηση 11

Δημιουργούμε μια γεννήτρια από την πολυμεταβλητή κανονική κατανομή διάστασης \(n\) με διάνυσμα μέσης τιμής \(\mu = (\mu_1, \ldots, \mu_n)\) και πίνακα συνδιακυμάνσεων \(C\). Κάνουμε την παραγοντοποίηση Cholesky, δηλαδή υπολογίζουμε ένα κάτω τριγωνικό πίνακα \(A\) τέτοιον ώστε \(C=A A^T\) κaι εκφράζουμε το διάνυσμα \(X\) ως \(X=\mu + A Z\), όπου \(Ζ=(Z_1, \ldots, Z_n)\) ανεξάρτητες και ισόνομες \(~{\cal N}(0,1)\).

Σημειώνεται ότι η συνάρτηση chol του R υπολογίζει έναν άνω τριγωνικό πίνακα \(R\) τέτοιον ώστε \(R^T R =C\). Επομένως ο πίνακας \(A\) που απαιτείται εδώ είναι ο ανάστροφος αυτού που επιστρέφει η συνάρτηση chol.

genmvnorm=function(N, m, C)
# Generates N vectors from the multivariate normal N(m,C). Returns a matrix Nxn, where each row contains one observation from the multivariate normal distribution. 
{
  A=chol(C)
  A=t(A)
  n=length(m)
  X=matrix(rep(0,N*n), N,n)
  for (j in 1:N)
    X[j,]= A%*%rnorm(n,0,1)+m
  return(X)
}

Έστω η διανυσματική τυχαία μεταβλητή \(X\) που ακολουθεί πολυμεταβλητή κανονική κατανομή με \[ \mu = (0, -10, 20, 12), \ \ C=\left ( \begin{array}{rrrr} 6&-9&8&9 \\ -9&23&-2&-20\\ 8&-2&22&5\\ 9&-20&5&18 \end{array} \right ) . \] Δημιουργούμε \(N=1000\) παρατηρήσεις από αυτή την κατανομή σε ένα πίνακα \(X\) και υπολογίζουμε τις δειγματικές μέσες τιμές και τον δειγματικό πίνακα συνδιακυμάνσεων μεταξύ των στηλών του \(X\):

m=c(0, -10, 20, 12)
C=matrix(c(6,-9,8,9,-9, 23, -2, -20, 8, -2, 22, 5, 9, -20, 5, 18), 4,4)
N=1000
X=genmvnorm(N,m,C)
sample_means=apply(X, 2, mean)
sample_means
## [1] -0.05491126 -9.94828071 19.92559889 11.92649772
sample_C=cov(X)
sample_C
##           [,1]       [,2]      [,3]       [,4]
## [1,]  5.791296  -8.293537  8.278863   8.377767
## [2,] -8.293537  21.707731 -1.378237 -18.752705
## [3,]  8.278863  -1.378237 23.797767   4.623628
## [4,]  8.377767 -18.752705  4.623628  16.810682

Άσκηση 12

Στην άσκηση αυτή θα δημιουργήσουμε μια γεννήτρια για δεσμευμένες κατανομές από την πολυμεταβλητή κανονική. Έστω \(X \sim {\cal N}(\mu, C)\) \(n\)-διάστατη τυχαία μεταβλητή που ακολουθεί πολυμεταβλητή κανονική κατανομή με διάνυσμα μέσων τιμών \(\mu\) και πίνακα συνδιακυμάνσεων \(C\). Θεωρούμε διαμέριση του \(X\) σε ένα υποσύνολο συνιστωσών \(X_1\) και το συμπλήρωμά του \(X_2\). Μας ενδιαφέρει η δεσμευμένη κατανομή \(X_1 | X_2=x_2\).

Από τη θεωρία της πολυμεταβλητής κανονικής κατανομής γνωρίζουμε ότι δοθέντος \(X_2=x_2\), η \(X_1\) ακολουθεί πολυμεταβλητή κανονική κατανομή με μέση τιμή \[ \mu_{1|2} = \mu_1 + C_{12}C_{22}^{-1}(x_2-\mu_2) \] και πίνακα συνδιακυμάνσεων \[ C{1|2} = C_{11} - C_{12}C_{22}^{-1}C_{21}, \] όπου \(\mu_1, \mu_2\) και \(C_{11}, C_{12}, C_{21}, C_{22}\) οι αντίστοιχες διαμερίσεις των \(\mu\) και \(C\): \[ \mu = \left ( \begin{array}{l} \mu_1\\ \mu_2 \end{array} \right ), \] \[ C = \left ( \begin{array}{ll} C_{11}&C_{12}\\ C_{21}&C_{22} \end{array} \right ). \]

Παρακάτω δημιουργούμε μια συνάρτηση που υπολογίζει τις παραμέτρους \(\mu_{1|2}, C_{1|2}\) και δημιουργεί \(N\) παρατηρήσεις από τη δεσμευμένη κατανομή. Οι συνιστώσες του διανύσματος \(X\) που ορίζουν το \(X_1\) δίνονται με τη μορφή ενός διανύσματος συνιστωσών \(f \subset \{1,\ldots, n\}\). Ορίζοντας \(g = \{1,\ldots, n\} - f\) το συμπλήρωμα του \(f\), οι διαμερίσεις των \(\mu, C\) υπολογίζονται άμεσα ως m[f], m[g] και C[f,f], C[f,g], C[g,f], C[g,g].

gencondmvnorm=function(N, m, C, f, x)
{
  # Generate N observations from the conditional distribution X[f] | X[-f]=x, where X is multivariate normal with mean m and covariance matrix C
  n=length(m)
  K=1:n
  g=K[-f]
  m1=m[f]
  m2=m[g]
  C11=C[f,f]
  C12=C[f,g]
  C21=C[g,f]
  C22=C[g,g]
  C22inv=solve(C22)
  m12=m1+C12%*%C22inv%*%(x-m2)
  S12=C11=C12%*%C22inv%*%C21
  X=genmvnorm(N, m12, S12)
  return(X)
}

Θεωρούμε τώρα το παράδειγμα της προηγούμενης άσκησης και έστω η δεσμευμένη κατανομή των \((X_1, X_3)\) δοθέντος ότι \(X_2=7, X_4=3\). Μπορούμε να προσομοιώσουμε μια παρατήρηση \((x_1, x_3)\) από αυτή τη δεσμευμένη κατανομή

x=gencondmvnorm(1,m,C, c(1,3), c(7,3))
x
##          [,1]     [,2]
## [1,] 6.305656 49.98046