2
$\begingroup$

Suppose at each value of $x$ I define a normal distribution

dist[x_] := NormalDistribution[Sin[4 x], Cos[10 x]^2]

Now I want to find the distribution

$$ \mathcal{D} = \int_{-1}^1 \text{dist}(x)\ dx$$

This is essentially what TransformedDistribution is meant for but it's usually for simple combinations of distributions (e.g TransformedDistribution[ x + y, {x \[Distributed] NormalDistribution[], y \[Distributed] NormalDistribution[]}]) not integrals over a continuous set of disributions.

Of course since $\text{dist}(x)$ is always normal I know that any sum of them is also normal which allows me to compute $\mathcal{D}$ exactly, but I wondering if theres some way to get Mathemtatica to output the exact result that uses all the power of TransformedDistribution.

Thanks

Exact solution should be

NormalDistribution[Integrate[Sin[4 x], {x, -1, 1}], 
 Sqrt[Integrate[Cos[10 x]^4, {x, -1, 1}]]]

$$\text{NormalDistribution}\left[0,\frac{1}{4} \sqrt{\frac{1}{10} (120+8 \sin (20)+\sin (40))}\right] $$

$\endgroup$
6
  • 1
    $\begingroup$ I'm not understanding. You're not transforming a random variable. You're just finding a new mean and standard deviation for a normal distribution. Might you be wanting the distribution of $Y$ given that $Y|X=x$ has a normal distribution with mean Sin[4x] and standard deviation Cos[10x]^2 and $X$ has a uniform distribution on $(-1,1)$ ? $\endgroup$
    – JimB
    May 24 at 20:57
  • $\begingroup$ I think of it as a transformation which is the sum of an infinite number of random variables. Similar to the example I gave summing two normal distrubutions, but now I'm taking the limit of a Reiman sum over a continuous set. $\endgroup$ May 24 at 22:11
  • $\begingroup$ Sorry I'm being so thick about understanding your notation. Any sum of normal random variables is normal. Not disputing that. But the sum (or average) of normal pdf's is not necessarily normal. Does dist(x) represent a normal pdf? $\endgroup$
    – JimB
    May 24 at 23:01
  • $\begingroup$ @JimB, perhaps I'm also misunderstood. My goal is to integrate a function $f[x]$, but at each value of $x$ I don't know exactly what $f[x]$ is, only what the distribution is over $f[x]$ (similar to a gaussian process). What I want is the exact probability distribution for this integral. My intuition is that is integration is just a linear operation over an infinite set of gaussians so the output should itself be gaussian. $\endgroup$ May 24 at 23:53
  • $\begingroup$ @JimB: $\mathcal{D} = \int_{-1}^1 \text{dist}(x)\, dx$ is OK. As far as I understand it, such integrals arise in stochastic processes. $\endgroup$
    – user64494
    May 25 at 4:43

2 Answers 2

4
$\begingroup$

This can be done as follows. Let us consider the integral sums for $ \int_{-1}^1 \text{dist}(x)\, dx$, where the interval $[-1,1]$ is split into $n$ intervals of equal lengths and the values of $\text{dist}(x)$ are taken at its right ends:

ClearAll["Global`*"];
f[n_] := TransformedDistribution[Sum[2/n*dist[-1 + 2*j/n], 
{j, 1, n}], Table[dist[-1 + 2*j/n] \[Distributed] 
NormalDistribution[Sin[4 *(-1 + 2*j/n)], 
Cos[10 *(-1 + 2*j/n)]^2], {j, 1, n}]]

The result of

Table[f[n], {n, 1, 5}]

{NormalDistribution[2 Sin[4], 2 Cos[10]^2], NormalDistribution[Sin[4], Sqrt[1 + Cos[10]^4]], NormalDistribution[(2 Sin[4])/3, Sqrt[ 8/9 Cos[10/3]^4 + (4 Cos[10]^4)/9]], NormalDistribution[Sin[4]/2, Sqrt[1/4 + Cos[5]^4/2 + Cos[10]^4/4]], NormalDistribution[(2 Sin[4])/5, Sqrt[(8 Cos[2]^4)/25 + (8 Cos[6]^4)/25 + (4 Cos[10]^4)/25]]}

