Obiettivo

imparare a generare da qualunque distribuzione di probabilità continua e discreta univariata.

Metodo dell’Inversione della fdr (ITM)

Vediamo il caso Y~f , U = F(Y) con f pdf della v.c. Y ed F la sua fdr. Vale che U=F(Y)~U(0,1) motivo per cui campiono n valori da U(0,1) e calcolo l’inversa F^-1 che mi riporta a y = F^-1(u)

# Y~ Exp(theta) F(y) = 1 - exp(-y/theta) = u
# y = -ln(1-u)*theta

theta <- 3
n <- 1000
u <- runif(n)
y <- -log(1-u)*theta

# confronto con momenti teorici di un'esp. theta
df <- data.frame(campionari=c(mean(y), var(y)), teorici=c(theta, theta^2))
df
##   campionari teorici
## 1   2.921891       3
## 2   8.758949       9
# Y~ Weibull(lambda, k) F(y) = 1 - exp(-(y/lambda)^k) = u
# y = (-1*ln(1-u))^(1/k)*lambda
lambda <- 3
k <- 2
y.weibull <- lambda*(-1*log(1-u))^(1/k)
hist(y.weibull, freq=F)
curve(dweibull(x, k, lambda), add=T)

# Y~ Bin(k,p) ; supporto = {0,1,..,k}
k = 10
p = .8
supporto = 0:k
Pi = pbinom(q=0:k, size=k, prob=p)
# campiono u da U(0,1) e fisso y=min(k| Pk>=u)
u = runif(n)
y = c()
for (i in 1:n) {
y[i] = supporto[min(which(Pi>=u[i]))]
}
hist(y, breaks=0:k, freq=F)

Metodo della Trasformazione di v.c. più semplici

L’idea è quella di sfruttare relazioni note nella teoria delle distribuzioni per ottenere valori da distribuzioni complesse a partire da una funzione di valori estratti da distribuzioni più semplici.

Ad esempio, la somma di N esponenziali di parametro THETA è una v.c. Y di distribuzione nota: Gamma(N, THETA)

# Genera 1000 valori per Y Gamma(alpha, beta) 
alpha <- 4
beta <- 5
n <- 1000
w <- matrix(rexp(n=alpha*n, rate=beta), ncol=alpha, nrow=n)
y <- apply(w,MARGIN=1,FUN=sum) 

hist(y, freq=F, main='', breaks=50, col='lightblue')
curve(dgamma(x,shape=alpha,rate=beta), add=T, lwd=3)

Oppure, il rapporto tra due somme di ampiezza (rispettivamente) A,(A+B) di esponenziali negativi di parametro 1 è una BETA(A,B)

# Genera 1000 valori per Y Beta(a,b)
a <- 2
b <- 4
n <- 1000
w1 <- matrix(rexp(n=a*n), ncol=a, nrow=n)
w2 <- cbind(w1,matrix(rexp(n=b*n), ncol=b, nrow=n))
num <- apply(w1, MARGIN=1, FUN=sum)
den <- apply(w2, MARGIN=1, FUN=sum)
y <- num/den
hist(y, breaks=50, freq=F)
curve(dbeta(x, a, b), add=T)

Potremmo anche essere interessati a generare delle MISTURE di probabilità. Una v.c. X si dice MISTURA delle variabili casuali A,B,C, con (rispettivamente) funzioni di ripartizione F1(),F2(),F3() se e solo se vale F(x) = p1 F1(x) + p2 F2(x) + p3 F3(x) , dove sum(p[i]) = 1.

# Genero 100000 valori da Y mistura di 3 Normali
# A ~ N(0,1), B ~ N(7,1) , C ~ N(2,0.5) 
# f.Y(y) = p1*f.A(y) + p2*f.B(y) + p3*f.C(y)

mu <- c(0,7,2)
sd <- c(1,1,1/2)
p <- c(1/3, 1/4, 5/12)

rnmix = function(n,mu,sd,p){ 
index = sample(1:3,size=n,replace=TRUE, prob=p)
rnorm(n,mean=mu[index],sd=sd[index]) 
}

y <- rnmix(n=100000, mu=mu, sd=sd, p=p) 
hist(y, breaks=100, freq=F)
curve(p[1]*dnorm(x,mean=mu[1], sd=sd[1]), add=T, lwd=2, col='blue')
curve(p[2]*dnorm(x,mean=mu[2], sd=sd[2]), add=T, lwd=2, col='blue')
curve(p[3]*dnorm(x,mean=mu[3], sd=sd[3]), add=T, lwd=2, col='blue')

Algoritmo Accettazione-Rifiuto (AR)

Voglio un campione da Y con pdf f(y) ‘diffile’, quindi scelgo una distribuzione candidata X con pdf g(x) da cui è più facile campionare e t.c. il rapporto f(x)/g(x) <= K per ogni x.

# Voglio generare R=1000 valori da Y~ Beta(2,3)
# Scelgo X~U(0,1) <-> g(x) = 1 
a <- 2; b<-3; R=1000
f <- function(y,a,b) {dbeta(y,a,b)}

# step 0. Calcolo K = max wrt x di f/g
K <- optimize(f,interval=c(0,1),a=a, b=b, maximum=T)$objective

# step 1. Campiono in modo indipendente u da U(0,1) e x=candidato da g(X).
u <- runif(R)
x <- runif(R)

