Probability/Transformation of Probability Densities

A PDF version of Transformation of Probability Densities is available. 271 kb, 14 pages (info)

A Wikibook showing how to transform the probability density of a continuous random variable in both the one-dimensional and multidimensional case. In other words, it shows how to calculate the distribution of a function of continuous random variables.

This Wikibook shows how to transform the probability density of a continuous random variable in both the one-dimensional and multidimensional case.
In other words, it shows how to calculate the distribution of a function of continuous random variables.
The first section formulates the general problem and provides its solution.
However, this general solution is often quite hard to evaluate and simplifications are possible in special cases, e.g., if the random vector is one-dimensional or if the components of the random vector are independent.
Formulas for these special cases are derived in the subsequent sections.
This Wikibook also aims to provide an overview of different equations used in this field and show their connection.

First, we need to remember the definition of the cumulative distribution function, cdf, FY→(y→){\displaystyle F_{\vec {Y}}({\vec {y}})} of a random vector: It measures the probability that each component of Y takes a value smaller than the corresponding component of y.
We will use a short-hand notation and say that two vectors are "less or equal" (≤) if all of their components are.

If n=1 and m=1, X is a continuously distributed random variable with the density ϱX{\displaystyle \varrho _{X}} and f:R→R{\displaystyle f:\mathbb {R} \rightarrow \mathbb {R} } is a Borel-measurable function.
Then also Y := f(X) is continuously distributed and we are looking for the density ϱY(y){\displaystyle \varrho _{Y}(y)}.

In the following, f should always be at least differentiable.

Let us first note that there may be values that are never reached by f, e.g. y<0 if f(x) = x2.
For all those y, necessarily ϱY(y)=0{\displaystyle \varrho _{Y}(y)=0}.

At first, we limit ourselves to f those derivative is never 0 (thus, f is a diffeomorphism).
Then, the inverse map f−1{\displaystyle f^{-1}} exists and f is either monotonically increasing or monotonically decreasing.

If now the derivative f′(xi)=0{\displaystyle f^{\prime }(x_{i})=0} does vanish at some positions xi{\displaystyle x_{i}}, i=1,…,N{\displaystyle i=1,\ldots ,N}, then we split the definition space of f using those position into N+1{\displaystyle N+1} disjoint intervals Ij{\displaystyle I_{j}}.
Equation 5 holds for the functions fIj{\displaystyle f_{I_{j}}} limited in their definition space to those intervals Ij{\displaystyle I_{j}}. We have

Now we substitute u = f(x) in the integral on the right hand site, i.e. x=f−1(u){\displaystyle x=f^{-1}(u)} and dudx=f′(x){\displaystyle {\frac {du}{dx}}=f^{\prime }(x)}.
The integral limits are then from -∞ to y and in the rule “dx=dxdudu{\displaystyle dx={\frac {dx}{du}}\,du}” we have dxdu=df−1(u)du{\displaystyle {\frac {dx}{du}}={\frac {d\,f^{-1}(u)}{du}}} due to the inverse function theorem.
Consequentially:

Following the same argument as in the last section, we can again derive equation 6.

This rule often misleads physics books to present the following viewpoint, which might be easier to remember, but is not mathematically sound:
If you multiply the probability density ϱX(x){\displaystyle \varrho _{X}(x)} with the “infinitesimal length” dx{\displaystyle dx}, then you will get the probability ϱX(x)dx{\displaystyle \varrho _{X}(x)\,dx} for X to lie in the interval [x, x+dx].
Changing to new coordinates y, you will get by substitution:

The intuitive interpretation of the last expression is: one integrates over all possible x-values and uses the delta “function” to pick all positions where y = f(x).
This formula is often found in physics books, possibly written as expectation value, ⟨…⟩{\displaystyle \langle \ldots \rangle }:

Let us consider the following specific example: let ϱX(x)=exp⁡[−0.5x2]2π{\displaystyle \varrho _{X}(x)={\frac {\exp[-0.5x^{2}]}{\sqrt {2\pi }}}} and f(x)=x2{\displaystyle f(x)=x^{2}}. We choose to use equation 6 (equations 5 and 7 lead to the same result). We calculate the derivative f′(x)=2x{\displaystyle f^{\prime }(x)=2x} and find all x for which f(x)=y, which are −y{\displaystyle -{\sqrt {y}}} and +y{\displaystyle +{\sqrt {y}}} if y>0 and none otherwise. For y>0 we have:

Random numbers are generated according to the standard normal distribution X ~ N(0,1). They are shown on the x-axis of figure (a). Many of them are around 0. Each of those numbers is mapped according to y = x², which is shown with grey arrows for two example points. For many generated realisations of X, a histogram on the y-axis will converge towards the wanted probability density function ρY shown in (c). In order to analytically derive this function, we start by observing that in order for Y to be between any v and v+Δv, X must have been between either −v+Δv{\displaystyle -{\sqrt {v+\Delta v}}} and −v{\displaystyle -{\sqrt {v}}}, or, between v{\displaystyle {\sqrt {v}}} and v+Δv{\displaystyle {\sqrt {v+\Delta v}}}. Those intervals are marked in figures (b) and (c). Because the probability is equal to the area under the probability density function, we can determine ρY from the condition that the grey shaded area in (c) must be equal to the sum of the areas in (b). The areas are calculated using integrals and it is useful to take the limit Δv→0{\displaystyle \Delta v\to 0} in order to get the formula noted in (c).

If we want to obtain random numbers according to a distribution with the pdf ϱZ{\displaystyle \varrho _{Z}}, we choose f as the inverse function of the cdf of Z, i.e. Y=f(X)=FZ−1(X){\displaystyle Y=f(X)=F_{Z}^{-1}(X)}. We can now show that Y will have the same distribution as the wanted Z, ϱY=ϱZ{\displaystyle \varrho _{Y}=\varrho _{Z}} by using equation 5 and the fact that f−1=FZ{\displaystyle f^{-1}=F_{Z}}:

Random numbers yi are generated from a uniform distribution between 0 and 1, i.e. Y ~ U(0, 1). They are sketched as colored points on the y-axis. Each of the points is mapped according to x=F-1(y), which is shown with gray arrows for two example points. In this example, we have used an exponential distribution. Hence, for x ≥ 0, the probability density is ϱX(x)=λe−λx{\displaystyle \varrho _{X}(x)=\lambda e^{-\lambda \,x}} and the cumulated distribution function is F(x)=1−e−λx{\displaystyle F(x)=1-e^{-\lambda \,x}}. Therefore, x=F−1(y)=−ln⁡(1−y)λ{\displaystyle x=F^{-1}(y)=-{\frac {\ln(1-y)}{\lambda }}}. We can see that using this method, many points end up close to 0 and only few points end up having high x-values - just as it is expected for an exponential distribution.

We will now investigate the case when a random vector X with known density ϱX→{\displaystyle \varrho _{\vec {X}}} is mapped to (scalar) random variable Y and calculate the new density ϱY(y){\displaystyle \varrho _{Y}(y)}.

The direct evaluation of this equation is sometimes the easiest way, e.g., if there is a known formula for the area or volume presented by the integral.
Otherwise one needs to solve a parameter-depended multiple integral.

If the components of the random vector X→{\displaystyle {\vec {X}}} are independent, then the probability density factorizes:

Figure (a) shows a sketch of random sample points from a uniform distribution in a circle of radius 1, subdivided into rings. In figure (b), we count how many of those points fell into each of the rings with the same width Δv. Since the area of the rings increases linearly with the radius, one can expect more points for larger radii. For Δv → 0, the normalized histogram in (b) will converge to the wanted probability density function ρY. In order to calculate ρY analytically, we first derive the cumulated distribution function FY, plotted in figure (d). FY(y) is the probability to find a point inside the circle of radius v (shown in grey in figure (c)). For v between 0 and 1, we find FY(v)=AvAtot=πv2π12=v2{\displaystyle F_{Y}(v)={\frac {A_{v}}{A_{tot}}}={\frac {\pi v^{2}}{\pi 1^{2}}}=v^{2}}. The slope of FY is the wanted probability density function ρY(y)=dFY(y)dy=2y{\displaystyle \rho _{Y}(y)={\frac {dF_{Y}(y)}{dy}}=2y}, in agreement with figure (b).

