CODEWITHSUNDEEP
r
user3609848
user3609848

Reputation: 141

Adding two random variables via convolution in R

I would like to compute the convolution of two probability distributions in R and I need some help. For the sake of simplicity, let's say I have a variable x that is normally distributed with mean = 1.0 and stdev = 0.5, and y that is log-normally distributed with mean = 1.5 and stdev = 0.75. I want to determine z = x + y. I understand that the distribution of z is not known a priori.

As an aside the real world example I am working with requires addition to two random variables that are distributed according to a number of different distributions.

Does anyone know how to add two random variables by convoluting the probability density functions of x and y?

I have tried generating n normally distributed random values (with above parameters) and adding them to n log-normally distributed random values. However, I wish to know if I can use the convolution method instead. Any help would be greatly appreciated.

EDIT

Thank you for these answers. I define a pdf, and try to do the convolution integral, but R complains on the integration step. My pdfs are Log Pearson 3 and are as follows 

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 complains at the last step since the f.t function is bounded and my integration limits are probably not correct. Any ideas on how to solve this?

Upvotes: 13

Views: 13195

Answers (2)

holdmygruyere
holdmygruyere

Reputation: 31

I was having trouble getting integrate() to work for different density parameters, so I came up with an alternative to @jlhoward's using Riemann approximation:

set.seed(1)

#densities to be convolved. could also put these in the function below
d1   <- function(x) dnorm(x,1,0.5) #
d2   <- function(y) dlnorm(y,1.5, 0.75) 

#Riemann approximation of convolution 
conv <- function(t, a, b, d) { #a to b needs to cover the range of densities above. d needs to be small for accurate approx.
  z <- NA
  x <- seq(a, b, d)
  for (i in 1:length(t)){
    print(i)
    z[i] <- sum(d1(x)*d2(t[i]-x)*d)
  }
  return(z)
}

#check against sampled convolution
X <- rnorm(1000, 1, 0.5)
Y <- rlnorm(1000, 1.5, 0.75)
Z <- X + Y

t <- seq(0, 50, 0.05) #range to evaluate t, smaller increment -> smoother curve

hist(Z, breaks = 50, freq = F, xlim = c(0,30))
lines(t, conv(t, -100, 100, 0.1), type = "s", col = "red")

plot

Upvotes: 2

jlhoward
jlhoward

Reputation: 59365

Here is one way.

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")

Same thing using package 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")

Upvotes: 14

Related Questions