# step 2. Se u<=f(x)/(K*g(x)) -> accetto il candidato e pongo y=x
#         Altrimenti rifiuto e torno al punto 1
test <- u<=f(x,a,b)/K
y <- x[test==T]

prop.accettazione <- length(y)/R
1/K 
## [1] 0.5625
hist(y, breaks=40, freq=F)
curve(f(x,a=2,b=3), from=0, to=1, add=T)

Nell’applicazione del metodo è cruciale scegliere una v.a. X ausiliaria che abbia supporto coincidente o contenuto a quello di Y; inoltre fissato R (num.simulazioni) il tasso di accettazione dei candidati x estratti da g() è circa pari a 1/K, con K=max(f/g).
Vediamo un altro esempio:

# genera 1000 determinazioni dalla v.c. Y con f(y)=3*y^2 con y in [0,1]
# scegli g(x) = 1 

f <- function(x) {3*x^2}
K <- optimize(f,interval=c(0,1),maximum=T)$objective
prop.acc. <- 1/K 

# circa 1/3 dei valori generati sarà accettato; genero 3000 oss.

x <- runif(3000)
u <- runif(3000)
test <- u<= f(x)/K
y <- x[test==T]
hist(y, breaks=50, freq=F)
curve(f, from=0, to=1, add=T)

Metodi ad-hoc per generare da v.c. Normale

1.Metodo basato sul TLC

nsim <- 5000
u <- matrix(data=runif(nsim*12), ncol=12, nrow=nsim)
z <- apply(u, 1, function(x) sum(x) - 6)
hist(z, freq=F)
curve(dnorm, from=-3, to=3, add=T)

2.Metodo delle COORDINATE POLARI o BOX-MUELLER

# Vogliamo generare osservazioni (z1,z2) iid da N(0,1)
# Applico trasf. (z1,z2) -> (r,theta) con
# r = (z1^2 + z2^2)^(1/2); R^2 ~ exp(1/2)
# theta = arctg(z2/z1); THETA ~ U(0,2*pi)
R = 10000
u1 <- runif(R)
u2 <- runif(R)

# campiono (r,theta)
r <- sqrt(-2*log(u1))
theta <- 2*pi*u2

# riapplico trasf. inversa per ottenere (z1,z2)
z1 <- r*cos(theta)
z2 <- r*sin(theta)

# controllo che z1 è N(0,1)
hist(z1,freq=F)
curve(dnorm, from=-4, to=4, add=T)

3. Rejection Polar Method

# voglio generare (z1,z2) iid N(0,1)
nsim = 5000
u1 <- runif(nsim); u2 <- runif(nsim)
# v1,v2 iid U(-1,1)
v1 <- 2*u1 - 1; v2 <- 2*u2 - 1
s <- v1^2 + v2^2
test <- s<1
s <- s[test==T]
# z1,z2
z1 <- v1[test==T]*sqrt(-2*log(s)/s)
z2 <- v2[test==T]*sqrt(-2*log(s)/s)

# test empirico z1 N(0,1)
hist(z1, freq=F)
curve(dnorm, from=-4, to=4, add=T)

c(mean(z1),var(z1))
## [1] -0.01443862  0.98933291

Metodo del rapporto di Uniformi

Voglio generare determinazioni da X con pdf f(x).
Considero funz. h(x) t.c. f = h/integrale(h).
Campiono (u,v) da due uniformi U(0,a), U(bminus,bplus),
con a = sup h, bminus = inf xh, bplus = sup xh
Se (u,v) : 0<=u<=sqrt(h(v/u)) -> x = v/u

# f Cauchy : f(x) = 1/pi * 1/(1+x^2)
# h(x) = 1/(1+x^2)
nsim = 5000
h <- function(x) {sqrt(1/(1+x^2))}
xh <- function(x) {x*sqrt(1/(1+x^2))}
a <- optimize(h, c(-5,5),maximum=T)$objective
bminus <- optimize(xh, c(-5,5))$objective
bplus <- optimize(xh, c(-5,5), maximum=T)$objective
# campiono (u,v)
u <- runif(nsim,0,a)
v <- runif(nsim,bminus,bplus)
# u <= sqrt(h(v/u)) è equivalente a u^2 <= 1 - v^2
test <- (u^2)<=(1-v^2)
x <- v/u
x1 <- x[test==T]
hist(x1,freq=F)
curve(pi*1/(1+x^2),add=T)

#f Normale(0,1)
#h prop. f 
nsim=5000
h <- function(x) {sqrt(exp(-x^2/2))}
xh <- function(x) {x*sqrt(exp(-x^2/2))}
a <- optimize(h, c(0,2), maximum=T)$objective
bminus <- optimize(xh, c(-5,0))$objective
bplus <- optimize(xh, c(0,2), maximum=T)$objective
# campiono (u,v) da U(0,a), U(bmeno,bpiù)
u <- runif(nsim, 0, a)
v <- runif(nsim, bminus, bplus)
# (u,v) in C(h) = {u<=sqrt(h(v/u))}
# u^2 <= exp(-0.5*(v/u)^2) 
# -4ln(u) >= (v/u)^2 
# v^2 <= -4*u^2*ln(u)
x <- v/u
test = (v^2) <= (-4*u^2*log(u))
x1 <- x[test==T]
hist(x1,freq=F)
curve(dnorm,add=T)