Προσομοίωση

1. Γεννήτριες ψευδοτυχαίων ομοιόμορφων αριθμών - Linear Congruential Generator

Η πιο κάτω συνάρτηση υλοποιεί τον αλγόριθμο Linear Congruential Generator.

lcg <- function(m,a,c,X0, N)
{ # Integer m>1 is the modulus
  # Integer 0<a<m is the multuplier
  # Integer 0<=c<m is the increment
  # X0 natural number or zero is the seed
  # N is the number of uniform samples required
  # returns sequence X1,...,XN of numbers which for good values of m,a,c 
  # behave like uniformly distributed independent random numbers
  x <- c()
  x[1] <- X0
  for(i in 1:N)
  {
    x[i+1] <- (a*x[i] + c)%%m
  }
  return(x[-1])
}

Για παράδειγμα για \(m=8, a=5, c=1, X_0=0, N=10\), παίρνουμε την εξής ακολουθία που συμφωνεί με τον υπολογισμό στη θεωρία.

x <- lcg(8,5,1,0,1000)
x

##    [1] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##   [35] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##   [69] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [103] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [137] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [171] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [205] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [239] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [273] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [307] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [341] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [375] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [409] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [443] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [477] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [511] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [545] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [579] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [613] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [647] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [681] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [715] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [749] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [783] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [817] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [851] 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4
##  [885] 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2
##  [919] 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0
##  [953] 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6 7 4 5 2 3 0 1 6
##  [987] 7 4 5 2 3 0 1 6 7 4 5 2 3 0

y <- lcg(2^32, 1103515245, 12345, 0, 100)
y

##   [1]      12345 3554416254 2802067456 1358247936 2138638336 1459132416
##   [7] 1445521408 2518349824 4044081152 3908327424 3000608768 1742518272
##  [13] 3372910592 4276885504  165046272 4093278272 1942545408 3657516032
##  [19] 3443663872 1679918080 4034749440 4057711616 2286593024 2646606848
##  [25] 3287933952 3973989376  727393280  314999808 2461951040 1193041920
##  [31]  743772160  722993152 2839031808 1267507200 3282063360   60964864
##  [37] 2073391160  706873344 2129500160 4286179328 1536686080  491046912
##  [43]  456742976 2446214016 2724457984 2869986816 2195537408  353867264
##  [49]  784272960 1664139520 1552964864 1405995264 3621569792 2332430336
##  [55] 1786052608 4191236096  984719360 3697364992 1540308992 3166597120
##  [61]   53862400  137400376  200189984 4103570912 3129896960 1959124992
##  [67] 1115721728 3701309440 3960451072 1783504896 3506688000  134201344
##  [73] 3563384896 1101921280 3789323264 1084816384 3971709952 3550215168
##  [79] 2285009920 1195834368 3791748096 2523528192 1038646272  432339968
##  [85] 3852958784  699111424 2522867712 2944792576 1570424832 3017234432
##  [91] 1065396224  754432000 4018767872  352262144  971860032 3780711168
##  [97] 3367454720 4261009408 1395730432 1444820992

Για να πάρουμε δείγμα από την ομοιόμορφη συνεχή στο διάστημα \([0,1]\), επιστρέφουμε \(U_n=\frac{X_n+1}{m+1}\).

lcg01 <- function(m, a, c, X0, N)
{ # Integer m>1 is the modulus
  # Integer 0<a<m is the multiplier
  # Integer 0<=c<m is the increment
  # X0 natural number or zero is the seed
  # N is the number of uniform samples required
  # returns sequence X1,...,XN of numbers which for good values of m,a,c 
  # behave like uniformly distributed independent random numbers
  x <- c()
  x[1] <- X0
  for(i in 1:N)
  {
    x[i+1] <- (a*x[i] + c) %% m
  }
  return((x[-1] + 1) / (m + 1))
}

x <- lcg01(8,5,1,0,1000)
x[1:30]

