# Syntax template for convolving two probability distributions

I'd appreciate some help with the syntax of the following type of problem. Once I have a template of how to do it in Mathematica, I'll be able to expand on it and work with it, but I'm new and don't know how to get there.

I'd like to get an expression for the sum of random variables taken from the following two distributions [or equivalently the convolution of the PDFs of them]:

1) Log ( ( negative binomial with parameters r and p) + 0.5)

2) Log normal distribution with parameters mu and sigma

So my eventual distribution would be of the sum of a random variable from (1) and a random variable from (2).

Ultimately, I'd like to get the distribution of a sum of random variables from several different distributions, but once I know how to do it for two, I should be able to extrapolate from there.

Thanks very much for any help

• You want to look up TransformedDistribution. :) – Quantum_Oli Mar 10 '17 at 18:27
• You really want the distribution of the sum of a discrete and a continuous random variable? I don't think you're going to get a closed-form solution but you can certainly get random samples from such a distribution when you set numerical values for r, p, mu, and sigma. (The next best thing to getting an explicit distribution function is the ability to generate zillions of samples from that distribution.) – JimB Mar 10 '17 at 18:58

## 2 Answers

@mikado and @Quantum_Oli give you the basic way to determine the distribution of the sum of two (independent) random variables and that's probably the main thrust of your question: how to do it in general.

But with your specific example, the best you can hope for using such an approach is the ability to construct random samples from the distribution of the sum. And despite what I said in my comment above, for this example you can actually get a specific form for the distribution function.

Essentially one conditions on each value of the discrete random variable (the negative binomial variable which takes on positive probabilities for all non-negative integers) and then weight the resulting conditional distributions on the negative binomial probabilities. Such a function could look like the following:

cdf[z0_, r_, p_, \[Mu]_, \[Sigma]_] := If[z0 < Log[1/2], 0,
p^r Sum[Binomial[r + i - 1, r - 1] (1 - p)^i Erfc[-Log[z0 - Log[i + 1/2]]/2^(1/2)]/2,
{i, 0, Floor[Exp[z0] - 1/2]}]]


To use a specific example generate some data:

r = 3;
p = 0.5;
μ = 0;
σ = 1;
d = TransformedDistribution[
Log[x1 + 0.5] + x2, {x1 \[Distributed] NegativeBinomialDistribution[r, p],
x2 \[Distributed] LogNormalDistribution[μ, σ]}]
n = 10000;
data = RandomVariate[d, n];
ed = EmpiricalDistribution[data];


Now plot the empirical CDF and the true CDF:

 Plot[{CDF[ed, z], cdf[z, r, p, μ, σ]}, {z, -1, 10},
PlotRange -> {Automatic, {0, 1}},
PlotStyle -> {{Cyan, Thickness[0.02]}, Red},
PlotLegends -> {"Empirical", "Exact"}] It is quite simple to define a function to combine two distributions by summation e.g.

sumdistribution[dist1_, dist2_] :=
TransformedDistribution[x + y, {Distributed[x, dist1], Distributed[y, dist2]}]


For example, this gives

sumdistribution[NormalDistribution[m1, s1], NormalDistribution[m2, s2]]
(* NormalDistribution[m1 + m2, Sqrt[s1^2 + s2^2]] *)


and

sumdistribution[CauchyDistribution[m1, s1], CauchyDistribution[m2, s2]]
(* CauchyDistribution[m1 + m2, s1 + s2] *)


These results can be used in the standards ways (compute PDF, CDF, means etc or generate random numbers). e.g.

PDF[CauchyDistribution[m1 + m2, s1 + s2], x]
(* 1/(π (s1 + s2) (1 + (-m1 - m2 + x)^2/(s1 + s2)^2)) *)


There will be a lot of possible sums where a simple form is not available. However, the resulting distribution may still be useful numerically.