Άσκηση 1.

Έστω \(X\) συνεχής τυχαία μεταβλητή με συνάρτηση πυκνότητας πιθανότητας \[ f(x) = x+\frac{1}{2}, \ \ 0 \leq x \leq 1. \]

Να περιγραφεί γεννήτρια τυχαίων αριθμών από αυτή την κατανομή χρησιμοποιώντας τη μέθοδο αντίστροφου μετασχηματισμού.

Απάντηση

Η τυχαία μεταβλητή παίρνει τιμές στο διάστημα \([0,1]\). Η αθροιστική συνάρτηση κατανομής \(F(x)\) για \(0 \leq x \leq1\) είναι \[ F(x) = \int_0^x f(y) dy = \frac{x^2+x}{2} . \] Σύμφωνα με τη μέθοδο του αντίστροφου μετασχηματισμού, αν \(U \sim U(0,1\), τότε αν λύσουμε την εξίσωση \(F(X)=U\) ως προς \(X\), η λύση \(X\) που προκύπτει έχει συνάρτηση κατανομής την \(F\). Επομένως η γεννήτρια αντίστροφου μετασχηματισμού εδώ ορίζεται ως εξής: Η εξίσωση \(F(X) = U\) γράφεται \[ \frac{X^2+X}{2} = U \Leftrightarrow X^2 + X - 2 U = 0. \] Η δευτεροβάθμια εξίσωση έχει δύο λύσεις : \[ X_1 = \frac{-1 - \sqrt{1+8U}}{2}, \ \ X_2 = \frac{-1 + \sqrt{1+8U}}{2}. \] Επειδή θέλουμε \(0 \leq X \leq 1\), η πρώτη λύση απορρίπτεται επομένως η γεννήτρια αντίστροφου μετασχηματισμού για την κατανομή \(F\) είναι: \[ X = \frac{-1 + \sqrt{1+8U}}{2}, \] όπου \(U \sim U(0,1)\).

Η γεννήτρια μπορεί να υλοποιηθεί με την παρακάτω συνάρτηση, η οποία δημιουργεί \(n\) παρατηρήσεις

generate_F=function(n){
    u=runif(n,0,1)
    x=(-1 + sqrt(1+8*u))/2
    return(x)
}

Δημιουργούμε 2000 παρατηρήσεις από την κατανομή και παίρνουμε το ιστόγραμμα και την εκτιμημένη συνάρτηση πυκνότητας:

set.seed(229)
X=generate_F(2000)
hist(X, probability = TRUE, 
     col = "lightblue", border = "white", 
     main = "Histogram with Density Curve", 
     xlab = "x")
lines(density(X), col="red")

Βλέπουμε ότι η συνάρτηση πυκνότητας προσεγγίζει αρκετά καλά τη δοσμένη \(f(x)\) στο διάστημα \([0,1]\). Για έναν πιο αυστηρό έλεγχο, μπορούμε να χρησιμοποιήσουμε τον έλεγχο Kolmogorov-Smirnov, όπου πρώτα πρέπει να οριστεί η συνάρτηση κατανομής της \(X\):

F=function(x)
  return((x^2+x)/2)
ks.test(X, F)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  X
## D = 0.017494, p-value = 0.5731
## alternative hypothesis: two-sided

Άσκηση 2.

Μια ασφαλιστική εταιρεία έχει 500 ασφαλισμένους καθένας από τους οποίους έχει πιθανότητα 5% να ζητήσει μια αποζημίωση στη διάρκεια της επόμενης χρονιάς. Υποθέστε ότι το ποσό μιας αποζημίωσης είναι τυχαία μεταβλητή που ακολουθεί εκθετική κατανομή με μέση τιμή 1200 ευρώ. Να εκτιμήσετε μέσω προσομοίωσης την πιθανότητα το συνολικό ποσό αποζημιώσεων που θα πρέπει να πληρώσει η εταιρεία την επόμενη χρονιά να υπερβαίνει τις 35000 ευρώ.

Απάντηση

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

Επίσης, έστω \(N\) ο αριθμός ασφαλισμένων που θα ζητήσουν αποζημίωση το επόμενο έτος. Η \(N\) είναι τυχαία μεταβλητή με \(N \sim B(n, p)\). Αν \(N=k\), τότε το συνολικό ποσό αποζημίωσης είναι \(S = X_1 + \cdots + X_k\), όπου \(X_1, \ldots, X_k\) ανεξάρτητες και ισόνομες τυχαίες μεταβλητές με κατανομή \(X \sim Exp(1/m)\).

Ζητείται να εκτιμηθεί η πιθανότητα \(\theta = P(S>D)\).

Για την εκτίμηση του \(\theta\) θα εφαρμόσουμε προσομοίωση Monte Carlo, όπου σε κάθε επανάληψη προσομοιώνεται πρώτα η τιμή της \(N\) και με βάση αυτή την τιμή, οι τιμές των ποσών αποζημίωσης. Έπειτα υπολογίζεται το συνολικό ποσό αποζημιώσεων \(S\) και το αποτέλεσμα της συνάρτησης είναι \(Z=0\) ή 1, αν στη συγκεκριμένη επανάληψη το συνολικό ποσό έχει υπερβεί την τιμή \(D\) ή όχι.

Αν \(Z_1, \ldots, Z_M\) είναι τα αποτελέσματα \(M\) επαναλήψεων της προσομοίωσης, η εκτίμηση του \(\theta\) είναι \(\theta = \bar{Z}_M\), δηλαδή το ποσοστό των επαναλήψεων όπου \(Z=1\).

Η παρακάτω συνάρτηση υλοποιεί μια επανάληψη της προσομοίωσης της \(Z\):

sim_insurance_default = function(n, p, m, D)
{
  k=rbinom(1, n, p)
  S=sum(rexp(k, 1/m))
  Z=as.numeric(S>D)
  return(Z)
}

Για να εκτιμήσουμε την τιμή του \(\theta\) με τα δεδομένα του προβλήματος, έχουμε

n=500
p=0.05
m=1200
D=35000
N=1000
z=rep(0,N)
for (i in 1:N)
{
  z[i]=sim_insurance_default(n,p,m,D)
}
thetaest=mean(z)
print(thetaest)
## [1] 0.282

Επιπλέον μπορούμε να ερευνήσουμε τη συμπεριφορά του \(\theta\) όταν μεταβάλλονται μια ή περισσότερες παράμετροι του προβλήματος. Για παράδειγμα έστω ότι θέλουμε να δούμε πώς μεταβάλλεται το \(\theta\) ως προς το μέγεθος του αποθεματικού \(D\). Μπορούμε να ενσωματώσουμε την παραπάνω διαδικασία σε μια νέα συνάρτηση, την οποία μετά θα καλέσουμε για τιμές του \(D\) στο διάστημα \([25000,50000]\):

default_prob_estimate=function(n, p, m, D, N)
{
  z=rep(0,N)
  for (i in 1:N)
  {
    z[i]=sim_insurance_default(n,p,m,D)
  }
  thetaest=mean(z)
  return(thetaest)
}
n=500
p=0.05
m=1200
N=1000
Dval=seq(25000, 50000, 500)
nD=length(Dval)
thetaval=rep(0, nD)
for (i in 1:nD)
{
  thetaval[i]=default_prob_estimate(n,p,m,Dval[i],N)
}
plot(thetaval~Dval, type="l", xlab="D (reserve)", ylab="Probability of Default")

Η καμπύλη των εκτιμήσεων δεν είναι πολύ ομαλή. Αυτό οφείλεται στο ότι για κάθε τιμή του \(D\) η προσομοίωση δίνει μια εκτίμηση της τιμής του \(\theta\). Αν αυξήσουμε τον αριθμό των επαναλήψεων σε 10000 για κάθε τιμή του \(D\), η καμπύλη που προκύπτει είναι πολύ ομαλότερη:

N=10000
thetaval=rep(0, nD)
for (i in 1:nD)
{
  thetaval[i]=default_prob_estimate(n,p,m,Dval[i],N)
}
plot(thetaval~Dval, type="l", xlab="D (reserve)", ylab="Probability of Default")
abline(0.05,0, col="red")

Παρατηρούμε ότι η πιθανότητα χρεωκοπίας μειώνεται όσο αυξάνει το αποθεματικό. Για παράδειγμα, για να μειωθεί η πιθανότητα χρεωκοπίας στο 5%, απαιτείται αποθεματικό περίπου 45000 ευρώ. Αυτή η τιμή στην ανάλυση κινδύνου στα οικονομικά ονομάζεται Value at Risk (VaR).

Άσκηση 3.

Η ηλικία των ατόμων ενός πληθυσμού ακολουθεί κατανομή \(Gamma(5, 0.1)\). Η αρτηριακή πίεση ενός ατόμου ηλικίας \(x\) ακολουθεί κανονική κατανομή με μέση τιμή \(9+x/12\) και τυπική απόκλιση \(0.04 x\).

(α) Να περιγραφεί γεννήτρια τυχαίων αριθμών για την αρτηριακή πίεση ενός τυχαίου ατόμου αυτού του πληθυσμού.

(b) Αν η αρτηριακή πίεση άνω του 15 θεωρείται υπερβολική, να εκτιμηθεί το μέσο επίπεδο υπέρβασης από τη φυσιολογική τιμή : (i) για όλο τον πληθυσμό και (ii) για άτομα άνω των 65 ετών.

Απάντηση

Θα μοντελοποιήσουμε το πρόβλημα παραμετρικά. Έστω \(X\) η ηλικία ενός τυχαίου ατόμου και \(Y\) η αρτηριακή του πίεση. Δίνεται η περιθώρια κατανομή της \(X\): \(X\sim Gamma(\alpha, \lambda)\) και η δεσμευμένη κατανομή της \(Y\): \(Y|X=x \sim N(\beta_0+ \beta_1 x, (d x)^2)\).

(α) Δημιουργούμε μια γεννήτρια τυχαίων αριθμών για το ζεύγος \((X, Y)\). Η συνάρτηση προσομοιώνει πρώτα την \(X\) από την περιθώρια κατανομή και με βάση την τιμή \(X=x\), προσομοιώνει μια τιμή για την \(Y\) από την αντίστοιχη δεσμευμένη:

age_bp_generator=function(a, l, b0, b1, d)
{
  x=rgamma(1, a, l)
  y=rnorm(1, b0+b1*x, d*x)
  return(c(x,y))
}

Τώρα δημιουργούμε 1000 τυχαίες παρατηρήσεις ηλικίας-αρτηριακής πίεσης με τις δοθείσες τιμές των παραμέτρων:

a=5
l=0.1
b0=9
b1=1/12
d=0.04
N=1000
xy=matrix(rep(0, 2*N),N, 2)
for (i in 1:N)
{
  xy[i,]=age_bp_generator(a, l, b0, b1, d)
}

Από το ψευδοτυχαίο δείγμα μπορούμε να πάρουμε αρκετές πληροφορίες. Για παράδειγμα, μπορούμε να εκτιμήσουμε τη μέση αρτηριακή πίεση ενός τυχαίου ατόμου, ή την πιθανότητα η αρτηριακή πίεση να υπερβαίνει την τιμή 15:

mean_bp=mean(xy[,2])
mean_bp
## [1] 13.04474
prob_high=mean(xy[,2]>15)
print(prob_high)
## [1] 0.202

Επίσης μπορούμε να δούμε το scatter plot του δείγματος :

plot(xy[,2]~xy[,1], cex=0.3, xlab="x", ylab="y")

Στο διάγραμμα φαίνεται η ετεροσκεδαστικότητα του προβλήματος, καθώς η διασπορά της \(Y\) αυξάνει με την τιμή της \(X\). Αν το δείγμα αυτό ήταν πραγματικό και θέλαμε να εκτιμήσουμε τις τιμές των \(\beta_0, \beta_1\), δεν θα μπορούσαμε να χρησιμοποιήσουμε το απλό γραμμικό μοντέλο, επειδή θα παραβιαζόταν η υπόθεση της ομοσκεδαστικότητας των καταλοίπων.

(β i) Για την εκτίμηση της μέσης υπέρβασης, μπορούμε να χρησιμοποιήσουμε το ψευδοτυχαίο δείγμα των 1000 παρατηρήσεων του προηγούμενου ερωτήματος, να ορίσουμε μια νέα μεταβλητή \(Z=\max(Y-15,0)\) και να υπολογίσουμε τη μέση τιμή της από το δείγμα:

Z=rep(1000,0)
for (i in 1:1000)
  Z[i]=max(xy[i,2]-15,0)
mean_excess = mean(Z)
print(mean_excess)
## [1] 0.4570077

(β ii) Για την εκτίμηση της μέσης υπέρβασης σε άτομα ηλικίας άνω των 65 θα μπορούσαμε από το προηγούμενο δείγμα να επιλέξουμε μόνο τις παρατηρήσεις που έχουν τιμή \(X \geq 65\) και να υπολογίσουμε τη μέση τιμή της \(Z\) μόνο για αυτές:

mean_excess_above_65 = mean(Z[xy[,1]>=65])
print(mean_excess_above_65)
## [1] 1.701421

Ένα μειονέκτημα που έχει αυτή η προσέγγιση είναι ότι δεν μπορούμε να ελέγξουμε εξ αρχής το μέγεθος δείγματος με βάση το οποίο υπολογίζεται η μέση τιμή, επειδή δεν είναι γνωστό εξ αρχής πόσα άτομα στο αρχικό δείγμα θα έχουν ηλικία άνω των 65. Έστω ότι θέλουμε να εκτιμήσουμε την επιθυμητή ποσότητα με βάση ένα δείγμα \(N\) ατόμων ηλικίας άνω των 65. Αυτό μπορούμε να το κάνουμε με δύο εναλλακτικούς τρόπους.

Κατά πρώτον θα μπορούσαμε να τροποποιήσουμε τον αλγόριθμο δημιουργίας \(N\) παρατηρήσεων που χρησιμοποιήσαμε στο ερώτημα (α), με τρόπο ώστε να δημιουργεί παρατηρήσεις από το ζεύγος \((X,Y)\) αλλά να διατηρεί μόνο αυτές που έχουν \(X\geq 65\) και τις άλλες να τις απορρίπτει, μέχρι να συμπληρωθεί ο απαιτούμενος αριθμός:

N=1000
xy_65=matrix(rep(0, 2*N),N, 2)
n=0
while (n<N)
{
  w=age_bp_generator(a, l, b0, b1, d)
  if (w[1] >=65)
  {
    n=n+1
    xy_65[n, ]=w
  }
}
Z_65=rep(0,N)
for (i in 1:N)
{
  Z_65[i] = max(xy_65[i,2]-15,0)
}
mean_excess_above_65_ii=mean(Z_65)
print(mean_excess_above_65_ii)
## [1] 1.923157

Εναλλακτικά, θα μπορούσαμε να τροποποιήσουμε τη συνάρτηση που δημιουργεί τις παρατηρήσεις \((X,Y)\) έτσι ώστε να δημιουργεί παρατηρήσεις με τιμές \(X\geq A\), όπου \(A\) ορίζεται παραμετρικά:

age_bp_generator_above_Α=function(a, l, b0, b1, d, Α)
{
  x=rgamma(1, a, l)
  while (x<A) x=rgamma(1, a, l)
  y=rnorm(1, b0+b1*x, d*x)
  return(c(x,y))
}

Στην παραπάνω συνάρτηση μπορούμε να αποδείξουμε ότι η τιμή της \(X\) που προκύπτει ακολουθεί τη δεσμευμένη κατανομή \(X|X\geq A\).

Επομένως μπορούμε να δημιουργήσουμε 1000 παρατηρήσεις από αυτή τη γεννήτρια με \(A=65\) και να υπολογίσουμε τη μέση τιμή της υπέρβασης:

N=1000
A=65
xy_65_iii=matrix(rep(0, 2*N),N, 2)
for (i in 1:N)
{
  xy_65_iii[i, ] = age_bp_generator_above_Α(a, l, b0, b1, d, Α)
}
Z_65_iii=rep(0,N)
for (i in 1:N)
{
  Z_65_iii[i] = max(xy_65_iii[i,2]-15,0)
}
mean_excess_above_65_iii=mean(Z_65_iii)
print(mean_excess_above_65_iii)
## [1] 1.846244

Άσκηση 4.

Έστω ένα μοντέλο παλινδρόμησης με εξαρτημένη μεταβλητή \(Y\) και ανεξάρτητη \(X\), όπου η δεσμευμένη κατανομή της \(Y\) δεδομένου ότι \(X=x\) είναι εκθετική με μέση τιμή \(0.4 x + 2\). Υποθέτουμε ότι αυτή η σχέση δεν είναι γνωστή, και γίνεται μια στατιστική μελέτη για τη διερεύνηση της συσχέτισης μεταξύ των δύο μεταβλητών. Για τη μελέτη συγκεντρώνεται ένα δείγμα μεγέθους 100 στο οποίο υπάρχουν οι τιμές της \(X=1,2,\ldots, 10\), από 10 φορές η καθεμιά. Ο αναλυτής δε γνωρίζει την πραγματική κατανομή των καταλοίπων. Υποθέτει ότι αυτή είναι κανονική με σταθερή διασπορά, έτσι ώστε να μπορεί να εφαρμόσει το μοντέλο γραμμικής παλινδρόμησης. Μετά τον υπολογισμό των εκτιμήσεων ελαχίστων τετραγώνων γίνεται υπολογισμός των καταλοίπων jacknife και έλεγχος κανονικότητας αυτών μέσω του ελέγχου .

(α) Χρησιμοποιώντας προσομοίωση για τις τιμές της \(Y\) στο δοσμένο δείγμα, να εκτιμηθεί η μέση απόκλιση των εκτιμήσεων ελαχίστων τετραγώνων \(\hat{\beta}_0, \hat{\beta}_1\) από τις αληθινές τιμές 2 και 0.4, αντίστοιχα.

(β) Να εκτιμηθεί η πιθανόητα απόρριψης του ελέγχου κανονικότητας.

Απάντηση

(α) Ορίζουμε το διάνυσμα τιμών της \(X\), και μια συνάρτηση η οποία εκτελεί μια επανάληψη της προσομοίωησης, συγκεκριμένα: (1) δημιουργεί ένα διάνυσμα Υ από την κατανομή που δίνεται, (2) εκτελεί το γραμμικό μοντέλο \(Y=\beta_0 + \beta_1 X + \epsilon\), με την υπόθεση ότι \(\epsilon \sim N(0,\sigma^2)\), (3) υπολογίζει τα κατάλοιπα jacknife και εκτελεί έναν έλεγχο Shapiro Wilk για κανονικότητα και (4) επιστρέφει τις τιμές των εκτιμήσεων ελαχίστων τετραγώνων \(\hat{\beta}_0, \hat{\beta}_1\) και την τιμή του pvalue του Shapiro Wilk test.

x=rep(1:10, 10)
simulate_lm=function() {
  y=rep(0,100)
  for (i in 1:100)
  {
    y[i]=rexp(1, rate=1/(0.4*x[i]+2))
  }
  lm.fit=lm(y~x)
  beta_hat=lm.fit$coefficients
  lm.jack=rstudent(lm.fit)
  sw=shapiro.test(lm.jack)
  p=sw$p.value
  return(c(beta_hat, p))
}

Τώρα εκτελούμε 1000 επαναλήψεις της προσομοίωσης και υπολογίζουμε τις μέσες τιμές των εκτιμήσεων και του pvalue. Επίσης εκτιμάμε την πιθανότητα απόρριψης στο επίπεδο 5%, υπολογίζοντας το ποσοστό των επαναλήψεων στις οποίες προκύπτει \(pvalue< 0.05\):

N=1000
results=matrix(rep(0, 3*N), N, 3)
for (i in 1:N)
{
  results[i,] = simulate_lm()
}
means=apply(results, 2, mean)
b0=means[1]
b1=means[2]
pval=means[3]
prob_reject=mean(results[,3]<0.05)
cat(sprintf("b0 = %.3f, b1 = %.3f, p-value = %.3f, prob_reject = %.3f\n", b0, b1, pval, prob_reject))
## b0 = 2.041, b1 = 0.394, p-value = 0.000, prob_reject = 1.000

Άσκηση 5.

Έστω το ολοκλήρωμα \[ A = \int_0^{\infty} \frac{\sqrt{x+4}}{x^3+1} dx. \]

(α) Να περιγράψετε μια μέθοδο εκτίμησης της τιμής του \(A\) με τη μέθοδο προσομοίωσης .

(β) Να υπολογίσετε τον απαιτούμενο αριθμό παρατηρήσεων έτσι ώστε η τιμή του \(A\) να εκτιμηθεί με διάστημα εμπιστοσύνης 95% και ακρίβεια \(\pm 0.01\).

(γ) Να βρεθεί το παραπάνω Δ.Ε.

Απάντηση

(α) Έστω \(X\) τυχαία μεταβλητή που ακολουθεί εκθετική κατανομή με παράμετρο 1. Η συνάρτηση πυκνότητας πιθανότητας της \(X\) είναι \(f(x) = e^{-x}, x\geq 0\). Το ολοκλήρωμα \(A\) μπορεί να γραφτεί \[ A = \int_0^{\infty} \frac{\sqrt{x+4}}{x^3+1} e^x e{-x} dx = \int_0^{\infty} \frac{\sqrt{x+4}}{x^3+1} e^x f(x) dx = E(Y), \] όπου η τυχαία μεταβλητή \(Y = \frac{\sqrt{X+4}}{X^3+1} e^X\). Επομένως μια μέθοδος εκτίμησης της τιμής του ολοκληρώματος είναι μέσω προσομοίωησης Monte Carlo, θεωρώντας το \(A\) ως μέση τιμή της \(Y\). Δημιουργούμε \(N\) παρατηρήσεις της \(X\) από την εκθετική κατανομή με παράμετρο 1, υπολογίζουμε τις αντίστοιχες τιμές της \(Υ\) και εκτιμούμε το \(A\) από τη μέση τιμή.

N=1000
x=rexp(N, 1)
y=sqrt(x+4)*exp(x)/(x^3+1)
A_estimate=mean(y)
cat(sprintf("Estimate of A = %.4f", A_estimate))
## Estimate of A = 2.6730

(β) Θέλουμε να εκτιμήσουμε τη μέση τιμή της τυχαίας μεταβλητής \(Y\) με ένα διάστημα εμπιστοσύνης 95% και ακρίβεια 0.01, δηλαδή το μισό μήκος του διαστήματος εμπιστοσύνης να είναι ίσο με \(h=0.01\) έτσι ώστε με πιθανότητα 95% το διάστημα \(\bar{Y} \plusminus h\) να περιέχει την πραγματική μέση τιμή.

Έστω \(\sigma\) η τυπική απόκλιση της \(Y\). Ο τύπος για το μισό μήκος είναι \[ h=\frac{\sigma z_{1-\alpha/2}}{\sqrt{N}}, \] όπου \(N\) το μέγεθος δείγματος. Λύνοντας ως προς \(N\) παίρνουμε ότι το απαιτούμενο μέγεθος δείγματος είναι \[ N_0 = \frac{\sigma^2 z^2_{1-\alpha/2}}{h^2}, \] Στον παραπάνω τύπο δεν γνωρίζουμε την διασπορά \(\sigma^2\). Αν χρησιμοποιήσουμε ως εκτίμηση της διασποράς τη δειγματική διασπορά του \(Y\) από το προηγούμενο ερώτημα μπορούμε να πάρουμε μια προσέγγιση για το απαιτούμενο μέγεθος δείγματος:

a=0.05
h=0.01
v=var(y)
z=qnorm(1-a/2, 0, 1)
N0=round(v*z^2/h^2)
print(N0)
## [1] 43300

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

x1=rexp(N0, 1)
y1=sqrt(x1+4)*exp(x1)/(x1^3+1)
y1_mean=mean(y1)
s1=sd(y1)
h=s1*qt(1-a/2, N0-1)/sqrt(N0)
L=y1_mean-h
U=y1_mean+h
cat(sprintf("Εκτίμηση A = %.4f, Διάστημα Εμπιστοσύνης = (%.4f, %.4f)", y1_mean, L, U))
## Εκτίμηση A = 2.6560, Διάστημα Εμπιστοσύνης = (2.6487, 2.6633)

Άσκηση 6.

Μια βιολογική μέτρηση \(X\) είναι τυχαία μεταβλητή που περιγράφεται από το μοντέλο \(X=(1-U^2)e^U\), όπου η \(U\) ακολουθεί ομοιόμορφη κατανομή στο διάστημα \((0,1)\). Έστω \(\theta=E(X)\) η μέση τιμή της \(X\).

Να εκτιμηθεί το \(\theta\) με ένα διάστημα εμπιστοσύνης 95% βασισμένο σε 1000 προσομοιωμένες παρατηρήσεις της \(X\).

Απάντηση

Η σημειακή εκτίμηση και το διάστημα εμπιστοσύνης υπολογίζονται ως εξής:

N=1000
a=0.05
u=runif(N, 0,1)
x=(1-u^2)*exp(u)
xmean=mean(x)
xs=sd(x)
h=xs*qt(1-a/2, N-1)/sqrt(N)
L=xmean-h
U=xmean+h
cat(sprintf("Εκτίμηση του θ = %.4f, Διάστημα Εμπιστοσύνης = (%.4f, %.4f)", xmean, L, U))
## Εκτίμηση του θ = 1.0139, Διάστημα Εμπιστοσύνης = (0.9952, 1.0326)

Άσκηση 7.

Έστω το ζεύγος τυχαίων μεταβλητών \((X,Y)\) με από κοινού συνάρτηση κατανομής πιθανότητας \(f(x,y) = k x, 0\leq x \leq y \leq 1\).

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

(β) Να εκτιμηθούν η μέση τιμή και η διασπορά των \(X,Y\) όπως επίσης και ο συντελεστής συσχέτισής τους.

Απάντηση

Η συνάρτηση πυκνότητας πιθανότητας είναι \[ f(x,y) = \left \{ \begin{array}{ll} k x, & 0\leq x \leq 1, 0 \leq y \leq 1, x \leq y \\ 0, & \mbox{διαφορετικά} \end{array} \right . \] Για τη δημιουργία γεννήτριας του ζεύγους \((X,Y)\) μέσω της μεθόδου MCMC ακολουθούμε τα εξής βήματα :

  1. Ορίζουμε μια Μαρκοβιανή διαδικασία με κατάσταση το διδιάστατο διάνυσμα \(s=(x,y)\) και χώρο καταστάσεων το σύνολο \(\{(x,y): 0 \leq x \leq 1, 0 \leq y \leq 1\}\).

Για ευκολία μπορούμε να πάρουμε τις καταστάσεις να είναι ανεξάρτητες από τη μια περίοδο στην άλλη και σε κάθε κατάσταση οι \(X,Y\) να είναι ανεξάρτητες με ομοιόμορφη κατανομή μεταξύ 0 και 1. Με αυτή την υπόθεση οι πιθανότητες μετάβασης ενός βήματος από την κατάσταση \((x,y)\) στην κατάσταση \((x', y')\) είναι \(q_{(x,y), (x',y')} = 1, 0\leq x' \leq 1, 0 \leq y' \leq 1\) ανεξάρτητα από την αρχική κατάσταση \((x,y)\).

  1. Τώρα ορίζουμε την τελική Μαρκοβιανή διαδικασία με τον ίδιο χώρο καταστάσεων, η οποία εξελίσσεται ως εξής: Ξεκινάμε από μια αυθαίρετη κατάσταση \((x,y)\) με \(0 \leq x \leq y \leq 1\). Δημιουργούμε ένα ζεύγος \((x', y')\) από ανεξάρτητες ομοιόμορφες παρατηρήσεις μεταξύ 0 και 1. Αυτή η παρατήρηση γίνεται δεκτή με πιθανότητα \[ \alpha =\min \left \{1, \frac{f(x', y') q_{(x,y), (x',y')}}{f(x, y) q_{(x',y'), (x,y)}} \right \} \] Αυτός ο τύπος απλοποιείται σημαντικά με τις κατανομές που έχουμε υποθέσει: Πρώτον, \(q_{(x,y), (x',y')} = q_{(x',y'), (x,y)} =1\). Επίσης, επειδή σε κάθε βήμα η τρέχουσα κατάσταση \((x,y)\) ικανοποιεί την \(0 \leq x \leq y \leq 1\), έχουμε \(f(x,y) = kx\). Τέλος, για τη νέα κατάσταση \((x', y')\) ισχύει \(0\leq x' \leq 1\) και \(0 \leq y' \leq 1\), όμως μπορεί να μην ισχύει \(x' \leq y'\). Επομένως \(f(x', y') = kx'\) αν \(x' \leq y'\) και 0 διαφορετικά. Μετά από τις παραπάνω απλοποιήσεις η πιθανότητα αποδοχής γίνεται \[ \alpha = \min \left \{1, \frac{x'}{x} \right \}, \] αν \(x' \leq y'\) και 0 διαφορετικά.