##  [1] 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667 0.3333333 0.4444444
##  [8] 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667 0.3333333
## [15] 0.4444444 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556 0.6666667
## [22] 0.3333333 0.4444444 0.1111111 0.2222222 0.7777778 0.8888889 0.5555556
## [29] 0.6666667 0.3333333

hist(x)

plot.ts(x, main="LCG Path m=8, a=5, c=1, x0=0")

y <- lcg01(2^32, 1103515245, 12345, 0, 1000)
y[1:30]

##  [1] 2.874527e-06 8.275770e-01 6.524072e-01 3.162417e-01 4.979405e-01
##  [6] 3.397307e-01 3.365617e-01 5.863490e-01 9.415860e-01 9.099784e-01
## [11] 6.986337e-01 4.057117e-01 7.853169e-01 9.957900e-01 3.842783e-02
## [16] 9.530406e-01 4.522841e-01 8.515818e-01 8.017905e-01 3.911364e-01
## [21] 9.394133e-01 9.447596e-01 5.323889e-01 6.162112e-01 7.655318e-01
## [26] 9.252665e-01 1.693594e-01 7.334161e-02 5.732176e-01 2.777767e-01

hist(y)

plot.ts(y, main="LCG Path m=2^32, a=1103515245, c=12345, x0=0")

Στη 1η περίπτωση επειδή η περίοδος είναι πολύ μικρή δεν φαίνονται να είναι τυχαίοι οι αριθμοί αν και το ιστόγραμμα δείχνει να μοιράζονται ομοίομορφα στο διάστημα \([0,1]\). Στη 2η περίπτωση φαίνονται όντως τυχαίοι.

ks.test(x, "punif")

## Warning in ks.test(x, "punif"): ties should not be present for the
## Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.11111, p-value = 3.782e-11
## alternative hypothesis: two-sided

ks.test(y, "punif")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.02262, p-value = 0.6856
## alternative hypothesis: two-sided

Βλέπουμε ότι με τις κατάλληλες επιλογές των m, a, c η ακολουθία που επιστρέφει ο αλγόριθμος LCG μπορεί να ξεγελάσει το Kolmogorov-Smirnov test και άρα αποτελεί υποκατάστατο για τυχαίο ομοιόμορφο δείγμα στο \([0,1]\).

Πιο κάτω διερευνούμε κατά πόσο τα στοιχεία της ακολουθίας φαίνονται να είναι ανεξάρτητα μεταξύ τους. Συγκεκριμένα παρουσιάζουμε διαγράμματα διασποράς των \(X_i\) ως προς τα \(X_{i+1}\) για διάφορες επιλογές των παραμέτρων m, a, c καθώς και όταν χρησιμοποιούμε την γεννήτρια ομοιόμορφων αριθμών της R.

par(mfrow=c(2,2))
u <- runif(1000)
x <- lcg01(81, 1, 8, 0, 1000)
y <- lcg01(1024, 401, 101, 0, 1000)
z <- lcg01(2^32, 1664525, 1013904223, 0, 1000)

plot(u[-1000], u[-2], xlab="Xn", ylab="X(n+1)")
plot(x[-1000], x[-2], xlab="Xn", ylab="X(n+1)")
plot(y[-1000], y[-2], xlab="Xn", ylab="X(n+1)")
plot(z[-1000], z[-2], xlab="Xn", ylab="X(n+1)")

Παρατηρούμε ότι για κακές επιλογές των παραμέτρων το αποτέλεσμα του LCG δεν μοιάζει να είναι ανεξάρτητο δείγμα. Στο πάνω δεξιά γράφημα φαίνεται ξεκάθαρα το ότι κάθε \(X_n\) μπορεί να ακολουθείται από ακριβώς ένα \(Χ_{n+1}\). Αν και το ίδιο ισχύει για τα δύο γραφήματα στο κάτω μέρος, η εξάρτηση είναι πολύ πιο πολύπλοκη και συγκεκριμένα στο κάτω δεξιά γράφημα δεν μπορεί να διακρηθεί η εξάρτηση.

Πιο κάτω εμφανίζουμε και τα Autocorrelation plots:

acf(u, lag.max=1000)

