# How to perform convolutions on this complicated pdf?

For a known pdf f[x], I need to calculate integrals over the probability density of the sum of Log[f[x]] over independent samples:

$$Q[t,n]\equiv P(\sum_{i=1}^{n}Log[f[x_{i}]]=t)$$

Since the sample terms in the sum are independent, Q[t,n] is just the pdf of values of the individual terms convolved n times: $$Q[t,n]\equiv q[t]*Q[t,n-1]$$ $$q[t]\equiv P(Log[f]=t)$$

But the pdf we need to self-convolve is difficult to represent even numerically:

$$q[t]=\sum_{\{x;Log[f[x]]=t\}}\frac{f[x]}{|\frac{dLog[f[x]]}{dx}|}=e^{2t}\sum_{\{x;f[x]=e^{t}\}}|\frac{df[x]}{dx}|^{-1}$$

The sum is over all the values of x for which Log[f[x]]=t, so calculating q[y] exactly would involve inverting $f[x]\rightarrow x[f]$ which is not generally possible. Also note the complicating absolute value in the denominator.

My Question is: For a known f[x], how can I represent q[y] in MMa so that I can calculate integrals of Q[n,t] with known Precision?

$$R_{n}[t_{1},t_{2}]\equiv P(t_{1}\leq\sum_{i=1}^{n}Log[f[x_{i}]]\leq t_{2})=\int_{t_{1}}^{t_{2}}Q[t,n]dt$$

EXAMPLE: Monte Carlo integration

My f[x]'s will almost always be Gaussian Mixtures, so consider a two-term example. The histograms below numerically approximate f[x] and q[t]. Note that the histogram of t has 3 singularities, corresponding to the 3 values of x where df/dx=0 and to the Abs[df/dx] in q[t]'s denominator.

gau[x_, v_] := Exp[-(x^2)/(2*v)]/Sqrt[2*Pi*v];
f[x_] := (2/3)*gau[x + 3, 2] + (1/3)*gau[x - 2, 2];
Plot[f[x], {x, -8, 8}]
p = ProbabilityDistribution[f[x], {x, -Infinity, Infinity}];
bigNumber=10^5
xx = RandomVariate[p, bigNumber];
Histogram[xx, 100]
t=Log[f[xx]];
Histogram[Select[t, # > -4 &], 100]


So a brute force way to approximate $R_{n}[t_{1},t_{2}]$ would be to simply generate many sets of n random x's and count how often the sum of logs of f[x] fell in the specified range.

R[n_, t1_, t2_] := Length[Select[Total[Log[f[RandomVariate[p, {n, bigNumber}]]]], (t1 <= #) && (# <= t2) &]]/bigNumber // N

R[3,-7,-5]

(*  0.59384  *)


But this method is obviously a bit slow and imprecise. Is there a more practical way?

• Just a thought: does The Convolution Theorem help at all? – jackskis Jun 22 '16 at 2:05
• Deriving an analytic solution looks to be quite hard. The best approach may depend on how many samples you need to sum, how far down the tails you need to integrate, and the accuracy you require. When summing large numbers of samples, without going too far down the tail, the central limit theorem might give a good enough approximation. Monte Carlo integration should give you reliable error bounds and is easy to program, but you may not have time to wait for the results. If you are really interested in the lower tail probabilities, I guess that specific analysis of that might be effective. – mikado Jun 27 '16 at 21:21
• @mikado Yep, that's pretty much the conclusions I came to. I posted the Question because I hoped someone knew some cool MMa wizardry I didn't know about. – Jerry Guern Jun 27 '16 at 22:21
• @mikado BTW, the CLT is a no-go because my n's are small and this calculation is unnecessary for large n. – Jerry Guern Jun 28 '16 at 2:00

## 1 Answer

I offer the following based on Monte Carlo integration, but trying to make use of the facilities Mathematica provides:

weights = {2/3, 1/3};
means = {-3, 2};
sd = Sqrt[2];
distributions = NormalDistribution[#, sd] & /@ means;
fdist = MixtureDistribution[weights, distributions];
qdist1 = TransformedDistribution[Log[PDF[fdist, f]],
Distributed[f, fdist]];
ntries = 10^5;
data[n_] := data[n] = Total /@ RandomVariate[qdist1, {ntries, n}]
qdist[n_] := EmpiricalDistribution[data[n]]
R[n_, t1_, t2_] := CDF[qdist[n], t2] - CDF[qdist[n], t1]
Rerr[n_, t1_, t2_] :=
Module[{p = R[n, t1, t2]}, Sqrt[p (1 - p)/ntries]]


We first define weights for a MixtureDistribution of NormalDistributions. We define the desired distribution of Q in terms of a TransformedDistribution. We cache random data generated from this distribution, create an EmpiricalDistribution object, then define the required probabilities R based on its CDF. We generate an estimate of the standard deviation of the estimated probability Rerr using results for the Binomial distribution.

ntries could be chosen to achieve the required accuracy, though this might be a problem in the tails of the distribution.

If a large number of values of R need to be computed, it might be worth using FunctionInterpolation to create a more rapidly evaluated version of CDF[qdist[n],t].

The EmpiricalDistribution might be replaced by other DataDistribution objects, e.g. HistogramDistribution would be better for looking at the PDF rather than CDF.

I get results as follows

Plot[CDF[qdist[1], x], {x, -6, 1}]


and

Plot[CDF[qdist[5], x], {x, -20, -8}]


• 1) This is a BIG improvement, the same basic method but you've demo'd several built-ins I didn't know about. Thank you. 2) I'm still hoping for a non-Monte Carlo method because of speed and precision. 3) Was there some strategic reason you calculated q[n_] by summing over n terms instead of calculating q[1] and convolving n times? – Jerry Guern Jun 29 '16 at 0:56
• @JerryGuern Point 2) Can you quantify your need for speed and precision? Point 3) a. The empirical CDF can be thought of as an enormous Piecewise expression. Convolving it with itself will increase the number of parts from n to n^2 (approximately). b. Calculating the random numbers involves comparatively trivial computation and only needs to be done once. c. As the distribution has finite variance, summing a number of terms will tend to make the distribution better behaved. – mikado Jun 29 '16 at 18:16
• Part 2) I can't really quantify, I was just hoping for something fast so I can do thinks like like optimize R[n,t1,t2] over parameters of f[x]. And when one tries to optimize using Monte Carlo, the result is affected by the statistical variation, so knowing the precision is helpful. Part 3) Ah, good points, thanks for that explanation. – Jerry Guern Jun 30 '16 at 21:06
• I've discovered that CDF and PDF of qdist1 can be evaluated for numerical values (though not very quickly). There may be a route to solve this by using FunctionInterpolation to find a rapidly evaluable approximation to the PDF. A numerical convolution might then be used. – mikado Jun 30 '16 at 21:58