Επομένως η νέα κατάσταση γίνεται δεκτή αν \(x' \leq y'\) και \(x' \geq x\), (επειδή τότε \(\alpha=1\)), γίνεται δεκτή με πιθανότητα \(x'/x\) αν \(x' \geq y'\) και \(x' \geq x\), και απορρίπτεται αν \(x' > y'\). Στην τελευταία περίπτωση η κατάσταση παραμένει στο \((x,y)\).

  1. Ξεκινώντας από μια αυθαίρετη αρχική κατάσταση \((x,y)\) με \(0 \leq x \leq y \leq 1\), επαναλαμβάνουμε την παραπάνω διαδικασία για ένα μεγάλο αριθμό βημάτων (π.χ. 40 βήματα) και καταγράφουμε την κατάσταση \((x,y)\) που έχει προκύψει στο τελευταίο βήμα. Αυτή η κατάσταση με καλή προσέγγιση ακολουθεί την επιθυμητή κατανομή.

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

mcmc_generator = function(x0, y0, n)
{
  if ( (x0>=0)&(x0<=y0)&(y0<=1)) initial_accept=TRUE else initial_accept=FALSE
  if (initial_accept==FALSE) 
  {
    print("Initial state (x,y) not acceptable")
    return(NULL)
  } else 
  {
    t=0
    x=x0
    y=y0
    for (t in 1:n)
    {
      x1=runif(1,0,1)
      y1=runif(1,0,1)
      if (x1 > y1) 
      {
        accept_new=FALSE
      }
      else if (x1>x)
      {
        accept_new=TRUE
      }
      else
      {
        u=runif(1,0,1)
        accept_new=(u<x1/x)
      }
      if (accept_new==TRUE) 
      {
        x=x1
        y=y1
      }
    }
  }
  return(c(x,y))
}

(β)

Δημιουργούμε 1000 παρατηρήσεις από τη γεννήτρια που σχεδιάσαμε στο (α) και σχηματίζουμε ένα scatterplot

N=1000
xy=matrix(rep(0,2*N),N,2)
for (i in 1:N)
{
  xy[i,]=mcmc_generator(0,0,40)
}
plot(xy[,2]~xy[,1], xlab="x", ylab="y", pch="x", cex=0.4)

Από το δείγμα μπορούμε εύκολα να εκτιμήσουμε μέσες τιμές, διασπορές και συσχέτιση

mx=mean(xy[,1])
my=mean(xy[,2])
vx=var(xy[,1])
vy=var(xy[,2])
rho=cor(xy[,1], xy[,2])
cat(sprintf("Mean(X) = %.3f, Mean(Y) = %.3f, Var(X) = %.3f, Var(Y) = %.3f, Correlation = %.3f", mx, my, vx, vy, rho))
## Mean(X) = 0.497, Mean(Y) = 0.741, Var(X) = 0.053, Var(Y) = 0.040, Correlation = 0.599