acf(x, lag.max=1000)

acf(y, lag.max=1000)

acf(z, lag.max=1000)

Τέλος, διερευνούμε την επίδραση της αρχικής συνθήκης \(Χ_0\).

lcg01(2^32, 1664525, 1013904223, X0=0, 10)

##  [1] 0.2360680 0.2785669 0.8195338 0.6678669 0.3840774 0.6218075 0.3437259
##  [8] 0.6400356 0.5077781 0.5817975

lcg01(2^32, 1664525, 1013904223, X0=0, 10)

##  [1] 0.2360680 0.2785669 0.8195338 0.6678669 0.3840774 0.6218075 0.3437259
##  [8] 0.6400356 0.5077781 0.5817975

lcg01(2^32, 1664525, 1013904223, X0=Sys.time(), 10)

##  [1] 0.75848959 0.11837989 0.51934885 0.37326421 0.84494024 0.39026531
##  [7] 0.60150061 0.03835695 0.33988502 0.35351772

lcg01(2^32, 1664525, 1013904223, X0=Sys.time(), 10)

##  [1] 0.7584900 0.7964029 0.7053922 0.2250348 0.2627183 0.3544751 0.8884294
##  [8] 0.2637895 0.4339731 0.3925671

runif(10)

##  [1] 0.9564193 0.4724768 0.4912508 0.0817008 0.8336688 0.6610708 0.8902970
##  [8] 0.9221510 0.4044853 0.4541261

runif(10)

##  [1] 0.46460075 0.98192125 0.03131526 0.21535477 0.26002716 0.12271118
##  [7] 0.11684848 0.30052653 0.96816678 0.17911271

set.seed(1)
runif(10)

##  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
##  [7] 0.94467527 0.66079779 0.62911404 0.06178627

runif(10)

##  [1] 0.2059746 0.1765568 0.6870228 0.3841037 0.7698414 0.4976992 0.7176185
##  [8] 0.9919061 0.3800352 0.7774452

set.seed(1)
runif(10)

##  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
##  [7] 0.94467527 0.66079779 0.62911404 0.06178627

runif(10)

##  [1] 0.2059746 0.1765568 0.6870228 0.3841037 0.7698414 0.4976992 0.7176185
##  [8] 0.9919061 0.3800352 0.7774452

2. Προσομοίωση Διακριτών Τυχαίων Μεταβλητών

2.1. Μέθοδος Αντρίστροφου Μετασχηματισμού

Έστω θέλουμε να προσομοιώσουμε από τυχαία μεταβλητή \(Χ\) που παίρνει 3 τιμές \(x_1, x_2, x_3\) (χωρίς βλάβη της γενικότητας \(1,2,3\)) με πιθανότητες \(p_1, p_2, p_3\). Η προσομοίωση μπορεί να γίνει με την εξής συνάρτηση:

my.discrete.sim <- function(n, p)
{
  # n is the number of samples required
  # p einai i sunartisi mazas pithanotitas
  # eleghoume oti p einai s.m.p tritimis
  if(length(p) != 3 | min(p) < 0 | sum(p) != 1)
  {
    stop('p vector of probabilities of length 3')
  }
  sim.vector <- c()
  for(j in 1:n)
  {
    u <- runif(1)
    ifelse(u < p[1], sim.vector[j] <- 1, ifelse(u < p[1] + p[2], sim.vector[j] <- 2, sim.vector[j] <- 3))
  }
  return(sim.vector)
}

Η συνάρτηση παίρνει ως ορίσματα το επιθυμητό πλήθος δειγμάτων \(n\) και το διάνυσμα πιθανοτήτων \(p\) και επιστρέφει \(n\) δείγματα από την αντίστοιχη τρίτιμη κατανομή.

x <- my.discrete.sim(1000, c(0.3, 0.5, 0.2))
hist(x)

barplot(table(x))

