.. _random_mixture: Random Mixture: affine combination of independent univariate distributions -------------------------------------------------------------------------- A multivariate random variable :math:\vect{Y} may be defined as an affine transform of :math:n independent univariate random variable, as follows: .. math:: :label: RandomMixtureFormula \displaystyle \vect{Y}=\vect{y}_0+\mat{M}\,\vect{X} where :math:\vect{y}_0\in\mathbb{R}^d is a deterministic vector with :math:d\in\{1,2,3\}, :math:\mat{M}\in\mathcal{M}_{d,n}(\mathbb{R}) a deterministic matrix and :math:(X_k)_{ 1 \leq k \leq n} are some independent univariate distributions. In such a case, it is possible to evaluate directly the distribution of :math:\vect{Y} and then to ask :math:\vect{Y} any request compatible with a distribution: moments, probability and cumulative density functions, quantiles (in dimension 1 only) ... **Evaluation of the probability density function of the Random Mixture** As the univariate random variables :math:X_i are independent, the characteristic function of :math:\vect{Y}, denoted :math:\phi_Y, is easily defined from the characteristic function of :math:X_k denoted :math:\phi_{X_k} as follows : .. math:: :label: CharactFuncY \displaystyle \phi_Y(u_1,\hdots,u_d)=\prod_{j=1}^de^{iu_j{y_0}_j}\prod_{k=1}^n\phi_{X_k}((M^tu)_k), \mbox{ for } \vect{u}\in\mathbb{R}^d | Once :math:\phi_Y evaluated, it is possible to evaluate the probability density function of :math:Y, denoted :math:p_Y : several techniques are possible, as the inversion of the Fourier transformation. This technique is not easy to implement. | Another technique is used, based on the Poisson sum formulation, defined as follows: .. math:: :label: PoissonSum \displaystyle \sum_{j_1\in\mathbb{Z}}\hdots\sum_{j_d\in\mathbb{Z}} p_Y\left(y_1+\frac{2\pi j_1}{h_1},\hdots,y_d+\frac{2\pi j_d}{h_d}\right)= \prod_{j=1}^d \frac{h_j}{2*\pi}\sum_{k_1\in\mathbb{Z}}\hdots\sum_{k_d\in\mathbb{Z}}\phi\left(k_1h_1,\hdots,k_dh_d\right)e^{-\imath(\sum_{m=1}^{d}k_m h_m y_m)} | By fixing :math:h_1,\hdots,h_d small enough, :math:\frac{2k\pi}{h_j} \approx +\infty and :math:p_Y(\hdots,\frac{2k\pi}{h_j},\hdots) \approx 0 because of the decreasing properties of :math:p_Y. Thus the nested sums of the left term of :eq:PoissonSum are reduced to the central term :math:j_1=\hdots=j_d = 0: the left term is approximatively equal to :math:p_Y(y). | Furthermore, the right term of :eq:PoissonSum is a series which converges very fast: only few terms of the series are enough to get machine-precision accuracy. Let us note that the factors :math:\phi_Y(k_1 h_1,\hdots,k_d,h_d), which are expensive to evaluate, do not depend on :math:y and are evaluated once only. | It is also possible to greatly improve the performance of the algorithm by noticing that equation  is linear between :math:p_Y and :math:\phi_Y. We denote :math:q_Y and :math:\psi_Y respectively the density and the characteristic function of the multivariate normal distribution with the same mean :math:\vect{\mu} and same covariance matrix :math:\vect{C} as the random mixture. By applying this multivariate normal distribution to the equation , we obtain by subtraction: .. math:: :label: algoPoisson \displaystyle p_Y\left(y\right) = \sum_{j\in\mathbb{Z}^d} q_Y\left(y_1+\frac{2\pi j_1}{h_1},\cdots,y_d+\frac{2\pi j_d}{h_d}\right)+ \frac{H}{2^d\pi^d}\sum_{|k_1|\leq N}\cdots\sum_{|k_d|\leq N} \delta_Y\left(k_1h_1,\cdots,k_dh_d\right)e^{-\imath(\sum_{m=1}^{d}k_m h_m y_m)} where :math:H = h_1\times\cdots\times h_d, :math:j=(j_1,\cdots,j_d), :math:\delta_Y:=\phi_Y - \psi_Y | In the case where :math:n \gg 1, using the limit central theorem, the law of :math:\vect{Y} tends to the normal distribution density :math:q, which will drastically reduce :math:N. The sum on :math:q will become the most CPU-intensive part, because in the general case we will have to keep more terms than the central one in this sum, since the parameters :math:h_1, \dots h_d were calibrated with respect to :math:p and not :math:q. The parameters :math:h_1, \dots h_d are calibrated using the following formula: .. math:: h_\ell = \frac{2\pi}{(\beta+4\alpha)\sigma_\ell} where :math:\sigma_\ell=\sqrt{\Cov{\vect{Y}}_{\ell,\ell}} and :math:\alpha, :math:\beta are respectively the number of standard deviations covered by the marginal distribution (:math:\alpha=5 by default) and :math:\beta the number of marginal deviations beyond which the density is negligible (:math:\beta=8.5 by default). The :math:N parameter is dynamically calibrated: we start with :math:N=8 then we double :math:N value until the total contribution of the additional terms is negligible. **Evaluation of the moments of the Random Mixture** The relation :eq:RandomMixtureFormula enables to evaluate all the moments of the random mixture, if mathematically defined. For example, we have: .. math:: \left\{ \begin{array}{lcl} \Expect{\vect{Y}} & = & \vect{y_0} + \mat{M}\Expect{\vect{X}} \\ \Cov{\vect{Y}} & = & \mat{M}\,\Cov{\vect{X}}\mat{M}^t \end{array}\right\} **Computation on a regular grid** The interest is to compute the density function on a regular grid. Purposes are to get an approximation quickly. The regular grid is of form: .. math:: \begin{aligned} \forall r\in\{1,\hdots,d\},\forall m\in\{0,\hdots,M-1\},\:y_{r,m}=\mu_r+b\left(\frac{2m+1}{M} - 1\right)\sigma_r \end{aligned} By denoting :math:p_{m_1,\hdots,m_d}=p_{\vect{Y}}(y_{1,m_1},\hdots,y_{d,m_d}): .. math:: \begin{aligned} p_{m_1,\hdots,m_d}= Q_{m_1,\hdots,m_d}+S_{m_1,\hdots,m_d} \end{aligned} for which the term :math:S_{m_1,\hdots,m_d} is the most CPU consuming. This term rewrites: .. math:: \begin{aligned} S_{m_1,\hdots,m_d}=&\frac{H}{2^d\pi^d}\sum_{k_1=-N}^{N}\hdots\sum_{k_d=-N}^{N}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) \label{Eq:S} \end{aligned} with: .. math:: \begin{aligned} \delta\left(k_1h_1,\hdots,k_dh_d\right)&=(\phi-\psi)\left(k_1h_1,\hdots,k_dh_d\right)\\ E_{m_1,\hdots,m_d}(k_1,\hdots,k_d)&=e^{-i\sum_{j=1}^d k_jh_j\left(\mu_j+b\left(\frac{2m_j+1}{M}-1\right)\sigma_j\right)} \end{aligned} The aim is to rewrite the previous expression as a :math:d- discrete Fourier transform, in order to apply Fast Fourier Transform (*FFT*) for its evaluation. We set :math:M=N and :math:\forall j \in\{1,\hdots,d\},\: h_j=\frac{\pi}{b\sigma_j} and :math:\tau_j=\frac{\mu_j}{b\sigma_j}. For convenience, we introduce the functions: .. math:: f_j(k) = e^{-i\pi (k+1)\left(\tau_j-1+\frac{1}{N}\right)} We use :math:k+1 instead of :math:k in this function to simplify expressions below. We obtain: .. math:: \begin{aligned} E_{m_1,\hdots,m_d}(k_1,\hdots,k_d)&=e^{-i\sum_{j=1}^{d} k_jh_jb\sigma_j\left(\frac{\mu_j}{b\sigma_j}+\frac{2m_j}{N}+\frac{1}{N}-1\right)}\notag\\ &=e^{-2i\pi\left(\frac{\sum_{j=1}^{d}k_j m_j}{N}\right)}e^{-i\pi\sum_{j=1}^{d} k_j\left(\tau_j-1+\frac{1}{N}\right)} \notag\\ &=e^{-2i\pi\left(\frac{\sum_{j=1}^{d}k_j m_j}{N}\right)} f_1(k_1-1) \times\hdots\times f_d(k_d-1) \label{Eq:E} \end{aligned} For performance reasons, we want to use the discrete Fourier transform with the following convention in dimension 1: .. math:: A_m = \sum_{k=0}^{N-1} a_k e^{-2i\pi\frac{km}{N}} which extension to dimensions 2 and 3 are respectively: .. math:: A_{m,n} = \sum_{k=0}^{N-1}\sum_{l=0}^{N-1} a_{k,l} e^{-2i\pi\frac{km}{N}} e^{-2i\pi\frac{ln}{N}}\\ .. math:: A_{m,n,p} = \sum_{k=0}^{N-1}\sum_{l=0}^{N-1}\sum_{s=0}^{N-1} a_{k,l,s} e^{-2i\pi\frac{km}{N}} e^{-2i\pi\frac{ln}{N}} e^{-2i\pi\frac{sp}{N}} We decompose sums of  on the interval :math:[-N,N] into three parts: .. math:: :label: decomposition-sum \begin{aligned} \sum_{k_j=-N}^{N}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) = & \sum_{k_j=-N}^{-1} \delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) \notag\\ & + \delta\left(k_1h_1,\hdots,0,\hdots,k_dh_d\right) E_{m_1,\hdots,0,\hdots,m_d}(k_1,\hdots,0,\hdots,k_d) \notag\\ & + \sum_{k_j=1}^{N}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) \end{aligned} If we already computed :math:E for dimension :math:d-1, then the middle term in this sum is trivial. To compute the last sum of equation, we apply a change of variable :math:k_j'=k_j-1: .. math:: \begin{aligned} \sum_{k_j=1}^{N}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) = & \sum_{k_j=0}^{N-1}\delta\left(k_1h_1,\hdots,(k_j+1)h_j,\hdots,k_dh_d\right) \times\notag\\ & \hspace*{3cm} E_{m_1,\hdots,m_d}(k_1,\hdots,k_j+1,\hdots,k_d) \end{aligned} Equation gives: .. math:: \begin{aligned} E_{m_1,\hdots,m_d}(k_1,\hdots,k_j+1,\hdots,k_d) &= e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N} +\frac{m_j}{N}\right)} f_1(k_1-1)\times\hdots\times f_j(k_j)\times\hdots\times f_d(k_d-1)\notag\\ &= e^{-2i\pi\left(\frac{m_j}{N}\right)} e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N}\right)} f_1(k_1-1)\times\hdots\times f_j(k_j)\times\hdots\times f_d(k_d-1) \end{aligned} Thus .. math:: \begin{aligned} \sum_{k_j=1}^{N}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}&(k_1,\hdots,k_d) = e^{-2i\pi\left(\frac{m_j}{N}\right)} \sum_{k_j=0}^{N-1}\delta\left(k_1h_1,\hdots,(k_j+1)h_j,\hdots,k_dh_d\right) \times\notag\\ & e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N}\right)} f_1(k_1-1)\times\hdots\times f_j(k_j)\times\hdots\times f_d(k_d-1) \label{Eq:j-sigma+} \end{aligned} To compute the first sum of equation, we apply a change of variable :math:k_j'=N+k_j: .. math:: \begin{aligned} \sum_{k_j=-N}^{-1}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}(k_1,\hdots,k_d) = & \sum_{k_j=0}^{N-1}\delta\left(k_1h_1,\hdots,(k_j-N)h_j,\hdots,k_dh_d\right) \times\notag\\ & \hspace*{3cm} E_{m_1,\hdots,m_d}(k_1,\hdots,k_j-N,\hdots,k_d) \end{aligned} Equation  gives: .. math:: \begin{aligned} E_{m_1,\hdots,m_d}(k_1,\hdots,k_j-N,\hdots,k_d) &= e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N} -m_j\right)} f_1(k_1-1)\times\hdots\times f_j(k_j-1-N)\times\hdots\times f_d(k_d-1) \notag\\ &= e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N}\right)} f_1(k_1-1)\times\hdots\times \overline{f}_j(N-1-k_j)\times\hdots\times f_d(k_d-1) \end{aligned} Thus: .. math:: \begin{aligned} \sum_{k_j=-N}^{-1}\delta\left(k_1h_1,\hdots,k_dh_d\right) E_{m_1,\hdots,m_d}&(k_1,\hdots,k_d) = \sum_{k_j=0}^{N-1}\delta\left(k_1h_1,\hdots,(k_j-N)h_j,\hdots,k_dh_d\right) \times\notag\\ & e^{-2i\pi\left(\frac{\sum_{l=1}^{d}k_l m_l}{N}\right)} f_1(k_1-1)\times\hdots\times \overline{f}_j(N-1-k_j)\times\hdots\times f_d(k_d-1) \label{Eq:j-sigma-} \end{aligned} To summarize: #. In order to compute sum from :math:k_1=1 to :math:N, we multiply by :math:e^{-2i\pi\left(\frac{m_1}{N}\right)} and consider :math:\delta((k_1+1)h,\hdots)f_1(k_1) #. In order to compute sum from :math:k_1=-N to :math:-1, we consider :math:\delta((k_1-N)h,\hdots)\overline{f}_1(N-1-k_1) .. topic:: API: - See :class:~openturns.RandomMixture .. topic:: Examples: - See :doc:/auto_probabilistic_modeling/distributions/plot_random_mixture_distribution - See :doc:/auto_probabilistic_modeling/distributions/plot_random_mixture_distribution_discrete .. topic:: References: - [abate1992]_