I am used to applying a specific change of variable method in the case of univariate quantity (see Univariate case section below).
I have difficulty transposing this method to the multivariate case (see Multivariate case section below), despite the solution being analog to the univariate case.
Is my approach a dead end? In particular problem seems to result from mapping to a lower dimensional space (see Discussion section below).
Note that the results I am seeking are well known in both cases (e.g. http://www2.stat.duke.edu/~st118/sta732/mvnormal.pdf, theorem 2). I am simply curious to find out if and how the method used for univariate case can be made to work in the multivariate setting...
Univariate case
Consider the following univariate random variable \begin{align} x &\sim \mathcal{N}(\mu, \sigma^2),\\ x &\in \mathbb{R}, \end{align} and the following relation \begin{align} y &= wx \end{align} where $w$ denotes a scalar weight value (deterministic).
Given the above, $p(y)$ can be derived as follows: \begin{align} p(y) &= \int_{\mathbb{R}} p(y, x)\ \text{d}x\\ &= \int_\mathbb{R} p(y| x)p(x)\ \text{d}x \\ &= \int_\mathbb{R} \delta(y - wx)\mathcal{N}(x|\mu, \sigma^2)\ \text{d}x\\ &\propto \int_\mathbb{R} \delta(wx - y)\exp\left(-\frac{1}{2}\frac{1}{\sigma^2} (x-\mu)^2 \right)\ \text{d}x\\ &\qquad \begin{array}{c c} v = wx & x = \frac{v}{w}\\ \text{d} v = w\text{d}x & \text{d}x = \frac{1}{w}\text{d}v \end{array}\\ p(y) &\propto \frac{1}{w}\int_\mathbb{R} \delta(v - y)\exp\left(-\frac{1}{2}\frac{1}{\sigma^2} \left(\frac{v}{w}-\mu\right)^2 \right)\ \text{d} v\\ &\propto \int_{\mathbb{R}} \delta(v - y)\exp\left(-\frac{1}{2}\frac{1}{\sigma^2w^2} \left(v-w\mu\right)^2 \right)\ \text{d} v\\ &= \exp\left(-\frac{1}{2}\frac{1}{\sigma^2w^2} \left(y-w\mu\right)^2 \right) \end{align} from which we obtain the well known result \begin{align} y &\sim \mathcal{N}\left(w\mu, w^2\sigma^2\right). \end{align}
Multivariate case
$\mathbf{x}$ is now a multivariate random variable \begin{align} \mathbf{x} &\sim \mathcal{N}(\mathbf{\mu}, \Sigma),\\ \mathbf{x} &\in \mathbb{R}^d,\qquad d \in \{2, 3, 4, ...\}, \end{align} and \begin{align} y &= \mathbf{w}^\text{T}\mathbf{x} \end{align} where $\mathbf{w}$ denotes a weight vector (as in univariate case, this quantity is deterministic).
Can we demonstrate the following result using a change of variable method similar to the one given for the univariate case? \begin{align} y &\sim \mathcal{N}(\ \mathbf{w}\mathbf{\mu}, \ \mathbf{w}^\text{T}\Sigma\mathbf{w}) \end{align}
My attempt so far \begin{align} p(y) &= \int_{\mathbb{R}^d} p(y, \mathbf{x})\ \underbrace{\text{d}\mathbf{x}}_{\text{d}x_1\text{d}x_2,\ ...,\ \text{d}x_d}\\ &= \int_{\mathbb{R}^d} p(y| \mathbf{x})p(\mathbf{x})\ \text{d}\mathbf{x} \\ &= \int_{\mathbb{R}^d} \delta(y - \mathbf{w}^\text{T}\mathbf{x})\mathcal{N}(\mathbf{x}|\mathbf{\mu}, \Sigma)\ \text{d}\mathbf{x}\\ &\propto \int_{\mathbb{R}^d} \delta(\mathbf{w}^\text{T}\mathbf{x} - y)\exp\left(-\frac{1}{2} (\mathbf{x}-\mathbf{\mu})^\text{T}\Sigma^{-1} (\mathbf{x}-\mathbf{\mu}) \right)\ \text{d}\mathbf{x}\\ &\qquad \begin{array}{c c} v = \mathbf{w}^\text{T}\mathbf{x} & \mathbf{x} = ?\\ \text{d} v = ? & \text{d}x = ? \end{array}\\ \end{align}
Discussion
It seems that the issues comes results from the mapping $y = \mathbf{w}^\text{T}\mathbf{x}$ not being bijective (not having an inverse).
According to Betancourt (https://betanalpha.github.io/assets/case_studies/probability_theory.html, section 4.2):
Deriving the pushforward probability density function for transformations that change the dimensionality of the space, such as marginalizations, are more challenging and require analytically integrating an appropriately reparameterized probability density function over the complementary spaces.
Unfortunately, I am not sure what ``complementary spaces'' mean here, and how to achieve this in practice...
A similar question was asked in Change of variables to lower-dimensional space
The matter is quite simple when you work with random variables rather than the densities: $X$ is multivariately normally distributed with mean vector $\mu$ and covariance matrix $\Sigma$. To avoid pathological cases let's make the technical assumption that the vector $X$ is built from $d$ uncorrelated normal variables $Z_1,...,Z_d$ by $$ X=LZ $$ where $LL^\top=\Sigma$ is the Cholesky decomposition. This ensures that $$ Y=\boldsymbol{w}^TX=w_1X_1+...+w_dX_d $$ is univariately normally distributed (because it is effecively a linear combination of uncorrelated normals $Z_i$). The mean of $Y$ is clearly $$ \mu_Y=w_1\mu_1+...+w_d\mu_d=\boldsymbol{w}^T\mu. $$ The variance of $Y$ is clearly $$ \sum_{i,j=1}^d w_iw_j\Sigma_{ij}=\boldsymbol{w}^\top\Sigma\boldsymbol{w}\,. $$ This shows $$ Y\sim{\cal N}(\boldsymbol{w}^T\mu,\boldsymbol{w}^\top\Sigma\boldsymbol{w}). $$