Ο έλεγχος Kolmogorov-Smirnov είναι ακατάλληλος για διακριτές κατανομές (ειδικά αν έχουν μόνο λίγες τιμές και άρα δεν μοιάζουν καθόλου με συνεχή κατανομή). Ένας τρόπος για να ελέγξουμε με στατιστικό έλεγχο κατά πόσο το δείγμα που πήραμε ακολουθεί την επιθυμητή κατανομή, είναι με το Chi-squared goodness of fit test. Αυτός ο έλεγχος συγκρίνει τις συχνότητες με τις οποίες εμφανίζονται οι διάφορες τιμές (εδώ 3) στο δείγμα, με κάποια συνάρτηση μάζας πιθανότητας. Η μηδενική υπόθεση είναι ότι το δείγμα προέρχεται από την επιθυμητή συνάρτηση μάζας πιθανότητας.

freq <- table(x) # dinei poses fores emfanistike kathe timi
freq

## x
##   1   2   3 
## 294 506 200

chisq.test(freq, p=c(0.3, 0.5, 0.2))

## 
##  Chi-squared test for given probabilities
## 
## data:  freq
## X-squared = 0.192, df = 2, p-value = 0.9085

chisq.test(freq, p=c(0.2, 0.7, 0.1))

## 
##  Chi-squared test for given probabilities
## 
## data:  freq
## X-squared = 197.95, df = 2, p-value < 2.2e-16

Προσομοίωση Διακριτής Ομοιόμορφης: H πιο κάτω συνάρτηση επιστρέφει \(n\) δείγματα από τη διακριτή ομοιόμορφη στο \(\{1,\dots,k\}\).

my.discrete.unif <- function(n, k)
{
  sim.vector <- c()
  for (j in 1:n) {
  u <- runif(1)
  sim.vector[j] <- floor(k*u) + 1
  }
  return(sim.vector)
}

Προγραμματιστικά καλύτερη υλοποίηση:

my.discrete.unif1 <- function(n,k)
{
  return(ceiling(k*runif(n)))
}

x <- my.discrete.unif1(100, 5)
table(x)

## x
##  1  2  3  4  5 
## 21 20 23 20 16

barplot(table(x))

chisq.test(table(x), p = rep(0.2,5))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(x)
## X-squared = 1.3, df = 4, p-value = 0.8614

Προσομοίωση Poisson: Η πιο κάτω συνάρτηση επιστρέφει \(n\) δείγματα από την κατανομή Poisson με ρυθμό lambda.

my.poisson <- function(n,lambda)
{
  u <- runif(n)
  sim.vector <- c()
  for(j in 1:n){
     i <- 0
     p <- exp(-lambda)
     F <- p
     while(u[j] > F){ 
       p <- lambda*p/(i+1)
       F <- F+p
       i<-i+1
    }
    sim.vector[j] <- i
  }
return(sim.vector)
}

x <- my.poisson(100, 5)
barplot(table(x))

Για να κάνουμε τον έλεγχο chi-square θα συγκρίνουμε τις συχνότητες με τις οποίες εμφανίστηκαν οι διάφορες τιμές, με ένα διάνυσμα πιθανοτήτων που αντιστοιχούν στη συνάρτηση μάζας πιθανότητας της Poisson για αυτές τις τιμές. Επειδή στο δείγμα έχουμε άνω όριο στις τιμές, ενώ η Poisson παίρνει ως τιμές όλους τους φυσικούς αριθμούς, πρέπει να δώσουμε και την πιθανότητα να έχουμε μεγαλύτερη τιμή σύμφωνα με την Poisson από το μέγιστο του δείγματος. Επίσης στο δείγμα πιθανόν κάποιες μικρές τιμές (πχ το 0) ή και κάποιες ενδιάμεσες τιμές να μην εμφανίζονται, συνεπώς σε τέτοιες περιπτώσεις πρέπει να προσθέσουμε την πληροφορία στο table(x) ότι δεν εμφανίζονται αυτές οι τιμές.

values=as.integer(names(table(x)))
probs <- dpois(values, lambda = 5) # i syn. mazas. pithanotitas stis times pou emfanizontai
rem <- 1 - sum(probs) # i syn. mazas pithanotitas stis times pou den emfanizontai
chisq.test(c(table(x),0), p = c(probs, rem)) # sygkirnw tis syxnotites pou paratirw, me tis antistoixes thewritikes pithanotites 

