Suppose $X\sim \bf{N}_m(\bf{\mu}, \bf{\Sigma})$ where the covariance matrix $\bf \Sigma$ is positive definite. Show that the distribution of $X^T{\bf{B}}X \sim \chi^2(\bf{tr(B\Sigma)})$where $\bf B$ satisfies    $\bf{B\Sigma B}=\bf{B}$ $\bf\mu^TB\mu=0$  Now, previously I had done this result, when $\Sigma$ was replaced by $\bf{I}$ and $B$ was idempotent. This one is where I am stuck. I have found that $\bf B\Sigma$ and $\bf \Sigma B$ are idempotent. Now, I tried to find some $U\sim \bf{N_{\text{tr}(B\Sigma)}}(0, I)$ such that $U^TU=X^T{\bf{B}}X$ then we would be done, but I can't seem to find such an $U$. Can someone help me? Thanks.
123557
Suppose $\underline{X}$ is an $m$-dimensional vector following multivariate Normal distribution i.e. $\underline{X}$~$N_m(\underline{\mu},\Sigma)$ where $\Sigma$ is positive definite. Let $B$ be a symmetric $m$ by $m$ matrix, having the following properties:  $B\Sigma B=B$  $\underline{\mu}'B\underline{\mu}=0$  $trace(B\Sigma)=k$, where $k&lt;m$   Prove that $\underline{X}'B\underline{X}$ follows chi-squared ($k$).  I proceeded a lot but got stuck in one last step where I think linear algebra is troubling me. Kindly have  a look at my &quot;solution&quot;.  The second condition can be written as $\underline{\mu}'B\Sigma B\underline{\mu}=0\implies B\underline{\mu}=\underline{0}$ as $\Sigma$ is Positive definite. Then after a little algebra and considering condition 2, we can write $\underline{X}'B\underline{X}=(\underline{X}-\underline{\mu})'B(\underline{X}-\underline{\mu})=(\underline{X}-\underline{\mu})'B\Sigma B(\underline{X}-\underline{\mu})$ But $\Sigma$ is positive definite hence we can have a symmetric matrix $\sqrt{\Sigma}$ defined such that $\sqrt{\Sigma}\sqrt{\Sigma}=\Sigma$. Using this, we can write $(\underline{X}-\underline{\mu})'B\Sigma B(\underline{X}-\underline{\mu})=(\sqrt{\Sigma}B(\underline{X}-\underline{\mu}))'(\sqrt{\Sigma}B(\underline{X}-\underline{\mu}))$ We want to show that this follows chi-squared ($k$). Now, it is obvious that as $(\underline{X}-\underline{\mu})$~$N(\underline{0},\Sigma)$, $\sqrt{\Sigma}B(\underline{X}-\underline{\mu})$~$N(\underline{0},\sqrt{\Sigma}B\Sigma B\sqrt{\Sigma})=N(\underline{0},\sqrt{\Sigma}B\sqrt{\Sigma})$ Now, $trace(\sqrt{\Sigma}B\sqrt{\Sigma})=trace(B\Sigma)=k$ and $\sqrt{\Sigma}B\sqrt{\Sigma}$ is idempotent (verified). So its eigenvalues are only $0$ or $1$. Since its trace is $k$, there are $k$ eigenvalues with value $1$ and other $m-k$ eigenvalues are $0$. Now here is the thing I want to prove: I want to show that$\sqrt{\Sigma}B\sqrt{\Sigma}$, which is the covariance matrix, is diagonal and has only $k$ $1$'s and remaining $m-k$ diagonal entries $0$. The reason is, then I will be able to say that $k$ components of $\sqrt{\Sigma}B(\underline{X}-\underline{\mu})$ are independent standard normals whose sum of squares is our quadratic form and thus it follows chi-squared ($k$). But I really could not show that $\sqrt{\Sigma}B\sqrt{\Sigma}$ has the desired property i.e. it is diagonal with rank $k$. But this is crucial, it is the finale. If this does not happen, then some components will be correlated with each other and we cannot say that our quadratic form has chi-squared ($k$). I would really appreciate if someone helped me to prove this final part.
34689
Definition. Suppose $\mathbf{y} \sim \mathcal{N}(\boldsymbol{\mu}, I_{n \times  n})$.      Then $$w = \mathbf{y}^{T}\mathbf{y} = \|\mathbf{y}\|^2 \sim \chi^{2}_{n}\left(\theta = \|\boldsymbol{\mu}\|^2/2 =\boldsymbol{\mu}^{T}\boldsymbol{\mu}/2 \right)\text{,}$$   i.e., the $\chi^2$ distribution with $n$ degrees of freedom and   noncentrality parameter $\theta$.   I am trying to (at the very least) find a proof of the following theorem:     Theorem. Suppose $\Sigma$ is a $n \times n$ positive definite and symmetric matrix, $A$ is a $n \times n$ symmetric matrix with rank $m$, and $(A\Sigma)^2 = A\Sigma$.      Then $\mathbf{y} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma) \implies \mathbf{y}^{T}A\mathbf{y} \sim \chi^2_m\left(\boldsymbol{\mu}^{T}A\boldsymbol{\mu}/2\right)\text{.}$   I don't see how this definition works. In my eyes, to even make the definition I gave above work with the theorem, we would have to "split" $A$ into powers of $1/2$ (whatever this means), i.e., look at $(A^{1/2}\mathbf{y})^TA^{1/2}\mathbf{y}$, unless I'm really missing something here.  This is just me doing self-reading, so any sources or anyone that can help me with a proof of this would be appreciated.
57420
Show that $(X_1,X_2)$ has a bivariate normal distribution with means $\mu_1, \mu_2$, variances $\sigma _1^2 $ and $\sigma _2^2$, and correlation coefficient $\rho $ if and only if every linear combination $t_1X_1+t_2X_2$ has a univariate normal distribution with mean $t_1\mu_1+t_2\mu_2$, and variances  $t_1^2\sigma _1^2+t_2^2\sigma _2^2+2t_1t_2\rho \sigma_{12}$ where $t_1, t_2$ are non zero constants  My work:  Joint mgf of $X_1 $and $X_2=e^{[t_1\mu_1+t_2\mu_2+\frac{1}{2}(t_1^2\sigma _1^2+t_2^2\sigma _2^2+2\rho \sigma_1 \sigma_1t_1t_2 )]} $  I am confused here to find the proper grouping to get the forward result.  For the backward result, I know I have to find the m.g.f. of $t_1X + t_2Y$, for any $t_1 $and $t_2$ real,and plug in the $E(t_1X + t_2Y)$ and $Var(t_1X + t_2Y)$.  I am struggling to get the correct backward result.  Second part:   Let $Y_i=X_i/\sigma_i, i=1,2$ Show that $Var(Y_1-Y_2)=2(1-\rho)$  I have found, $Var(Y_1)=\sigma_1, Var(Y_2)=\sigma_2$ since $X_1/\sigma_1=\sigma _1$  How can I prove that $Var(Y_1-Y_2)=2(1-\rho)$? Thank you so much for your help!
88502
Question  Suppose we have a multinomial distribution with $N$ possible outcomes, with probabilities $p_1,\ldots,p_N$. We sample this $n$ times, and denote the observed frequency of the $i$th outcome as $\xi_i$. In [1] the author claims that the distribution of the $\xi_i$ in the limit of large $n$ is:  $$p(\xi_1,\ldots,\xi_N)\propto\exp\left(-\frac{n}{2}\sum_{i=1}^N\frac{(\xi_i-p_i)^2}{p_i}\right).\;\;\;\;\;(1)$$  We can see immediately that this must be an approximation, as it assigns nonzero probabilities for $\xi_1+\cdots+\xi_N&gt;1$. However we can see that these have vanishing probability in the limit $n\rightarrow\infty$. My question is how do we derive (1) from the multinomial distribution, and show that they match in the $n\rightarrow\infty$ limit?  My thoughts  My first thought would be to appeal to the central limit theorem. The multinomial distribution has mean $\mu_i=p_i$ and covariance matrix $\Sigma_{ij}=\delta_{ij}p_i-p_ip_j$, so we would expect this in the large $n$ limit to be described by a multivariate Gaussian with mean $\mu$ and covariance $\frac{1}{n}\Sigma$. However, things are complicated by the fact that the multinomial covariance is singular (since $\xi_N$ is determined by the other $\xi_i$s), and so the multivariate Gaussian is not defined.  To address this, we may try and consider only the first $\xi_1,\ldots,\xi_{N-1}$, which have a non-singular covariance matrix and hence well-defined multivariate Gaussian distribution. Let's take the Binomial distribution $N=2$. The frequency $\xi_1$, this has mean $p_1$ and variance $p_1(1-p_1)$, so this would be described the the Gaussian: $$\propto\exp\left(-\frac{n}{2}\frac{(\xi_1-p_1)^2}{p_1(1-p_1)}\right).\;\;\;\;\;(2)$$ The expression (1) gives: $$\propto\exp\left(-\frac{n}{2}\left(\frac{(\xi_1-p_1)^2}{p_1}+\frac{(\xi_2-p_2)^2}{p_2}\right)\right).\;\;\;\;\;(3)$$ If we substitute $\xi_2\rightarrow 1-\xi_1$, $p_2\rightarrow 1-p_1$ into (3), we can verify that this gives the same answer as (2). I have verified that this also works for $N=4$.   I'm sure that if I just bashed out the algebra for general $N$ we would get agreement between the central limit theorem and (1) when we restrict the latter to $\xi_1+\cdots+\xi_N=1,p_1+\cdots+p_N=1$. However, how can we start with the multinomial distribution and derive (1) as a limit which is valid everywhere? One idea would be to say that (1) goes to zero as $n\rightarrow\infty$ when you are not on that plane, however I am a bit uncomfortable with this as it goes to zero everywhere except the mean as $n\rightarrow\infty$, so I don't know if that argument is good enough.  [1] Wootters, William K. "Statistical distance and Hilbert space." Physical Review D 23.2 (1981): 357.
38201
If $X \in \mathbb{C}^{m \times n}$ and $Y \in \mathbb{C}^{n \times m}$ ($m \geq n$), how to prove the following?  $$\sigma (XY) = \sigma(YX) \cup \underbrace{\left \{ 0, ..., 0 \right \}}_{m-n}$$  Here, $\sigma$ denotes the set of eigenvalues/spectrum.
9536
Given  $X \sim \mathcal{N}(\mu,\,\sigma^{2})$, what is the covariance between $X$ and $e^X$?
123566
Given that $a$~$\mathcal{N}(\mu, \Sigma)$ and $\mu$~$\mathcal{N}(0, \sigma^2 I)$, Is $p(a|\Sigma, \sigma)=\int p(a|\mu, \Sigma)p(\mu|\sigma)d\mu$ again a Normal distribution? And how to estimate the mean and variance respectively? Thanks a lot!
16510
Let $X_{1},\ldots,X_{n}\stackrel{\text{ i.i.d }}{\sim}N(\mu,\mu^{2})$.  $T=\left(\sum_{i=1}^{n}X_{i},\sum_{i=1}^{n}X_{i}^{2}\right)$ is a sufficient statistic for $\mu$. Also $T$ is minimal sufficient.  However, the question is: Prove that $T$ is not complete, hint: use $f(x,y)=x-2y^{2}/(n+1)$.  To prove that $T$ is not complete, I must find an example of a Borel function, $f$, for which $E_{\theta}(f(T))=0$ does not imply $f(T)=0$. However, I do not know how to interpret the 'hint', since I only can conclude that $f(T)=0$, whenever using $x=\sum_{i=1}^{n}X_{i}$ and $y=\sum_{i=1}^{n}X_{i}^{2}$ and calculating $E_{\theta}(f(T))=0$.
90165
I am trying to understand how 'completing the square' method is used to determine mean and co-variance of multivariate Gaussian distribution but I am not able to make much sense of the following step $$−((x−μ)^TΣ^-(x − μ))/2 = −(x^TΣ^−x + x^TΣ^-μ + const)/2$$ where x and μ are vectors and Σ is co-variance matrix.  What does this step indicate and how is it derived?   P.S. I am not well versed with linear algebra which might be a reason for not being able to comprehend this fully, if you can point me to certain properties from linear algebra that are being used here, that will be very helpful
50986
This question regards the problem of . Vectors and matrices will be denoted in bold.  Premises. Let $N,K$ be given integers, with $K \gg N &gt; 1$. The transpose of matrix $\mathbf{A}$ will be denoted with $\mathbf{A}^T$. Suppose the following statistical model holds $$ (*) \quad \mathbf{y} = \mathbf{Hx + n}, \quad \mathbf{n} \sim \mathcal{N}_{K}(\mathbf{0}, \mathbf{C})  $$ where $\mathbf{y} \in \mathbb{R}^{K \times 1}$ are the observables, $\mathbf{H} \in \mathbb{R}^{K \times N}$ is a known full-rank matrix, $\mathbf{x} \in \mathbb{R}^{N \times 1}$ is a deterministic vector of unknown parameters (which we want to estimate) and finally $\mathbf{n} \in \mathbb{R}^{K \times 1}$ is a disturbance vector (noise) with a known (positive definite) covariance matrix $\mathbf{C} \in \mathbb{R}^{K \times K}$. The  (ML) estimate of $\mathbf{x}$, denoted with $\hat{\mathbf{x}}_{ML}$, is given by $$ (1) \quad \hat{\mathbf{x}}_{ML} = (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{y} $$ and this is also the standard formula of Generalized Linear Least Squares (GLLS). Consider the standard formula of Ordinary Least Squares (OLS) for a linear model, i.e. $$ (2) \quad \hat{\mathbf{x}}_{OLS} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T \mathbf{y} $$ As a final note on notation, $\mathbf{I}_K$ is the $K \times K$ identity matrix and $\mathbf{O}$ is a matrix of all zeros (with appropriate dimensions). Now, for the problem at hand, assume that $\mathbf{C}^{-1} = \mathbf{I}_K + \mathbf{X}$, where $\mathbf{X} \in \mathbb{R}^{K \times K}$ is a symmetric, invertible matrix (and, for the formalism to make sense, $\mathbf{X}$ is such that $\mathbf{I}_K + \mathbf{X}$ is invertible and the inverse is positive definite). In this case, it can be proven (using matrix inversion lemma) that $$ (3) \quad (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T + \mathbf{Q} $$ where the expression for the matrix $\mathbf{Q} \in \mathbb{R}^{N \times K}$ can be found (I will omit it here). Now, finally,  Proposition 1. If $\mathbf{H}^T\mathbf{X} = \mathbf{O}_{N,K}$, then equation $(1)$ degenerates in equation $(2)$, i.e., there exists no difference between GLLS and OLS.  The proof is straigthforward and is valid even if $\mathbf{X}$ is singular. Now, my question is  Question: Can an equation similar to eq. $(3)$ (which "separates" an OLS-term from a second term) be written when $\mathbf{X}$ is a singular matrix?  I am not interested in a closed-form of $\mathbf{Q}$ when $\mathbf{X}$ is singular. But I do am interested in understanding the concept beyond that expression: what is the actual role of $\mathbf{Q}$? In which space does it operate?  I found this problem during a numerical implementation where both OLS and GLLS performed roughly the same (the actual model is $(*)$), and I cannot understand why OLS is not strictly sub-optimal. I found this slightly counter-intuitive, since you know a lot more in GLLS (you know $\mathbf{C}$ and make full use of it, why OLS does not), but this is somehow "useless" if some conditions are met. What are these conditions?  As a final note, I am rather new to the world of Least Squares, since I generally work within a ML-framework (or MMSE in other cases) and never studied the deep aspects of GLLS vs OLS, since, in my case, they are just intermediate steps during the derivation of MLE for a given problem.  If the question is, in your opinion, a bit too broad, or if there is something I am missing, could you please point me in the right direction by giving me references? Preferably well-known books written in standard notation.
72001
For the Quadratic Form $X^TAX; X\in\mathbb{R}^n, A\in\mathbb{R}^{n \times n}$ (which simplifies to $\Sigma_{i=0}^n\Sigma_{j=0}^nA_{ij}x_ix_j$), I tried to take the derivative wrt. X ($\Delta_X X^TAX$) and ended up with the following:  The $k^{th}$ element of the derivative represented as  $\Delta_{X_k}X^TAX=[\Sigma_{i=1}^n(A_{ik}x_k+A_{ki})x_i] + A_{kk}x_k(1-x_k)$  Does this result look right? Is there an alternative form?  I'm trying to get to the $\mu_0$ of Gaussian Discriminant Analysis by maximizing the log likelihood and I need to take the derivative of a Quadratic form. Either the result I mentioned above is wrong (shouldn't be because I went over my arithmetic several times) or the form I arrived at above is not the terribly useful to my problem (because I'm unable to proceed).  I can give more details about the problem or the steps I put down to arrive at the above result, but I didn't want to clutter to start off. Please let me know if more details are necessary.  Any link to related material is also much appreciated.
118473
To be specific, suppose we have $(n,1)$ random vector $x \sim N(\mu, \Sigma)$ where $\mu$ is $(n,1)$ and $\Sigma$ is $(n,n)$. Define: \begin{align*} Y &amp; = x'Ax + b'x + c \end{align*}  Then what is the following (for $t \in \mathbb{R}$)?  \begin{align*} E(e^{tY}) \end{align*}
61786
Same as . Let $X$ be a square integrable random variable on $(\Omega,\mathcal{F},P)$. Let $\mathcal{G}$ be a sub-$\sigma$-algebra of $\mathcal{F}$. Define the conditional variance of $X$ given $\mathcal{G}$ by $$\operatorname{var}(X\mid\mathcal{G})=\Bbb{E}[(X-\Bbb{E}[X\mid\mathcal{G}])^2\mid\mathcal{G}]$$ Prove the formula  $$\operatorname{var}(X)=\operatorname{var}(\Bbb{E}[X\mid\mathcal{G}])+\Bbb{E}[\operatorname{var}(X\mid\mathcal{G})]$$  1 I don't understand part of the presented solution.  2 Is my attempt also correct?  3 Alternate solutions?  1 I don't understand this part of the solution wherein  $$\operatorname{var}(X\mid \mathcal{G})=\Bbb{E}[(X-\Bbb{E}[X\mid\mathcal{G}])^2 \mid \mathcal{G}]$$  $ =\mathbb E[ X^2 \mid \mathcal G] - \mathbb E [X \mid \mathcal G]^2 $. I think it is making use of the fact that $\mathbb E [X \mid \mathcal G]$ is $\mathcal G$-measurable in saying:  $\mathbb E[XE[X\mid\mathcal G]\mid\mathcal G] = \mathbb E[X\mid\mathcal G] \ E[X\mid\mathcal G] = E[X\mid\mathcal G]^2$, but isn't $E[X\mid\mathcal G]$ supposed to be bounded? How to show this?  My attempt:  $E(X^2) &lt; \infty$  $\to E(X) &lt; \infty$  $\to E(X\mid\mathcal G) &lt; \infty$, but I don't think .  2 My own attempt: I tried evaluating $$E(\operatorname{var}(X\mid\mathcal{G}))=E(\Bbb{E}[(X-\Bbb{E}[X\mid\mathcal{G}])^2\mid\mathcal{G}])$$  Is it correct to say that $E(\Bbb{E}[(X-\Bbb{E}[X\mid\mathcal{G}])^2\mid\mathcal{G}]) = E([(X-\Bbb{E}[X\mid\mathcal{G}])^2])$?  After expanding the square I encounter $E[XE[X\mid\mathcal G]]$ and rewrite  $= E[E[XE[X\mid\mathcal G]]\mid\mathcal G]$  $ = E[E[X\mid\mathcal G] \ E[X\mid\mathcal G]]$ (*)  $ = E[E[X\mid\mathcal G]^2]$  (*) Assuming this is all correct, I am still unsure about this part. Is it really bounded?  3 Is there a way to go about this without saying $E[X\mid\mathcal G]$ is bounded?  4 Is this wrong? It's from Stochastic Calculus class notes, though the problem is for an Advanced Probability class (and the notes don't say anything about bounded weirdly).
36151
It's a standard result that given $X_1,\cdots ,X_n $ random sample from $N(\mu,\sigma^2)$, the random variable $$\frac{(n-1)S^2}{\sigma^2}$$ has a chi-square distribution with $(n-1)$ degrees of freedom, where $$S^2=\frac{1}{n-1}\sum^{n}_{i=1}(X_i-\bar{X})^2.$$ I would like help in proving the above result. Thanks.
33225
Let $X_1$ and $X_2$ are independent $N(0, \sigma^2)$ which means (mean = 0, variance = $\sigma^2$) random variables. What is the distribution of $X_1^2 + X_2^2$?  My approach is that  $X_1\sim N(0, \sigma^2)$ and $X_2\sim N(0, \sigma^2)$  Then $X_1^2$ and $X_2^2$ have chi-squared distribution with 1 degree of freedom. (I am not sure the degree of freedom and not sure how to show it as well(please help on this))  Then I found the moment-generating function for $X_1^2$ and $X_2^2$;$$m_{X_1^2} = (1-2t)^{-1/2}$$ and $$m_{X_2^2} = (1-2t)^{-1/2}$$  So the moment generating function for $X_1^2 + X_2^2$ is $$m_{X_1^2}(t) m_{X_2^2}(t) = (1-2t)^{-2/2}$$  So $X_1^2 + X_2^2$ has a chi-squared distribution with 2 degrees of freedom. Is this correct?
4619
Consider the linear regression model:   $$Y_i=r(x_i)+\varepsilon_i\equiv\sum_{j = 1}^p x_{ij} \beta _j + \varepsilon _i,\quad i=1,\ldots,n.$$  where $x_1,\ldots,x_n\in \mathbb{R}^p$ are fixed, $E(\varepsilon_i)=0$,  $\operatorname{Var}(\varepsilon_i)=\sigma^2$. Denote $Y=(Y_1,\ldots,Y_n)^T$, $\beta=(\beta_1,\ldots,\beta_p)^T$, $X=(x_{ij})_{n\times p}$. As is known, $\hat \beta  = \arg\min\limits_{\beta  \in \mathbb{R}^p} (Y - X\beta)^T (Y - X\beta)=(X^TX)^{-1}X^TY$ if the matrix $X^TX$ is invertible, and so an estimator of $r(x)$ at $x=(x_1,\ldots,x_p)\in\mathbb{R}^p$ is given by $\hat r_n(x)=x^T\hat\beta$.  Let $$\hat{\sigma}^2 = \frac{1}{\sigma^2}\sum_{i = 1}^n (Y_i - {\hat r}_n (x_i))^2.$$ I am stucking the problem: $\hat\sigma^2$ has $\chi^2_{n-p}$ distribution. How to prove the statement?
70859
We have a multivariate normal vector $\mathbf{X} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$, where $\mathbf{X} =  \left[ \begin{matrix} X_1 \\ X_2 \\ \dot\\\\ \dot\\\\ X_n \end{matrix} \right]$,    $\boldsymbol{\mu} =  \left[ \begin{matrix} \mu_1 \\ \mu_2 \\ \dot\\\\ \dot\\\\ \mu_n \end{matrix} \right]$, and  $\boldsymbol{\Sigma} =  \begin{bmatrix} \sigma_1^2 &amp; \sigma_1\sigma_2 &amp; \dots &amp; \sigma_1\sigma_n   \\ \sigma_2\sigma_1 &amp; \sigma_2^2 &amp; \dots &amp; \sigma_2\sigma_n  \\ \dot\\\\ \dot\\\\ \sigma_n\sigma_1 &amp; \sigma_n\sigma_2  &amp; \dots &amp; \sigma_n^2 \end{bmatrix} $ is the covariance matrix.  Y is a random scalar of the form $Y = \beta_0 + \boldsymbol{\beta^TX} + \epsilon$, where $\epsilon \sim \mathcal{N}(0, \sigma_\epsilon^2)$  How can we go about deriving the joint probability density function(pdf) of the random vector $(Y, \boldsymbol{X^T})$?  The marginal multivariate pdf of $\boldsymbol{X}$ is given by  $ pdf(\boldsymbol{x}) = \frac{1}{\sqrt{(2\pi})^n. \begin{vmatrix} \boldsymbol{\Sigma} \end{vmatrix}}.exp\left[ \frac{-1}{2}\boldsymbol{(x-\mu)}^T \boldsymbol{\Sigma^{-1}}\boldsymbol{(x-\mu)} \right]$  The conditional probability density function of $Y \,\mid \, \boldsymbol{X=x}$ is given by  $pdf(y\mid \boldsymbol{X=x}) = \frac{1}{2\pi\sigma_\epsilon}.exp\left[  \frac{-1}{2\sigma_\epsilon^2}(y - \beta_0 - \boldsymbol{\beta^Tx})^2  \right]   $  To compute the joint probability density function, we can use the fact that   $pdf(y, \boldsymbol{x})=  pdf(y\mid \boldsymbol{X=x}).pdf(\boldsymbol{x})$  I am unable to do the witty manipulation required in the algebra to arrive at the result the joint distribution will have a probability density function of a normal distribution. All questions previously asked seem to address the problem of deriving the conditional density function given the joint probability density function, for example,   Edit:  Multiplying the pdfs, we get $pdf(y, \boldsymbol{x}) = \frac{1}{\sqrt{(2\pi})^n. \begin{vmatrix} \boldsymbol{\Sigma} \end{vmatrix}}. \frac{1}{2\pi\sigma_\epsilon}.exp\left[ \frac{-1}{2}\boldsymbol{(x-\mu)}^T \boldsymbol{\Sigma^{-1}}\boldsymbol{(x-\mu)} + \frac{-1}{2\sigma_\epsilon^2}(y - \beta_0 - \boldsymbol{\beta^Tx})^2 \right]$   Like the comment below mentions, most manipulations dealing with normal pdfs involve completing squares; but since this specific algebraic expression contains $\boldsymbol{x}$ in both terms within the exp term, I'm unsure on how to go about completing the squares.