2014-05-09 10 views
9

Vorrei calcolare la convoluzione di due distribuzioni di probabilità in R e ho bisogno di aiuto. Per semplicità, diciamo che ho una variabile x che è normalmente distribuita con media = 1.0 e stdev = 0.5, e che è log-normalmente distribuita con media = 1.5 e stdev = 0.75. Voglio determinare z = x + y. Capisco che la distribuzione di z non è nota a priori.Aggiunta di due variabili casuali tramite convoluzione in R

Per inciso, l'esempio del mondo reale su cui sto lavorando richiede l'aggiunta di due variabili casuali distribuite in base a diverse distribuzioni.

Qualcuno sa come aggiungere due variabili casuali convolutando le funzioni di densità di probabilità di x e y?

Ho provato a generare n valori casuali distribuiti normalmente (con i parametri precedenti) e li ho aggiunti a n valori logici distribuiti normalmente. Tuttavia, vorrei sapere se posso usare invece il metodo di convoluzione. Qualsiasi aiuto sarebbe molto apprezzato.

EDIT

Grazie di queste risposte. Definisco un pdf e provo a fare l'integrale di convoluzione, ma R si lamenta del passaggio di integrazione. Le mie PDF sono Log Pearson 3 e sono i seguenti

dlp3 <- function(x, a, b, g) { 
p1 <- 1/(x*abs(b) * gamma(a)) 
p2 <- ((log(x)-g)/b)^(a-1) 
p3 <- exp(-1* (log(x)-g)/b) 
d <- p1 * p2 * p3 
return(d) 
} 

f.m <- function(x) dlp3(x,3.2594,-0.18218,0.53441) 
f.s <- function(x) dlp3(x,9.5645,-0.07676,1.184) 

f.t <- function(z) integrate(function(x,z) f.s(z-x)*f.m(x),-Inf,Inf,z)$value 
f.t <- Vectorize(f.t) 
integrate(f.t, lower = 0, upper = 3.6) 

R lamenta all'ultimo passo in quanto la funzione f.t è delimitata ed i miei limiti di integrazione non sono probabilmente corrette. Qualche idea su come risolvere questo?

+1

Vi suggerisco di controllare il [pacchetto Distr] (http://cran.r-project.org/web/packages/distr/index.html) o almeno prendere un rapido sguardo a Il [vignetta ] (http://cran.r-project.org/web/packages/distr/vignettes/newDistributions.pdf). Penso che sia esattamente quello che stai cercando. Anche se la strategia che hai usato per generare valori casuali è perfettamente valida. – MrFlick

risposta

12

Ecco un modo.

f.X <- function(x) dnorm(x,1,0.5)  # normal (mu=1.5, sigma=0.5) 
f.Y <- function(y) dlnorm(y,1.5, 0.75) # log-normal (mu=1.5, sigma=0.75) 
# convolution integral 
f.Z <- function(z) integrate(function(x,z) f.Y(z-x)*f.X(x),-Inf,Inf,z)$value 
f.Z <- Vectorize(f.Z)     # need to vectorize the resulting fn. 

set.seed(1)        # for reproducible example 
X <- rnorm(1000,1,0.5) 
Y <- rlnorm(1000,1.5,0.75) 
Z <- X + Y 
# compare the methods 
hist(Z,freq=F,breaks=50, xlim=c(0,30)) 
z <- seq(0,50,0.01) 
lines(z,f.Z(z),lty=2,col="red") 

stessa cosa usando il pacchetto distr.

library(distr) 
N <- Norm(mean=1, sd=0.5)   # N is signature for normal dist 
L <- Lnorm(meanlog=1.5,sdlog=0.75) # same for log-normal 
conv <- convpow(L+N,1)    # object of class AbscontDistribution 
f.Z <- d(conv)     # distribution function 

hist(Z,freq=F,breaks=50, xlim=c(0,30)) 
z <- seq(0,50,0.01) 
lines(z,f.Z(z),lty=2,col="red") 
+0

Per quanto riguarda le distribuzioni dlp3, sei sicuro che la tua equazione sia corretta? La tua funzione 'dlp3 (...)' diverge. Prova 'x <- seq (0,100, .1); plot (x, dlp3 (x, 3, -0.2,0.5)) '. In effetti, può essere mostrato seguire ~ x^4 * log (x)^2 per x grandi. – jlhoward