## Warning in chisq.test(c(table(x), 0), p = c(probs, rem)): Chi-squared
## approximation may be incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  c(table(x), 0)
## X-squared = 15.881, df = 10, p-value = 0.1031

2.2. Mέθοδος Απόρριψης

Θέλουμε να προσομοιώσουμε \(n=100\) δείγματα από διακριτή κατανομή με 5 τιμές, έστω \(\{1,2,3,4,5\}\), με αντίστοιχο διάνυσμα πιθανοτήτων \(p=(0.05, 0.25, 0.45, 0.15, 0.1)\). Στην πιο κάτω συνάρτηση χρησιμοποιούμε τη μέθοδο απόρριψης με χρήση της διακριτής ομοιόμορφης, \(q=(0.2,0.2,0.2,0.2,0.2)\). Είδαμε στη θεωρία ότι \(M=2.25\) και άρα περιμένουμε ότι η μέθοδος θα χρειάζεται κατά μέσο όρο \(nM\) επαναλήψεις για να παράξει τα \(n\) δείγματα.

my.disc.reject <- function(n){

  p <- c(.05,.25,.45,.15,.10)
  sim.vector=c()
  reps=0

  for(j in 1:n){
    u1 <- runif(1)
    y <- floor(5*u1) + 1 # omoiomorfi sto {1,...,5}
    u2 <- runif(1)
    reps <- reps + 1 # metroume toso apotiximenes oso 
    # kai epitiximenes epanalipseis
     
    while(u2 >= p[y]/0.45){ #M=2.25, q_Y=0.2 gia oles tis times tou Y
      u1 <- runif(1)
      y <- floor(5*u1) + 1
      u2 <- runif(1)
      reps <- reps + 1
    }
 
    sim.vector[j] <- y
  }

return(list(sim.vector=sim.vector, no.reps=reps))
}

sim <- my.disc.reject(100)
x <- sim$sim.vector
r <- sim$no.reps
barplot(table(x))

r # Mn=2.25*n

## [1] 253

chisq.test(table(x), p = c(.05,.25,.45,.15,.10))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(x)
## X-squared = 1.0667, df = 4, p-value = 0.8995

3. Προσομοίωση Συνεχών Τυχαίων Μεταβλητών

3.1. Μέθοδος Αντρίστροφου Μετασχηματισμού

\(\textbf{Προσομοίωση Εκθετικής Κατανομής:}\) Η πιο κάτω συνάρτηση παίρνει ρυθμό \(\lambda\) και επιθυμητό πλήθος δειγμάτων \(n\) και επιστρέφει, n ανεξάρτητα δείγματα από την κατανομή Exp(\(\lambda\)).

my.exponential <- function(n, lambda)
{
  if(lambda<=0){
    stop("lambda needs to be positive")
  }
  return(-log(runif(n))/lambda)
}


sim.vector<- my.exponential(100, 2) #dokimaste n=1000,10000 i alla lambda
sim.vector[1:10]

##  [1] 0.48546587 0.01751487 0.22227204 0.10775182 0.04220074 1.41831879
##  [7] 0.36944599 0.18056902 1.55420968 0.13929225

Πιο κάτω ελέγχουμε με διάφορους τρόπους (γραφικούς και με στατιστικό έλεγχο), ότι το δείγμα που πήραμε ακολουθεί πράγματι εκθετική κατανομή με παράμετρο \(\lambda=2\).

hist(sim.vector)

dens=density(sim.vector)
x=dens$x

plot(dens, main="Density Estimation", ylim=c(0,2))
lines(x,dexp(x,rate=2), col="blue") 

legend(1.5,1,legend=c("sample", "εκθετική"),
       col=c("black", "blue"), lty=1, cex=0.8)

