2010-08-25 9 views

risposta

8

Questo è abbastanza semplice utilizzando il " distr" pacchetto:

library(distr) 

A <- Exp(rate=3) 
B <- Gammad(shape=2, scale=3) 

conv <- 0.5*(A+B) 

plot(conv) 
plot(conv, to.draw.arg=1) 

Edit da JD lungo

grafico risultante è simile al seguente: alt text

+0

È abbastanza carino. Non avevo familiarità con il pacchetto distr. –

1

Io non sono un programmatore di R, ma potrebbe essere utile sapere che per le variabili casuali indipendenti con i PDF f (x) ed f 2 (x), il PDF della somma delle due variabili è data dalla convoluzione f * f (x) dei due PDF ingresso.

2

Ecco un tentativo di convoluzione (a cui si riferisce @Jim Lewis) in R. Si noti che probabilmente ci sono modi molto più efficienti per farlo.

lower <- 0 
upper <- 20 
t <- seq(lower,upper,0.01) 
fA <- dexp(t, rate = 0.4) 
fB <- dgamma(t,shape = 8, rate = 2) 
## C has the same distribution as (A + B)/2 
dC <- function(x, lower, upper, exp.rate, gamma.rate, gamma.shape){ 
    integrand <- function(Y, X, exp.rate, gamma.rate, gamma.shape){ 
    dexp(Y, rate = exp.rate)*dgamma(2*X-Y, rate = gamma.rate, shape = gamma.shape)*2 
    } 
    out <- NULL 
    for(ix in seq_along(x)){ 
    out[ix] <- 
     integrate(integrand, lower = lower, upper = upper, 
       X = x[ix], exp.rate = exp.rate, 
       gamma.rate = gamma.rate, gamma.shape = gamma.shape)$value 
    } 
    return(out) 
} 
fC <- dC(t, lower=lower, upper=upper, exp.rate=0.4, gamma.rate=2, gamma.shape=8) 
## plot the resulting distribution 
plot(t,fA, 
    ylim = range(fA,fB,na.rm=TRUE,finite = TRUE), 
    xlab = 'x',ylab = 'f(x)',type = 'l') 
lines(t,fB,lty = 2) 
lines(t,fC,lty = 3) 
legend('topright', c('A ~ exp(0.4)','B ~ gamma(8,2)', 'C ~ (A+B)/2'),lty = 1:3) 
+0

Non preoccupatevi di generazione dei dati dalle funzioni, basta tracciare le funzioni. Dovresti usare curve() invece di lines() con add = TRUE. – John

5

Se stai cercando un grafico veloce, di solito faccio l'approccio di simulazione veloce e sporco. Faccio un po disegna, Slam una densità gaussiana sulla pareggi e trama che cattivo ragazzo:

numDraws <- 1e6 
gammaDraws <- rgamma(numDraws, 2) 
expDraws <- rexp(numDraws) 
combined <- .5 * (gammaDraws + expDraws) 
plot(density(combined)) 

output dovrebbe apparire un po 'come questo:

alt text

Problemi correlati