If f is injective, then it can be easier to introduce additional helper coordinates Ym+1 to Yn, then do the Rn→Rn{\displaystyle \mathbb {R} ^{n}\to \mathbb {R} ^{n}} transformation from section Invertible Transformation of a Random Vector and finally integrate out all helper coordinates of the so-obtained density.

Remark: If the joint pdf ϱY→(y1,y2){\displaystyle \varrho _{\vec {Y}}(y_{1},y_{2})}, i.e. the conditional distribution, is not of interest and one is only interested in the marginal distribution with ϱY1(y1)=∫RϱY→(y1,y2)dy2{\displaystyle \varrho _{Y_{1}}(y_{1})=\int _{\mathbb {R} }\varrho _{\vec {Y}}(y_{1},y_{2})\,dy_{2}}, then it is possible to calculate the density as described in section Mapping of a Random Vector to a Random Variable for the mapping Y1 = 1 X1 + 2 X2 + 3 X3 (and likewise for Y2 = 4 X1 + 5 X2 + 6 X3).

In order to show some possible applications, we present the following questions, which can be answered using the techniques outlines in this Wikibook. In principle, the answers could also be approximated using a numerical random number simulation: generate several realizations of X→{\displaystyle {\vec {X}}}, calculate Y→=f(X→){\displaystyle {\vec {Y}}=f({\vec {X}})} and make a histogram of the results. However, many such random numbers are needed for reasonable results, especially for higher-dimensional random vectors. Gladly, we can always calculate the resulting distribution analytically using the above formulas.

Suppose the velocity components of an ideal gas (Vx, Vy, Vz) are identically, independently normally distributed as in the last example. What is the probability density ϱV{\displaystyle \varrho _{V}} of V=Vx2+Vy2+Vz2{\displaystyle V={\sqrt {V_{x}^{2}+V_{y}^{2}+V_{z}^{2}}}}? (The answer is known as Maxwell-Boltzmann distribution.)

Suppose we do not know the exact value of X and Y, but we can assign a probability distribution to each of them. What is the distrubution of the derived property Z = X2 / Y and what is the mean value and standard deviation of Z? (To tackle such problems, linearisation around the mean values are sometimes used and both X and Y are assumed to be normally distributed. However, we are not limited to such restrictions.)

Suppose we consider the value of one gram gold, silver and platinum in one year from now as independent random variables G, S and P, respectively. Box A contains 1 gram gold, 2 gram silver and 3 gram platinum. Box B contains 4, 5 and 6 gram, respectively. Thus, (AB)=(123456)(GSP){\displaystyle {\begin{pmatrix}A\\B\end{pmatrix}}={\begin{pmatrix}1&2&3\\4&5&6\end{pmatrix}}{\begin{pmatrix}G\\S\\P\end{pmatrix}}}. What is the value of the contents in box A (or box B) in one year from now? (The answer is given in an example above.) Note that A and B are correlated.

Note that the above examples assume the distribution of X→{\displaystyle {\vec {X}}} to be known. If it is unknown, or if the calculation is based on only a few data points, methods from mathematical statistics are a better choice to quantify uncertainty.

Correlated random numbers can be obtained by first generating a vector of uncorrelated random numbers and then applying a function on them.

In order to obtain random numbers with covariance matrix CY, we can use the following know procedure: Calculate the Cholesky decomposition CY = A AT. Generate a vector x→{\displaystyle {\vec {x}}} of uncorrelated random numbers with all var(Xi) = 1. Apply the matrix A: Y→=AX→{\displaystyle {\vec {Y}}=A{\vec {X}}}. This will result in correlated random variables with covariance matrix CY = A AT.

With the formulas outlined in this Wikibook, we can additionally study the shape of the resulting distribution and the effect of non-linear transformations. Consider, e.g., that X is uniform distributed in [0, 2π], Y1 = sin(X) and Y2 = cos(X). In this case, a 2D plot of random numbers from (Y1, Y2) will show a uniform distribution on a circle. Although Y1 and Y2 are stochastically dependent, they are uncorrelated. It is therefore important to know the resulting distribution, because ϱY→(y1,y1){\displaystyle \varrho _{\vec {Y}}(y_{1},y_{1})} has more information than the covariance matrix CY.