qqplot(qexp(ppoints(10000), rate=2), sim.vector)
#H ppoints(m) dinei m isapexonta simeia metaxy 0 kai 1 
#ara h qexp(ppoints(m),...) dinei ta thewrtitika quantiles (peripou)
#pou xreiazomaste gia to qqplot
abline(0,1)

# i qqline opws einai den epitrepei na perasoume parametrous, opws rate=2
# Mporoume na allaksoume ton orismo tis synartisis qqline.exp gia na 
# tin kanoume na dexetai kai rythmo, deste paragrafo 9.5.2 pio katw
# https://www.stat.berkeley.edu/~rice/Stat135/The.R.Project.R.D.9.8.8.pdf

ks.test(sim.vector, "pexp", rate=2)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.vector
## D = 0.10508, p-value = 0.2194
## alternative hypothesis: two-sided

\(\textbf{Προσομοίωση}\) Gamma(\(\alpha, \beta\)) \(\textbf{με}\) \(\alpha\) \(\textbf{ακέραιο θετικό:}\)

H πιο κάτω συνάρτηση προσομοιώνει \(n\) δείγματα από κατανομή Gamma(\(\alpha, \beta\)), με \(\alpha\) ακέραιο θετικό και \(\beta>0\) με χρήση της μεθόδου αντίστροφου μετασχηματισμού.

my.intgamma <- function(n, a, b)
{
  if(floor(a)!=a | a<=0 | b<=0)
  { 
    stop('a needs to be positive integer, b>0')
  }
  m <- matrix(runif(n*a), nrow=a) #se kathe mia apo tis n stiles tou,
  #aftos o pinakas ehei tis a omoiomorfes t.m. pou xreiazomai gia na 
  #paraksw ena deigma apo ti Gamma(a,b)
  s=apply(m,2,prod) #ypologizw ginomeno kathe stilis, efarmozw tin
  #synartisi product, kata stiles (to 2 auto leei), ston pinaka m
  return(-log(s)/b)
}

sim.vector <- my.intgamma(1000, 2, 3)
sim.vector[1:10]

##  [1] 0.4709756 0.1264983 0.4767356 0.3561808 0.3892885 0.4255618 0.8783046
##  [8] 0.1327797 0.0152201 1.0597637

Πιο κάτω ελέγχουμε με διάφορους τρόπους (γραφικούς και με στατιστικό έλεγχο), ότι το δείγμα που πήραμε ακολουθεί πράγματι κατανομή Gamma(2,3).

hist(sim.vector)

dens=density(sim.vector)
x=dens$x

plot(dens, main="Density Estimation", ylim=c(0,1.2))
lines(x,dgamma(x,shape=2, rate=3), col="blue") 

legend(1.5,1,legend=c("sample", "Gamma(2,3)"),
       col=c("black", "blue"), lty=1, cex=0.8)

qqplot(qgamma(ppoints(10000), shape=2, rate=3), sim.vector)
abline(0,1)

ks.test(sim.vector, "pgamma", shape=2, rate=3)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.vector
## D = 0.039389, p-value = 0.08982
## alternative hypothesis: two-sided

3.2. Mέθοδος Απόρριψης

Πιο κάτω υλοποιούμε τη μέθοδο απόρριψης για προσομοίωση από συνεχή τυχαία μεταβλητή, με πεδίο τιμών το \((0,1)\) και συνάρτηση πυκνότητας πιθανότητας \(f_X(x)=20x(1-x)^3, x\in(0,1)\) (είναι η κατανομή Beta(2,4)). Θα προτείνουμε δείγματα από την ομοιόμορφη, \(g(x)=1, x\in(0,1)\). Είδαμε στη θεωρία ότι \(Μ=135/64\).

my.cts.reject <- function(n){

  sim.vector=c()
  reps=0
  M=135/64

  for(j in 1:n){
    y <- runif(1)
    u <- runif(1)
    reps <- reps + 1 # metroume toso apotiximenes oso 
    # kai epitiximenes epanalipseis
     
    while(u > 20*y*(1-y)^3/M){ #g(y)=1
      y <- runif(1)
      u <- runif(1)
      reps <- reps + 1
    }
 
    sim.vector[j] <- y
  }

return(list(sim.vector=sim.vector, no.reps=reps))
}