confirms that those sums are normal distributions. Now

FindSequenceFunction[Table[Mean[f[n]], {n, 1, 10}], n];

(2 Sin[4])/n

 DiscreteLimit[%,n->Infinity]

0

demonstrates that the mean of the limit distribution (which is normal according to certain theorems of probability theory) equals $0$. Unfortunately,

FindSequenceFunction[Table[Variance[f[n]], {n, 4, 14}], n]

returns the input. In view of it we have to somewhat help Mathematica. Making use of the known property of the variances of the sum of independent normal distributions, we obtain

DiscreteLimit[2/n*Sum[Variance[
NormalDistribution[Sin[4 *(-1 + 2*j/n)], 
 Cos[10 *(-1 + 2*j/n)]^2]], {j, 1, n}], n -> Infinity]

1/160 (120 + 8 Sin[20] + Sin[40])

Now

Simplify[Integrate[Cos[10 x]^4,{x,-1,1}]-%]

0

finishes the work.

$\endgroup$
4
  • $\begingroup$ Would you explain why you didn't use something like TransformedDistribution[Sum[2/n*x[j], {j, 1, n}], Table[x[j] \[Distributed] NormalDistribution[Sin[4*(-1 + 2*j/n)], Cos[10*(-1 + 2*j/n)]^2], {j, 1, n}]] for the definition of f. What advantage is there using dist[-1 + 2*j/n] in place of x[1] ? The resulting distributions are identical. $\endgroup$
    – JimB
    May 24 at 23:40
  • $\begingroup$ @JimB: I follow the notation of the question under consideration. The notation x[j]=-1 + 2*j/n is reserved for the end-points of the partition. $\endgroup$
    – user64494
    May 25 at 4:23
  • $\begingroup$ Thanks. I'm having a very slow-to-understand day. $\endgroup$
    – JimB
    May 25 at 4:27
  • $\begingroup$ @JimB: I am not a quick thinker too. $\endgroup$
    – user64494
    May 25 at 4:49
2
$\begingroup$

If (from your comment) you know that a random variable $Y$ given the value of $X=x$ has a normal distribution with mean $\sin(4x)$ and standard deviation $\cos^2(10x)$ for $-1<x<1$:

$$Y|X=x \sim N(\sin(4x),\cos^2(10x))$$

Now suppose $X\sim U(-1,1)$ (i.e., $X$ has a uniform distribution between -1 and 1). The (unconditional) probability density function for $Y=y$ is

$$p(y)=\int_{-1}^1 \frac{e^{-\frac{(y-\sin (4 x))^2}{2\cos^4(10x)}}}{\sqrt{2 \pi }\cos^2(10x)} \times\frac{1}{2}dx$$

Other than the multiplier of $1/2$, this is what you show with

$$ \mathcal{D} = \int_{-1}^1 \text{dist}(x)\ dx$$

if $\text{dist}(x)$ is the conditional pdf of $Y|X=x$. (If $\text{dist}(x)$ is something else, then I don't know what that is without further explanation from you.)

There doesn't seem to be a nice closed-form for $p(y)$. But a numerical evaluation is the following:

Plot[NExpectation[PDF[NormalDistribution[Sin[4 x], Cos[10 x]^2], y],
  x \[Distributed] UniformDistribution[{-1, 1}]], {y, -4, 4},
 AxesLabel -> {"Y", "Probability density function"}]

PDF of Y

A large random sample produces essentially the same pdf:

n = 1000000;
xx = RandomVariate[UniformDistribution[{-1, 1}], n];
parms = {Sin[4 #], Cos[10 #]^2} & /@ xx;
yy = RandomVariate[NormalDistribution[#[[1]], #[[2]]], 1][[1]] & /@ parms;
Histogram[yy, "FreedmanDiaconis", "PDF"]

Histogram of random sample from the distribution of Y

In short, it seems that you might be asking about a weighted average of pdfs rather than a weighted average of random variables.

$\endgroup$
1
  • 1
    $\begingroup$ Please, pay your attention to my comment to the question. $\endgroup$
    – user64494
    May 25 at 4:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.