Probability density function for white Gaussian noise in many signal processing text books and lectures we find that if we assume that the noise is white Gaussian then the probability density function itself takes the Gaussian form (see here for example) when trying to estimate parameters through the maximum-likelihood estimation method.
I do not understand this leap, why just because the noise is Gaussian the parameters themselves are Gaussian distributed parameters? I do not see how the white Gaussian noise fits into the probability density function at all! It seems we are always just guessing that the probability density function is normally distributed. Am I wrong? Or can anyone help me understand this or point me in a direction that does? Thank you very much.
 A: As specified in the comments:

what I do not understand is how a linear model with Gaussian noise produces Gaussian data

This is because the family of normal distributions is closed under linear transformations: simply put, once you've got a normally distributed random variable, you can't make it not normal by addition or multiplication with scalars. Let $X \sim \mathcal{N}(0, 1)$. Then for any constants $a,b$:
$$ a X + b = Y \sim \mathcal{N}(b, a^2)$$
In the stochastic process setting, this $Y$ is the data, $X$ is noise, and $b$ is defined by the fixed effects (what's sometimes called the DC offset in DSP, or intercept if this was a basic regression model). Apply the above equation and you'll get the needed distribution of $Y$.
A: There are many questions and answers on the sister site dsp.SE dealing with white noise and white Gaussian noise etc. The following is a somewhat adapted version of part of one answer on dsp.SE that I wrote.
(Continuous-time) white noise is a mythical process that is unobservable in all its glory in nature (probably just as well since it is infinitely powerful and would lead to an immediate solution to the energy crisis).  We poor mortals can only observe white noise through some kind of device that necessarily limits what we can observe -- kind of like watching a solar eclipse through special glasses--  and thus what we observe is a pale imitation of the real thing. Well, it has been observed that if an observation device is modeled as a linear filter with transfer function $H(f)$, then (with an open circuit at the filter input), the filter output is a wide-sense-stationary (also called a weakly stationary) random process with power spectral density $K|H(f)|^2$.  This is consistent with an assumption that the input to the filter is a white noise process with autocorrelation function $K\delta(t)$ (where $\delta(t)$ denotes a Dirac delta or impulse)  and power spectral density $S(f) = K, -\infty < f < \infty$ if we simply plug in $K$ for the input power spectral density in the power spectral density equation
$$S_{\text{output}}(f) = |H(f)|^2 S_{\text{input}}(f).$$
Never mind that many mathematicians will cringe at the cavalier treatment where we are ignoring that the above formula
implicitly assumes that the input process is a finite power process (which white noise is definitely not); but the final result is correct even if the process of arriving at the result is not. Notice however, that all of the above tells us very little about the probability distributions of the random variables constituting the output process (or the input process for that matter).
Now, (continuous-time) white Gaussian noise is also a mythical process with the extra property that the output process of a linear filter with transfer function not only has power spectral density $K|H(f)|^2$ but also that the the output process is a Gaussian process, which means, among lots of other things, that all the random variables comprising the process are Gaussian random variables and that any finite subset of the variables has a joint Gaussian distribution. Of course, wide-sense-stationary Gaussian processes are also strictly stationary.  Now, the standard theory of Gaussian processes in linear systems says that if the input ta linear system is a Gaussian process, the output is also a Gaussian process, and hence we bestow the adjective Gaussian on white noise processes that result in Gaussian processes when they pass through linear filters, but it is not appropriate to reverse-engineer the output process being a Gaussian process to say that all the random variables in the input white noise process also are Gaussian random variables unless we want to stretch the definition of a Gaussian random variable to include the case of the variance being infinite. Skeptical beginners should try writing out the pdf of a hypothetical $\mathcal N(0,\infty)$ random variable......
Finally, turning to discrete-time random processes which I think is what the OP really wants to know about, remember that one cannot sample the mythical beast called continuous-time white noise -- it does not exist in nature -- and the sampler is necessarily a device that observes the random process for a very short but nonzero time $\varepsilon$, and thus the sample $X[n]$ is actually something proportional to $\int_{nT-\varepsilon/2}^{nT+\varepsilon/2}X_t \mathrm dt$ which has variance $\sigma^2\varepsilon$ if $\{X_t\}$ is a white noise process.  So,

A discrete-time white noise process is a collection of zero-mean independent identically distributed random variables $X[n]$.
A discrete-time white Gaussian noise process is a collection of zero-mean independent identically distributed Gaussian random variables $X[n]$.

Yes, many DSP and statistics texts (as well as Wikipedia's definition of a discrete-time white noise process) and many people with much higher reputation than me on dsp.SE and stats.SE say that uncorrelatedness suffices for defining a white noise process, and in the case of white Gaussian noise it does because Gaussianity brings in the jointly Gaussian property: a discrete-time Gaussian random process is defined as a sequence of random  variables $\{X[n]\colon n \in \mathbb Z\}$ such that any set of $M\geq 1$ random variables $X[n_1], X[n_2], \ldots, X[n_M]$ enjoys a jointly Gaussian distribution, and so for white Gaussian noise, uncorrelatedness implies independence. However, for arbitrary white noises, it is best to insist on independence and not on just zero correlation. For the edification of all these important people who insist that uncorrelatedness is adequate, I present a discrete-time process in which every random variable is a Gaussian random variable, any two random variables are uncorrelated but are not necessarily independent, and not all sets of variables in the process enjoy a jointly Gaussian distribution. In short, the process defined bellow is not a discrete-time white Gaussian noise process as per anybody's standard definition. And why should all this matter in the least? Well, in typical applications we apply various mathematical operations on processes, and if $X[0]$ and $X[1]$ are uncorrelated Gaussian random variables and we cannot rely on $X[0]+X[1]$ also being a Gaussian random variable, things have come to pretty pass, and it is not a world I want to live in.

Example: Let $X$ be a $N(0,1)$ random variable and $B$ a discrete random variable that takes on values $+1$ and $-1$ with equal probability $\frac 12$ and independent of $X$. Set $Y = BX$ and note that $E[Y]=E[BX]=E[B]E[X]=0$. Furthermore,
$E[XY] = E[X^2B] = E[X^2]E[B] = 0$, and so $X$ and $Y$ are uncorrelated random variables. But what is the distribution of $Y$? Well,
\begin{align}
P(Y \leq a) &= P(Y\leq a \mid B=+1)P(B=+1) + P(Y\leq a \mid B=-1)P(B=-1)\\
&= \frac 12 P(BX\leq a \mid B=+1) + \frac 12 P(BX\leq a \mid B=-1)\\
&= \frac 12 P(X\leq a) + \frac 12 P(X\geq -a)\\
&= \frac 12 \Phi(a) +  \frac 12 \Phi(a)\\
&= \Phi(a),
\end{align}
that is, $Y$ is also a $N(0,1)$ random variable!!  But $X$ and $Y$ are not jointly Gaussian random variables. Note that conditioned on the value of $X$ being $\alpha$, $Y$ is a discrete random variable that takes on values $\pm\alpha$ with equal probability: with joint Gaussianity, $Y$ would have been a Gaussian random variable. 
With this as background, let $\{X[2n]\colon n \in \mathbb Z\}$ be a set of independent identically distributed zero-mean Gaussian random variables, that is, a standard discrete-time white Gaussian noise process on the even integers. Let $\{B[n]\colon n \in \mathbb Z\}$ be an independent process where the $B[n]$'s are independent discrete random variables  that take on values $+1$ and $-1$ with equal probability $\frac 12$. Set
$X[2n+1] = X[2n]B[n]$ and note that each pair $(X[2n],X[2n+1])$ is a pair of uncorrelated zero-mean Gaussian random variables that are not jointly Gaussian.  Now let's look at the random process
$\{X[m]\colon m \in \mathbb Z\}$ in which all the random variables are zero-mean Gaussian with the same variance. Any pair of random variables is uncorrelated: $X[2n]$ and $X[2n+1]$ by construction and all the more distant pairs because of independence. But, not all pairs of random variables have a jointly Gaussian distribution and so this is not a white Gaussian  noise process in the usual sense of the term; ymmv.