Πιο κάτω ελέγχουμε με διάφορους τρόπους (γραφικούς και με στατιστικό έλεγχο), ότι το δείγμα που πήραμε ακολουθεί πράγματι κατανομή Beta(2,4).

n=10000
sim<-my.cts.reject(n)
sim.vector<-sim$sim.vector
hist(sim.vector)

dens=density(sim.vector)
x=dens$x

plot(dens, main=" Density Estimation", ylim=c(0,2.2))
lines(x,dbeta(x,shape1=2, shape2=4), col="blue") 

legend(0.7,1,legend=c("sample", "Beta(2,4)"),
       col=c("black", "blue"), lty=1, cex=0.8)

qqplot(qbeta(ppoints(n), shape1=2, shape2=4), sim.vector)
abline(0,1)

ks.test(sim.vector, "pbeta", shape1=2, shape2=4)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.vector
## D = 0.0064871, p-value = 0.794
## alternative hypothesis: two-sided

r<-sim$no.reps
r #Mn peripou 2.11*n

## [1] 21028

3.3. Mέθοδος Box-Muller για Προσομοίωση από Κανονική Κατανομή

Η πιο κάτω συνάρτηση υλοποιεί τον αλγρόριθμο Box-Muller. Παίρνει ως όρισμα τον επιθυμητό αριθμό δειγμάτων \(n\), τη μέση τιμή \(m\) και την τυπική απόκλιση \(sd\) και επιστρέφει \(n\) ανεξάρτητα δείγματα από την \(N(m, sd^2)\). Αν δεν δοθούν τιμές για τη μέση τιμή και την τυπική απόκλιση από τον χρήστη, τότε αυτές παίρνουν την τιμή 0 και 1 αντίστοιχα.

my.box.muller <- function(n, mean = 0, sd = 1)
{
  # Produces n independent samples from N(m,sd^2)
  # n = number of samples
  # m = mean, takes default value 0
  # s = standard deviation, takes default value 1
  k <- ceiling(n/2) # B-M paragei kathe fora 2 aneksartites N(0,1)
  theta <- 2 * pi * runif(k)
  R <- sqrt(- 2 * log(runif(k)))
  return(mean + sd * c(R * cos(theta), R * sin(theta))[1:n]) 
  # Kratoume ta n prwta, afou an 
  # n perittos tha eixame n+1 deigmata
}

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

n=1000
sim.vector<-my.box.muller(n)
hist(sim.vector)

dens=density(sim.vector)
x=dens$x

plot(dens, main="Density Estimation", ylim=c(0,0.45))
lines(x,dnorm(x), col="blue") 

legend(2.2,0.3,legend=c("sample", "N(0,1))"),
       col=c("black", "blue"), lty=1, cex=0.8)

qqnorm(sim.vector)
qqline(sim.vector)

ks.test(sim.vector, "pnorm")

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.vector
## D = 0.038701, p-value = 0.1
## alternative hypothesis: two-sided

n=1000
sim.vector<-my.box.muller(n, 2, 2)
hist(sim.vector)

dens=density(sim.vector)
x=dens$x

plot(dens, main="Density Estimation", ylim=c(0,0.22))
lines(x,dnorm(x, 2 ,2), col="blue") 

legend(6,0.15,legend=c("sample", "N(2,4)"),
       col=c("black", "blue"), lty=1, cex=0.8)

qqplot(qnorm(ppoints(n), mean=2, sd=2), sim.vector)
abline(0,1)

ks.test(sim.vector, "pnorm", mean=2, sd=2)

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sim.vector
## D = 0.01957, p-value = 0.8384
## alternative hypothesis: two-sided

Πιο κάτω συγκρίνουμε το χρόνο εκτέλεσης της συνάρτησης που γράψαμε (αλγόριθμος Box-Muller) με τον αντίστοιχο χρόνο εκτέλεσης της rnorm που μας παρέχει η R. Βλέπουμε ότι έχουν παρόμοιο κόστος.

n <- 10000000
tic <- Sys.time()
sim.vector <- my.box.muller(n)
Sys.time()-tic

## Time difference of 0.985533 secs

tic <- Sys.time()
sim.vector <- rnorm(n)
Sys.time()-tic

## Time difference of 0.719965 secs

4. Mίξη Κατανομών

Η πιο κάτω συνάρτηση παράγει \(n\) δείγματα από την μικτή κατανομή \[a_1 P_1 + a_2 P_2 + a_3 P_3,\] όπου \(P_1=Ν(m_1,s_1^2), P_2=N(m_2, s_2^2), P_3=N(m_3, s_3^2)\) (ανεξάρτητες) και \(a_1+a_2+a_3=1, a_i\in[0,1]\). Παίρνει ως ορίσματα το \(n\) και τα διανύσματα των βαρών \(a=(a_1, a_2, a_3)\), των μέσων \(m=(m_1, m_2, m_3)\) και των τυπικών διασπορών \(sd=(s_1, s_2, s_3)\).

Προσοχή στο ότι η κατανομή που προκύπτει δεν είναι η κατανομή μιας τυχαίας μεταβλητής \[a_1 X_1+a_2 X_2+a_3X_3,\] όπου \(Χ_i \sim P_i\) ανεξάρτητες, η οποία είναι κανονική.

my.mixture<-function(n, a, m, sd)
{
  sim.vector<-c()
  u<-runif(n)
  for(j in 1:n)
  {
    ifelse(u[j]<a[1], sim.vector[j]<-rnorm(1, m[1],sd[1]),
           ifelse(u[j]<a[1]+a[2], sim.vector[j]<-rnorm(1, m[2],sd[2]), 
           sim.vector[j]<-rnorm(1, m[3],sd[3])))
  }
  return(sim.vector)
}

n=10000
a=c(0.1, 0.7, 0.2)
m=c(1, 2, 5)
sd=c(0.01, 0.6, 0.1)

sim<-my.mixture(n, a, m, sd)
hist(sim, breaks=100, main="Mixture of 3 Gaussians - Histogram")

plot(density(sim), main="Mixture of 3 Gaussians - Density Plot")

5. Προσμοίωση Διαδικασίας Poisson

H πιο κάτω συνάρτηση παράγει 1 πραγμάτωση διαδικασίας Poisson, \(\{N(t), t\in[0,T]\}\), με σταθερό ρυθμό \(\lambda\). Παίρνει ως όρισμα τον ρυθμό lambda και το χρόνο \(T\) και επιστρέφει μια λίστα που περιέχει ως μεταβλητές τον τελικό αριθμό γεγονότων \(Ν(Τ)\) και τους χρόνους στους οποίους είχαμε γεγονότα (υποσύνολο του \([0,T]\)).

my.hpoisproc<-function(lambda, T)
{
  t <- 0
  N <- 0
  s <- c()
  while(t <= T){
    u <- runif(1)
    t <- t-log(u)/lambda
    N <- N+1
    s[N] <- t
  }
  N <- N - 1
  s <- s[1:N]
  return(list( N=N, times=s))
}

Πιο κάτω κάνουμε γράφημα μιας διαδικασίας Poisson με ρυθμό \(\lambda\) στο διάστημα \([0,T]\). Εφόσον κάθε φορά που έχουμε άφιξη, η διαδικασία \(N(t)\) ανεβαίνει κατά 1, ενώ σε όλους τους άλλους χρόνους είναι σταθερή, αρκεί να κάνουμε γράφημα της ακολουθίας \((1,2,\dots, N)\) ως προς τους χρόνους αφίξεων. Χρησιμοποιούμε το όρισμα type=“s” στην εντολή plot, για να έχουμε σκαλοπάτια .

T<-10
sim<-my.hpoisproc(1/2, T)
times <- sim$times
N <- sim$N
plot(c(0,times,T),c(0:N,N),type="s", main="Simulated Poisson Process",xlab="t",ylab="N(t)")