Maximum likelihood estimation (MLE) approach

Since \({\bf y}\sim \mathcal{N}_n({\bf X}{\boldsymbol \beta}, \sigma^2{\bf I}_n)\), the likelihood function is given as \[\begin{eqnarray*} L({\boldsymbol \beta},\sigma^2)&=&f({\bf y}|{\boldsymbol \beta},\sigma^2)\\ &= &(2\pi)^{-\frac{n}{2}}\lvert {\boldsymbol \Sigma}\lvert^{-\frac{1}{2}} \exp\left\{-\frac{1}{2}({\bf y}-{\bf X}{\boldsymbol \beta})^{{\bf T}}{\boldsymbol \Sigma}^{-1}({\bf y}-{\bf X}{\boldsymbol \beta})\right\} , \end{eqnarray*}\]

where \({\boldsymbol \Sigma}=\sigma^2 {\bf I}_n\). Then the log likelihood can be written as \[\begin{eqnarray*} l({\boldsymbol \beta},\sigma^2)&=&\log L({\boldsymbol \beta},\sigma^2)\\ &=&-\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(\sigma^2) - \frac{1}{2\sigma^2}({\bf y}-{\bf X}{\boldsymbol \beta})^{{\bf T}}({\bf y}-{\bf X}{\boldsymbol \beta}). \end{eqnarray*}\] Note that \(l({\boldsymbol \beta},\sigma^2)\) is a concave function in \(({\boldsymbol \beta},\sigma^2)\). Hence, the maximum likelihood estimator (MLE) can be obtained by solving the following equations: \[\begin{eqnarray*} \frac{\partial l({\boldsymbol \beta},\sigma^2)}{\partial {\boldsymbol \beta}} &=&- \frac{1}{2\sigma^2}(-2{\bf X}^{{\bf T}} {\bf y}+ 2{\bf X}^{{\bf T}} {\bf X}{\boldsymbol \beta})=0; \\ \frac{\partial l({\boldsymbol \beta},\sigma^2)}{\partial \sigma^2} &=&-\frac{n}{2}\frac{1}{\sigma^2} + \frac{1}{2}\frac{1}{(\sigma^2)^2} \mathbf{({\bf y}-{\bf X}{\boldsymbol \beta})^T({\bf y}-{\bf X}{\boldsymbol \beta})}=0. \end{eqnarray*}\] Therefore, the MLEs of \({\boldsymbol \beta}\) and \(\sigma^2\) are given as \[\begin{eqnarray*} \hat{\beta} &=& ({\bf X}^{{\bf T}} {\bf X})^{-1}{\bf X}^{{\bf T}} {\bf y};\\ \hat{\sigma}^2 &=& \frac{({\bf y}-{\bf X}\hat{{\boldsymbol \beta}})^{{\bf T}}({\bf y}-{\bf X}\hat{{\boldsymbol \beta}})}{n}=\frac{ \|{\bf y}-\hat{{\bf y}}\|^2}{n}, \end{eqnarray*}\] where \(\hat{{\bf y}}={\bf X}\hat{{\boldsymbol \beta}}\).

Distribution of \(\hat{\boldsymbol \beta}\) and \(\hat\sigma^2\)

Note that if \({\bf z}\sim \mathcal{N}_k(\mu,\Sigma)\), then \(A{\bf z}\sim \mathcal{N}_m(A\mu,A\Sigma A^{{\bf T}})\), where \(A\) is an \(m\times k\) matrix. Since \({\bf y}\sim \mathcal{N}_n({\bf X}{\boldsymbol \beta}, \sigma^2{\bf I}_n)\) and \(\hat{{\boldsymbol \beta}} =({\bf X}^{{\bf T}} {\bf X})^{-1}{\bf X}^{{\bf T}}{\bf y}\), we have \[\begin{eqnarray*} \hat{{\boldsymbol \beta}} &\sim& \mathcal{N}_p\left( ({\bf X}^{{\bf T}} {\bf X})^{-1}{\bf X}^{{\bf T}}{\bf X}{\boldsymbol \beta}, \sigma^2 ({\bf X}^{{\bf T}} {\bf X})^{-1}{\bf X}^{{\bf T}}{\bf X}({\bf X}^{{\bf T}} {\bf X})^{-1}\right)\\ &=&\mathcal{N}_p\left({\boldsymbol \beta}, \sigma^2 ({\bf X}^{{\bf T}} {\bf X})^{-1}\right). \end{eqnarray*}\]

Note that \({\bf y}-\hat{{\bf y}}=({\bf I}_n-{\bf X}({\bf X}^{{\bf T}}{\bf X})^{-1}{\bf X}^{{\bf T}}){\bf y}\), where \({\bf I}_n-{\bf X}({\bf X}^{{\bf T}}{\bf X})^{-1}{\bf X}^{{\bf T}}\) is an idempotent matrix with rank \((n-p)\).

\(\color{red}{\text{Can you prove that ${\bf I}_n-{\bf X}({\bf X}^{{\bf T}}{\bf X})^{-1}{\bf X}^{{\bf T}}$ s an idempotent matrix of rank $(n-p)$ ?}}\)

\(\color{red}{\text{How to prove $r({\bf I}_n-{\bf H})=tr({\bf I}_n-{\bf H})$? }}\)

The eigenvalue of idempotent matrix is either 1 or 0, hence the rank of it is the sum of eigenvalues, which equals the trace of matrix;

\(\color{red}{\text{How to prove trace of matrix is the sum of eigenvalues?}}\)

It requires characteristic polynomial.

Another way to prove this is in this link

From Lemma 1, we have that \[\begin{eqnarray} \tag{1} n\frac{\hat{\sigma}^2}{\sigma^2}\sim \chi^2(n-p), \end{eqnarray}\] where \(\chi^2(n-p)\) denotes the chi-squared distribution with degrees of freedom \(n-p\).

\({\color{red}{\text{Can you prove Eq.}}}\)(1)?

\(\sigma^2\) is known

Suppose \(\sigma^2\) is known. We define the prior distribution of \({\boldsymbol \beta}\) by \({\boldsymbol \beta}\sim \mathcal{N}_p({\boldsymbol 0},\sigma^2\nu{\bf I}_p)\). Then the posterior density of \({\boldsymbol \beta}\) can be obtained by \[\begin{eqnarray*} \pi ({\boldsymbol \beta}|{\bf y}) &\propto& f({\bf y}|{\boldsymbol \beta})\pi ({\boldsymbol \beta}) \\ &\propto& \exp\left(-\frac{1}{2\sigma^2}\| {\bf y}-{\bf X}{\boldsymbol \beta}\|^2\right)\times \exp\left(-\frac{1}{2\sigma^2\nu}\|{\boldsymbol \beta}|^2\right)\\ &=&\exp\left[-\frac{1}{2\sigma^2}\left\{ ({\bf y}-{\bf X}{\boldsymbol \beta})^{{\bf T}}({\bf y}-{\bf X}{\boldsymbol \beta}) + \frac{1}{\nu}{\boldsymbol \beta}^{{\bf T}} {\boldsymbol \beta}\right\}\right]\\ &\propto&\exp\left\{-\frac{1}{2\sigma^2}\left( - 2 {\boldsymbol \beta}^{{\bf T}} {\bf X}^{{\bf T}} {\bf y}+ {\boldsymbol \beta}^{{\bf T}} {\bf X}^{{\bf T}} {\bf X}{\boldsymbol \beta}+ \frac{1}{\nu}{\boldsymbol \beta}^{\bf T}{\boldsymbol \beta} \right)\right\}\\ &=&\exp\left\{-\frac{1}{2\sigma^2}\left( - 2 {\boldsymbol \beta}^{{\bf T}} \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right) \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right)^{-1} {\bf X}^{{\bf T}} {\bf y}+ {\boldsymbol \beta}^{{\bf T}} \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right){\boldsymbol \beta}\right)\right\}\\ &\propto&\exp\left\{-\frac{1}{2\sigma^2}({\boldsymbol \beta}-\tilde{{\boldsymbol \beta}})^{{\bf T}} \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right)({\boldsymbol \beta}-\tilde{{\boldsymbol \beta}}) \right\} \end{eqnarray*}\] where \(\tilde{{\boldsymbol \beta}}=\left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}\right)^{-1} {\bf X}^{{\bf T}} {\bf y}\).

This implies that \[\begin{eqnarray} {\boldsymbol \beta}|{\bf y}\sim \mathcal{N}_p\left(\left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right)^{-1} {\bf X}^{{\bf T}} {\bf y},\sigma^{2} \left({\bf X}^{{\bf T}} {\bf X} +\frac{1}{\nu}I_p\right)^{-1}\right).\tag{2} \end{eqnarray}\] The Bayesian estimate \(\hat{\boldsymbol \beta}_{Bayes} = \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}I_p\right)^{-1} {\bf X}^{{\bf T}} {\bf y}\). It is worth noting that if \(\nu\to \infty\), the posterior distribution converges to the distribution of \(\hat{{\boldsymbol \beta}}_{MLE}\sim \mathcal{N}_p\left({\boldsymbol \beta}, \sigma^2 ({\bf X}^{{\bf T}} {\bf X})^{-1}\right)\).

\(\sigma^2\) is unknown

In general, \(\sigma^2\) is unknown. Then, we need to assign a reasonable prior distribution for \(\sigma^2\). We consider the inverse gamma distribution, which is called a , as follows: \(\sigma^2\sim \mathcal{IG}(a_0,b_0)\) with the density function \[\begin{eqnarray} \pi(\sigma^2)=\frac{b_0^{a_0}}{\Gamma(a_0)}(\sigma^2)^{-a_0-1}\exp\left(-\frac{b_0}{\sigma^2}\right), \end{eqnarray}\] where \(a_0>0\) and \(b_0>0\). In addition, we need to introduce prior for \({\boldsymbol \beta}|\sigma^2\):

\[{\boldsymbol \beta}|\sigma^2 \sim \mathcal{N}_p({\boldsymbol 0},\sigma^2\nu{\bf I}_p).\]

\(\color{red}{\text{Today, we derive the joint posterior distribution} \pi({\boldsymbol \beta},\sigma^2|{\bf y})\propto f({\bf y}|{\boldsymbol \beta},\sigma)\pi({\boldsymbol \beta}|\sigma^2)\pi(\sigma^2).}\)

Show

• 1. \(\pi(\sigma^2|{\boldsymbol \beta},{\bf y}) = \mathcal{IG}(a^*,b^*)\)
• 2. \(\pi({\boldsymbol \beta}|{\bf y}) \sim\) t-distribution with \(t^*\)

Note: the density of a multiple t-distribution with \({\boldsymbol \Sigma},{\boldsymbol \mu}\) and \(\nu\) is (3) \[\begin{eqnarray*} \frac{\Gamma[(\nu+p)/2]}{\Gamma(\nu/2)\nu^{p/2}\pi^{p/2}|{\boldsymbol \Sigma}|^{1/2}}\left[1+\frac{1}{\nu}({\bf X}- {\boldsymbol \mu})^{{\bf T}} \Sigma^{-1}({\bf X}-{\boldsymbol \mu})\right]^{-(\nu+p)/2}. \end{eqnarray*}\]

\(\sigma^2\) is known

\(\sigma^2\) is unknown

Suppose \(\sigma^2\) is unknown, we cannot use \(\pi({\boldsymbol \beta}|{\bf y},\sigma^2)\) directly.

But we know \(\pi({\boldsymbol \beta},\sigma^2|{\bf y}) = \pi({\boldsymbol \beta}|{\bf y},\sigma^2)\pi(\sigma^2|{\bf y})\) and we already know

\({\boldsymbol \beta}|{\bf y},\sigma^2 \sim \mathcal{N}_p\left(({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p)^{-1}{\bf X}^{{\bf T}}{\bf y}, ({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p)^{-1}\sigma^2\right)\)

So \[\begin{eqnarray*} &&\pi(\sigma^2|{\bf y}) = \int \pi(\sigma^2, {\boldsymbol \beta}| {\bf y}) d{\boldsymbol \beta} \propto \int f({\bf y}| \sigma^2, {\boldsymbol \beta})\pi({\boldsymbol \beta}|\sigma^2)\pi(\sigma^2)d{\boldsymbol \beta} =\int f({\bf y}|\sigma^2,{\boldsymbol \beta})\pi({\boldsymbol \beta}|\sigma^2)d{\boldsymbol \beta}\times \pi(\sigma^2)\\ &\propto& \int (\sigma^2)^{-\frac{n}{2}} \exp\left[-\frac{1}{2\sigma^2}({\bf y}- {\bf X}{\boldsymbol \beta})^{{\bf T}} ({\bf y}- {\bf X}{\boldsymbol \beta})\right](\sigma^2)^{-\frac{p}{2}}\exp(-\frac{1}{2\sigma^2\nu}{\boldsymbol \beta}^{{\bf T}}{\boldsymbol \beta}) d{\boldsymbol \beta}\times \pi(\sigma^2)\\ &\propto& \int \exp\left[-\frac{1}{2\sigma^2}({\bf y}^{{\bf T}}{\bf y}-2 {\boldsymbol \beta}^{{\bf T}}{\bf X}^{{\bf T}}{\bf y}+ {\boldsymbol \beta}^{{\bf T}}{\bf X}^{{\bf T}}{\bf X}{\boldsymbol \beta}+ \frac{1}{\nu}{\boldsymbol \beta}^{{\bf T}}{\boldsymbol \beta})\right]d{\boldsymbol \beta}\left((\sigma^2 )^{-\frac{1}{2}(n + p)}\pi(\sigma^2)\right)\\ &=& \int \exp \left[-\frac{1}{2\sigma^2}\left({\boldsymbol \beta}^{{\bf T}}({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p) {\boldsymbol \beta}-2{\boldsymbol \beta}^{{\bf T}}({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p)({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p)^{-1}{\bf X}^{{\bf T}} {\bf y}+ {\bf y}^{{\bf T}}{\bf y}\right)\right]d{\boldsymbol \beta}\\ &&\left((\sigma^2)^{-\frac{1}{2}(n + p)}\pi(\sigma^2)\right)\\ &=&\int\exp \left[-\frac{1}{2\sigma^2}({\boldsymbol \beta}- \tilde {\boldsymbol \beta})^{{\bf T}}{\bf A}({\boldsymbol \beta}- \tilde {\boldsymbol \beta}) + {\bf y}^{{\bf T}}{\bf y}-\tilde {\boldsymbol \beta}^{{\bf T}}{\bf A}\tilde {\boldsymbol \beta}\right]d{\boldsymbol \beta} \left((\sigma^2)^{-\frac{1}{2}(n + p)}\pi(\sigma^2)\right)\\ &=&\int\exp \left[-\frac{1}{2\sigma^2}({\boldsymbol \beta}- \tilde {\boldsymbol \beta})^{{\bf T}}{\bf A}({\boldsymbol \beta}- \tilde {\boldsymbol \beta}) \right]d{\boldsymbol \beta}\times \exp\left[-\frac{1}{2\sigma^2}\left({\bf y}^{{\bf T}}{\bf y}-{\bf y}^{{\bf T}} {\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}{\bf y}\right)\right]\left((\sigma^2)^{-\frac{1}{2}(n + p)}\pi(\sigma^2)\right)\\ &=& (2\pi)^{\frac{p}{2}} \lvert\sigma^2{\bf A}^{-1}\rvert^{\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^2} \left({\bf y}^{{\bf T}}(I_n - {\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}\right)\right]\left((\sigma^2)^{-\frac{1}{2}(n + p)} \pi(\sigma^2)\right)\\ &\propto&(\sigma^2)^{\frac{p}{2} - \frac{1}{2}(n+p)}(\sigma^2)^{-a_0 -1} \exp \left[-\frac{1}{\sigma^2}\left(b_0 + \frac{1}{2} {\bf y}^{{\bf T}}(I_n - {\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}\right) \right] \\ &=& (\sigma^2)^{-(\frac{1}{2}n + a_0)-1}\exp \left[-\frac{1}{\sigma^2} \left(b_0 + \frac{1}{2} {\bf y}^{{\bf T}}(I_n - {\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}\right)\right] \\ \\ &=& (\sigma^2)^{-a^*-1}\exp(-\frac{b^*}{\sigma^2}) \end{eqnarray*}\] where \({\bf A}= ({\bf X}^{{\bf T}}{\bf X}+\frac{1}{\nu}I_p),\quad \tilde {\boldsymbol \beta}= {\bf A}^{-1}{\bf x}^{{\bf T}}{\bf y}, \quad a^* =\frac{1}{2}n+ a_0,\quad b^*=b_0 + \frac{1}{2} {\bf y}^{{\bf T}}(I_n - {\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}\).

This is proportional to the pdf of \(\mathcal{IG}(a^*, b^*)\). Hence we have \[\begin{eqnarray*} E({\boldsymbol \beta}|y) &=& \int{\boldsymbol \beta}\pi({\boldsymbol \beta}|{\bf y})d{\boldsymbol \beta}= \int{\boldsymbol \beta}[ \int\pi({\boldsymbol \beta},\sigma^2|{\bf y})d\sigma^2 ]d{\boldsymbol \beta}\\ &=& \int\int{\boldsymbol \beta}[\pi({\boldsymbol \beta}|{\bf y},\sigma^2)\pi(\sigma^2|{\bf y})d\sigma^2]d{\boldsymbol \beta}, \end{eqnarray*}\] Similarly, \[\begin{eqnarray*} Var({\boldsymbol \beta}|y) &=& \int({\boldsymbol \beta}-E({\boldsymbol \beta}|y))\pi({\boldsymbol \beta}|{\bf y})d{\boldsymbol \beta}= \int({\boldsymbol \beta}-E({\boldsymbol \beta}|y))[ \int\pi({\boldsymbol \beta},\sigma^2|{\bf y})d\sigma^2 ]d{\boldsymbol \beta}\\ &=&\int \int({\boldsymbol \beta}-E({\boldsymbol \beta}|y))[\pi({\boldsymbol \beta}|{\bf y},\sigma^2)\pi(\sigma^2|{\bf y})d\sigma^2]d{\boldsymbol \beta}, \end{eqnarray*}\]

Multivariable Normal Conditional Distribution

• The conditional posterior of \(\beta_j\) is obtained as follows: \[\begin{eqnarray*} \pi(\beta_j|{\boldsymbol \beta}_{-j},\sigma^2,{\bf y}) &\propto& \pi(\beta_j,{\boldsymbol \beta}_{-j}|\sigma^2,{\bf y})\\ &=& \pi({\boldsymbol \beta}|\sigma^2,{\bf y})\\ &\propto&f({\boldsymbol \beta},{\bf y}|\sigma^2)\\ &=& f({\bf y}|{\boldsymbol \beta},\sigma^2)\pi({\boldsymbol \beta}|\sigma^2)\\ &\propto&\exp \left(-\frac{1}{2\sigma^2}\parallel {\bf y}-{\bf X}{\boldsymbol \beta}\parallel^2 \right) \exp\left(-\frac{1}{2\sigma^2\nu}\parallel{\boldsymbol \beta}\parallel^2\right)\\ &\propto&\exp\left[-\frac{1}{2\sigma^2}\left(({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j}-x_j\beta_j)^{{\bf T}} ({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j}-x_j\beta_j)+\frac{1}{\nu}\beta_j^{{\bf T}}\beta_j \right) \right]\\ &\propto&\exp \left[-\frac{1}{2\sigma^2} (-{\bf y}^{{\bf T}}x_j \beta_j + {\boldsymbol \beta}_{-j}^{{\bf T}}{\bf X}_{-j}^{{\bf T}}x_j \beta_j- \beta_j^{{\bf T}}x_j^{{\bf T}}{\bf y}+ \beta_j^{{\bf T}}x_j^{{\bf T}} {\bf X}_{-j}{\boldsymbol \beta}_{-j}+ \beta_j^{{\bf T}}x_j^{{\bf T}}x_j\beta_j + \frac{1}{\nu}\beta_j^{{\bf T}}\beta_j)\right]\\ &\propto&\exp \left[-\frac{1}{2\sigma^2} (-2\beta_j^{{\bf T}}x_j^{{\bf T}}{\bf y}+ 2\beta_j^{{\bf T}}x_j^{{\bf T}} {\bf X}_{-j}{\boldsymbol \beta}_{-j}+\beta_j^{{\bf T}}x_j^{{\bf T}}x_j \beta_j + \frac{1}{\nu}\beta_j^{{\bf T}}\beta_j)\right]\\ &=&\exp \left[-\frac{1}{2\sigma^2} \left(\beta_j^{{\bf T}}(x_j^{{\bf T}}x_j +\frac{1}{\nu})\beta_j -2\beta_j^{{\bf T}}x_j^{{\bf T}}({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j})\right)\right]\\ &=&\exp \left[-\frac{1}{2\sigma^2} (x_j^{{\bf T}}x_j +\frac{1}{\nu})\left(\beta_j^{{\bf T}}\beta_j -2\beta_j^{{\bf T}}(x_j^{{\bf T}}x_j +\frac{1}{\nu})^{-1}x_j^{{\bf T}}({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j})\right)\right]\\ &=&\exp \left[-\frac{1}{2\sigma^2}(x_j^{{\bf T}}x_j +\frac{1}{\nu})(\beta_j-\tilde \beta_j)^2 \right] \end{eqnarray*}\] where \(\tilde \beta_j=(x_j^{{\bf T}}x_j +\frac{1}{\nu})^{-1}x_j^{{\bf T}}({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j})\), \({\bf X}_{-j}\) is a submatrix of \({\bf X}\) without the \(j^{th}\) column, and \({\boldsymbol \beta}_{-j}\) is a subvector of \({\boldsymbol \beta}\) without the \(j^{th}\) element. Hence

\[ \beta_j|{\boldsymbol \beta}_{-j},\sigma^2,{\bf y}\sim \mathcal{N}\left((x_j^{{\bf T}}x_j +\frac{1}{\nu})^{-1} x_j^{{\bf T}}({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j}),\sigma^2\left[x_j^{{\bf T}}x_j +\frac{1}{\nu}\right]^{-1}\right) \]

Setup

• Consider

\[ f({\bf y}|{\boldsymbol \beta},\sigma^2) = (2\pi\sigma^2)^{-\frac{n}{2}}\exp\left(-\frac{1}{2\sigma^2} \parallel{\bf y}-{\bf X}{\boldsymbol \beta}\parallel^2\right), \] \[ \pi({\boldsymbol \beta}|\sigma^2) = (2\pi\sigma^2\nu)^{-\frac{p}{2}}\exp\left(-\frac{1}{2\sigma^2\nu} \parallel{\boldsymbol \beta}\parallel^2\right). \]

Generate 30 samples from \[ y = 1 + 1\times x_1 + 2 \times x_2 + 0 \times x_3 + 0 \times x_4 + 10\times x_5 + \epsilon \] where \(\epsilon \sim \mathcal{N}(0,2),x_i \sim \mathcal{N}(0,1).\beta=(1,1,2,0,0,10)\)

1. use frequentist way method to estimate \(\beta\)

2. use Bayesian MC

3. use Bayesian MCMC

Computation Algorithm

• In frequentist way, that is OLS method, we know the estimate of \({\boldsymbol \beta}\) is \(\hat {\boldsymbol \beta}= ({\bf X}^{{\bf T}}{\bf X})^{-1}{\bf X}^{{\bf T}}{\bf y}\). The mean square of error is \(\parallel\hat{\boldsymbol \beta}-{\boldsymbol \beta}\parallel^2/6\).

• In Bayesian MC way, we assume that

\(\sigma^2 \sim \mathcal{IG}(a_0, b_0)\)

\({\boldsymbol \beta}|\sigma^2 \sim \mathcal{N}(0, \sigma^2\nu I_p)\)

Then we have proved that

\(\sigma^2|{\bf y}\sim\mathcal{IG}(a^*, b^*)\), where \(a^* = \frac{n}{2}+ a_0,b^*=b_0+\frac{1}{2}{\bf y}^{{\bf T}} (I_n-{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}\)

\({\boldsymbol \beta}|{\bf y},\sigma^2 \sim\mathcal{N}(A^{-1}{\bf X}^{{\bf T}}{\bf y}, {\bf A}^{-1}\sigma^2)\), where \({\bf A}=({\bf X}^{{\bf T}}{\bf X}+ \frac{1}{\nu}I_p)\)

Finally, we estimate \(\hat{\boldsymbol \beta}\) by computing \(E({\boldsymbol \beta}|{\bf y})=\sum_{m=1}^{M} {\boldsymbol \beta}^{(m)} /M\). The mean square of error is \(\parallel\hat{\boldsymbol \beta}-{\boldsymbol \beta}\parallel^2/6\).

• In Bayesian MCMC way, we have proved that

• 1. \(\sigma^2|{\bf y},{\boldsymbol \beta}\sim\mathcal{IG}(a^*, b^*)\), where \(a^* = \frac{n+p}{2}+ a_0,b^*=b_0+\frac{1}{2}\parallel {\bf y}- {\bf X}{\boldsymbol \beta}\parallel^2+ \frac{1}{2\nu}\parallel {\boldsymbol \beta}\parallel^2\)
• 2. \({\boldsymbol \beta}|{\bf y},\sigma^2 \sim\mathcal{N}({\bf A}^{-1}{\bf X}^{{\bf T}}{\bf y}, {\bf A}^{-1}\sigma^2)\),where \({\bf A}=({\bf X}^{{\bf T}}{\bf X}+ \frac{1}{\nu}I_p)\)
• 3. \(\beta_j|{\boldsymbol \beta}_{-j},\sigma^2,{\bf y}\sim \mathcal{N}\left((x_j^{{\bf T}}x_j +\frac{1}{\nu})^{-1}x_j^{{\bf T}}({\bf y}-{\bf X}_{-j}{\boldsymbol \beta}_{-j}),\sigma^2\left[x_j^{{\bf T}}x_j +\frac{1}{\nu}\right]^{-1}\right)\)

Then, we use Gibbs sampling method to get the \({\boldsymbol \beta}^{(m)}=(\beta_0^{(m)},\beta_1^{(m)},\ldots,\beta_5^{(m)})\). Set the initial values \((\beta_0^{(0)},\beta_1^{(0)},\ldots,\beta_5^{(0)},\sigma^{2(0)})\).Then iterate the following steps for \(t=0,1,2,\ldots\) as follows:

• 1. Generate \(\sigma^{2(t+1)}\) from \(\pi(\sigma^2|{\bf y},\beta_0^{(t)},\beta_1^{(t)},\ldots,\beta_5^{(t)})\).
• 2. Generate \(\beta_0^{(t+1)}\) from \(\pi(\beta_0|{\bf y},{\boldsymbol \beta}_{-0}^{(t)},\sigma^{2(t+1)})\).
• 3. update \({\boldsymbol \beta}^{(t)}\) with \({\boldsymbol \beta}^{(t)}=(\beta_0^{(t+1)},\beta_1^{(t)},\ldots,\beta_5^{(t)})\)
• 4. Generate \(\beta_1^{(t+1)}\) from \(\pi(\beta_1|{\bf y},{\boldsymbol \beta}_{-1}^{(t)},\sigma^{2(t+1)})\).
• 5. update \({\boldsymbol \beta}^{(t)}\) with \({\boldsymbol \beta}^{(t)}=(\beta_0^{(t+1)},\beta_1^{(t+1)},\beta_2^{(t)},\ldots,\beta_5^{(t)})\) and similarly, generate \(\beta_2^{(t+1)},\ldots,\beta_5^{(t+1)}\) and update \({\boldsymbol \beta}^{(t)}\) after generating \(\beta_i^{(t+1)}\) every time,

where \({\boldsymbol \beta}_{-i}^{(t)}\) is a subvector of \({\boldsymbol \beta}\) without \(i^{th}\) element.

We setup the total sample size is 5000, and burning period is 2000.

Finally, we estimate \(\hat{\boldsymbol \beta}\) by computing \(E({\boldsymbol \beta}|{\bf y})=\sum_{m=2001}^{M} {\boldsymbol \beta}^{(m)} /M\). The mean square of error is \(\parallel\hat{\boldsymbol \beta}-{\boldsymbol \beta}\parallel^2/6\).

Simulation Results

Model Selection

Consider we have k candidate models, \(\mathcal{M}=\{M_1, M_2, \ldots, M_k\}\). A Bayesian way to select the best model is to find the model which maximizes the model posterior probability given data, that is

\[\hat M = \arg \max_{M\in\mathcal{M}}\pi(M|{\bf y})\] where \[\pi(M|{\bf y}) = \frac{f({\bf y}|M)\pi(M)}{\sum_{\tilde M\in\mathcal{M}}f({\bf y}|\tilde M)\pi(\tilde M)}\] in which, the denominator is free of \(M\). note that the prior \(\pi(M)=1/k\). Then we have \[\pi(M|{\bf y}) = \frac{f({\bf y}|M)}{\sum_{w\in\mathcal{M}}f({\bf y}|\tilde M)}\] This leads to \[\hat M = \arg\max_{M\in\mathcal{M}} f({\bf y}|M)\] where \(f({\bf y}|M) = \int f({\bf y}|{\boldsymbol \theta}_M,M)\pi({\boldsymbol \theta}_M)d{\boldsymbol \theta}_M\), where \({\boldsymbol \theta}_M\) is the parameter under model \(M\). We call \(f({\bf y}|M)\) as marginal likelihood under model \(M\),\(f({\bf y}|{\boldsymbol \theta}_M,M)\) likelihood under model \(M\) and \(\pi({\boldsymbol \theta}_M)\) prior under model M. Under \(M\) \[\begin{eqnarray*} f({\bf y}|M) &=& \int f({\bf y}|{\boldsymbol \theta},M)\pi({\boldsymbol \theta}|M)d{\boldsymbol \theta}\\ &=& \int\int f({\bf y}|{\boldsymbol \beta},\sigma^2,M)\pi({\boldsymbol \beta}|\sigma^2,M)\pi(\sigma^2|M)d{\boldsymbol \beta}d\sigma^2\\ &=& \int \int (2\pi)^{-\frac{n}{2}}(\sigma^2)^{-\frac{n}{2}}\exp\left[-\frac{1}{2\sigma^2}({\bf y}- {\bf X}{\boldsymbol \beta}_M)^{{\bf T}}({\bf y}-{\bf X}{\boldsymbol \beta}_M)\right]\times\\ &&\quad\quad (2\pi)^{-\frac{p_M}{2}}(\sigma^2\nu)^{-\frac{p_M}{2}}\exp(-\frac{1}{2\sigma^2\nu}{\boldsymbol \beta}_M^{{\bf T}}{\boldsymbol \beta}_M)\times\\ &&\quad\quad \frac{b_1^{a_1}}{\Gamma(a_1)}(\sigma^2)^{(-a_1-1)}\exp(-\frac{b_1}{\sigma^2})d{\boldsymbol \beta}_M d\sigma^2\\ &=& \frac{b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p_M}{2}}}(2\pi)^{(-\frac{n+p_M}{2})}\int \int \sigma^{2(\frac{-(n+p_M+2a_1+2)}{2})}\times\\ &&\exp\left[-\frac{1}{2\sigma^2} \left({\bf y}^T{\bf y}- 2{\boldsymbol \beta}_M^T{\bf X}^T{\bf y}+ {\boldsymbol \beta}_M^T{\bf X}^T{\bf X}{\boldsymbol \beta}_M + \frac{1}{\nu}{\boldsymbol \beta}_M^T{\boldsymbol \beta}_M +2b_1 \right)\right]d{\boldsymbol \beta}_M d\sigma^2\\ \end{eqnarray*}\] Note that \[\begin{eqnarray*} &&{\bf y}^T{\bf y}- 2{\boldsymbol \beta}_M^T{\bf X}^T{\bf y}+ {\boldsymbol \beta}_M^T{\bf X}^T{\bf X}{\boldsymbol \beta}_M + \frac{1}{\nu}{\boldsymbol \beta}_M^T{\boldsymbol \beta}_M +2b_1\\ &=&{\boldsymbol \beta}_M^{{\bf T}} \left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}{\bf I}_{p_M}\right){\boldsymbol \beta}_M - 2 {\boldsymbol \beta}_M^{{\bf T}} \left({\bf X}^{{\bf T}} {\bf X}+ \frac{1}{\nu}{\bf I}_{p_M}\right)\left({\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}{\bf I}_{p_M}\right)^{-1} {\bf X}^{{\bf T}} {\bf y}+ {\bf y}^T{\bf y}+2b_1\\ &=&{\boldsymbol \beta}_M^{{\bf T}}{\bf A}{\boldsymbol \beta}_M - 2{\boldsymbol \beta}_M^{{\bf T}}{\bf A}\tilde{{\boldsymbol \beta}_M} + \tilde{{\boldsymbol \beta}_M}^{{\bf T}} {\bf A}\tilde{{\boldsymbol \beta}_M} - \tilde{{\boldsymbol \beta}_M}^{{\bf T}}{\bf A}\tilde{{\boldsymbol \beta}_M}+ {\bf y}^{{\bf T}}{\bf y}+2b_1\\ &=& ({\boldsymbol \beta}_M- \tilde{{\boldsymbol \beta}_M})^T{\bf A}({\boldsymbol \beta}-\tilde{{\boldsymbol \beta}_M}) + {\bf y}^{{\bf T}}{\bf y}- {\bf y}^{{\bf T}}{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}{\bf y}+ 2b_1\\ &=& ({\boldsymbol \beta}_M- \tilde{{\boldsymbol \beta}_M})^T{\bf A}({\boldsymbol \beta}_M-\tilde{{\boldsymbol \beta}_M}) + {\bf y}^{{\bf T}}({\bf I}_n-{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}+ 2b_1 \end{eqnarray*}\] Therefore, \[\begin{eqnarray*} f({\bf y}|M) &=&\frac{b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p{_M}}{2}}}(2\pi)^{\left(-\frac{n+p{_M}}{2}\right)}\int \int \sigma^{2\left(\frac{-(n+p_M+2a_1+2)}{2}\right)}\times \\ &&\exp\left[-\frac{1}{2\sigma^2}\left(({\boldsymbol \beta}_M- \tilde{{\boldsymbol \beta}_M})^T{\bf A}({\boldsymbol \beta}_M-\tilde{{\boldsymbol \beta}_M}) + {\bf y}^{{\bf T}}({\bf I}_n-{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}+ 2b_1\right)\right]d{\boldsymbol \beta}_M d\sigma^2\\ &=& \frac{b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p_M}{2}}}(2\pi)^{\left(-\frac{n+p_M}{2}\right)}\int\sigma^{2\left(\frac{-(n+p_M+2a_1+2)}{2}\right)}\times\\ &&(2\pi)^{\frac{p_M}{2}}\lvert{\bf A}\rvert^{-\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^2}\left( {\bf y}^{{\bf T}}({\bf I}_n-{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}+ 2b_1\right)\right]d\sigma^2\\ &=& \frac{b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p_M}{2}}}(2\pi)^{(-\frac{n}{2})}\lvert{\bf A}\rvert^{-\frac{1}{2}}\int \sigma^{2(-a_M^*-1)}\exp(-\frac{b_M^*}{\sigma^2})d\sigma^2\\ &=& \frac{b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p_M}{2}}}(2\pi)^{(-\frac{n}{2})}\lvert{\bf A}\rvert^{-\frac{1}{2}} \Gamma(a_M^*)b_M^{*(-a_M^*)}\\ &=&(2\pi)^{(-\frac{n}{2})}\frac{\Gamma(a_M^*)b_M^{*(-a_M^*)}b_1^{a_1}}{\Gamma(a_1)\nu^{\frac{p_M}{2}}}\lvert{\bf A}\rvert^{-\frac{1}{2}} \end{eqnarray*}\] where \[\begin{eqnarray*} {\bf A}&=&{\bf X}^{{\bf T}} {\bf X}+\frac{1}{\nu}{\bf I}_{p_M}, \quad\tilde{{\boldsymbol \beta}_M} = {\bf A}^{-1} {\bf X}^{{\bf T}} {\bf y}\\ a_M^* &=& \frac{n+p_M}{2}+a_1,\quad b_M^*=\frac{1}{2}{\bf y}^{{\bf T}}({\bf I}_n-{\bf X}{\bf A}^{-1}{\bf X}^{{\bf T}}){\bf y}+ b_1 \end{eqnarray*}\] After Laplace approximation, \(f({\bf y}|M)\) is equivalent to BIC for large n.

MCMC

Minghui Chen (IWDE)

Goal is to compute \[ m({\bf y}) = \int f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})d{\boldsymbol \theta} \]

\[\begin{eqnarray} \pi({\boldsymbol \theta}|{\bf y})=\frac{f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})}{\int f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})d{\boldsymbol \theta}} \end{eqnarray}\] Let \(g({\boldsymbol \theta})\) be a density function of \({\boldsymbol \theta}\)

Then,

\[\begin{eqnarray*} &&\int g({\boldsymbol \theta})d{\boldsymbol \theta}= 1 \\ \Leftrightarrow && \int g({\boldsymbol \theta})\frac{\pi({\boldsymbol \theta}|{\bf y})}{\pi({\boldsymbol \theta}|{\bf y})}d{\boldsymbol \theta}=1\\ \Leftrightarrow && \int g({\boldsymbol \theta})\frac{\pi({\boldsymbol \theta}|{\bf y})}{\frac{f({\bf y}|{\boldsymbol \theta}) \pi({\boldsymbol \theta})}{m({\bf y})}}d{\boldsymbol \theta}=1\\ \Leftrightarrow && m({\bf y})\int \frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})}\pi({\boldsymbol \theta}|{\bf y}) d{\boldsymbol \theta}=1\\ \Leftrightarrow && m({\bf y}) = \left[\int \frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})}\pi({\boldsymbol \theta}|{\bf y}) d{\boldsymbol \theta}\right]^{-1}\\ &&= E\left[\frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})}|{\bf y}\right]^{-1}\\ &&\approx \left[\sum_{m=1}^{M}\frac{g({\boldsymbol \theta}^{(m)})}{f({\bf y}|{\boldsymbol \theta}^{(m)})\pi({\boldsymbol \theta}^{(m)})}/M\right]^{-1} \end{eqnarray*}\] where \({\boldsymbol \theta}^{(m)}\) is from \(\pi({\boldsymbol \theta}|{\bf y})\).

There maxmimum of \(m({\bf y})\) is equivalent to maximize \[ E\left[\frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})}|{\bf y}\right]^{-1} \]

Now let’s use Minghui Chen’s method to maximize the \(f({\bf y}|M_i)\).

\(f({\bf y}|M_i)=\int f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)d{\boldsymbol \theta}\), where \({\boldsymbol \theta}=(\beta_0,\beta_1,\sigma^2)\) \[\begin{eqnarray} \pi({\boldsymbol \theta}|{\bf y},M_i)&=&\frac{f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)}{\int f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)d{\boldsymbol \theta}}\\ &=&\frac{f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)}{f({\bf y}|M_i)} \end{eqnarray}\] Let \(g({\boldsymbol \theta})\) be a density function of \({\boldsymbol \theta}\)

Then

\[\begin{eqnarray*} &&\int g({\boldsymbol \theta})d{\boldsymbol \theta}= 1 \\ \Leftrightarrow && \int g({\boldsymbol \theta})\frac{\pi({\boldsymbol \theta}|{\bf y},M_i)}{\pi({\boldsymbol \theta}|{\bf y},M_i)}d{\boldsymbol \theta}=1\\ \Leftrightarrow && \int g({\boldsymbol \theta})\frac{\pi({\boldsymbol \theta}|{\bf y},M_i)}{\frac{f({\bf y}|{\boldsymbol \theta},M_i) \pi({\boldsymbol \theta}|M_i)}{f({\bf y}|M_i)}}d{\boldsymbol \theta}=1\\ \Leftrightarrow && f({\bf y}|M_i)\int \frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)}\pi({\boldsymbol \theta}|{\bf y},M_i) d{\boldsymbol \theta}=1\\ \Leftrightarrow && f({\bf y}|M_i) = \left[\int \frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)}\pi({\boldsymbol \theta}|{\bf y},M_i) d{\boldsymbol \theta}\right]^{-1}\\ &&= E\left[\frac{g({\boldsymbol \theta})}{f({\bf y}|{\boldsymbol \theta},M_i)\pi({\boldsymbol \theta}|M_i)}|{\bf y},M_i\right]^{-1}\\ &&\approx \left[\sum_{n=1}^{N}\frac{g({\boldsymbol \theta}^{(n)})}{f({\bf y}|{\boldsymbol \theta}^{(n)},M_i)\pi({\boldsymbol \theta}^{(n)}|M_i)}/N\right]^{-1} \end{eqnarray*}\] where \({\boldsymbol \theta}^{(n)}\) is from \(\pi({\boldsymbol \theta}|{\bf y},M_i), i=1,2,\ldots, k\).

The issue here is a choice of \(g({\boldsymbol \theta})\) to minimize MC error.

Recall \[ \pi({\boldsymbol \theta}|{\bf y}) \propto f({\bf y}|{\boldsymbol \theta})\pi({\boldsymbol \theta})\approx g({\boldsymbol \theta}) \] Let \(g({\boldsymbol \theta})\) be a normal density such that \(E_g({\boldsymbol \theta})=\) posterior mean \(\approx\sum_{m=1}^{M}{\boldsymbol \theta}^{(m)}/M = \bar {\boldsymbol \theta}\), \(V_g({\boldsymbol \theta})=\) posterior variance matrix \(\approx \sum_{m=1}^{M}({\boldsymbol \theta}^{(m)}-\bar {\boldsymbol \theta})({\boldsymbol \theta}^{(m)}-\bar {\boldsymbol \theta})^{{\bf T}}/M\), where \({\boldsymbol \theta}^{(m)}\) is a MCMC sample from \(\pi({\boldsymbol \theta}|{\bf y})\).

For given model M, using MCMC, generate \({\boldsymbol \theta}_M^{(k)}\) from \(\pi({\boldsymbol \theta}_M|{\bf y}, M)\), then compute \(\bar {\boldsymbol \theta}_M\) and \(V_M\).

Define \(g({\boldsymbol \theta}|M) = \mathcal{N({\boldsymbol \theta};\bar {\boldsymbol \theta}_M, V_M)}\).

Compute

\[\begin{eqnarray} \tag{4}f({\bf y}|M)\approx \left[\sum_{k=1}^{K}\frac{g({\boldsymbol \theta}^{(k)}|M)}{f({\bf y}|{\boldsymbol \theta}^{(k)},M)\pi({\boldsymbol \theta}^{(k)}|M)}/K\right]^{-1} \end{eqnarray}\]

Bayesian Simulation Case Study

Setup

• Consider

\[ f({\bf y}|{\boldsymbol \beta},\sigma^2) = (2\pi\sigma^2)^{-\frac{n}{2}}\exp\left(-\frac{1}{2\sigma^2} \parallel{\bf y}-{\bf X}{\boldsymbol \beta}\parallel^2\right), \] \[ \pi({\boldsymbol \beta}|\sigma^2) = (2\pi\sigma^2\nu)^{-\frac{p}{2}}\exp\left(-\frac{1}{2\sigma^2\nu} \parallel{\boldsymbol \beta}\parallel^2\right). \]

Generate 30 samples from \[ y = 1 + 1\times x_1 + 2 \times x_2 + 0 \times x_3 + 0 \times x_4 + 10\times x_5 + \epsilon \] where \(\epsilon \sim \mathcal{N}(0,2),x_i \sim \mathcal{N}(0,1).\beta=(1,1,2,0,0,10)\)

Therefore, the true model, say \(M_{ture}\), is \[ y = 1 + 1\times x_1 + 2 \times x_2 + 10\times x_5 + \epsilon \]

• 1. use Bayesian MC method to find true model
• 2. use frequentist way method to find true model

Our candidate models are as follows: \[ M_1: y = \beta_0 + \epsilon; \] \[ M_2: y = \beta_0 + \beta_1X_1+ \epsilon; \] \[ \ldots \] \[ M_{32}: y = 1 + 1\times x_1 + 2 \times x_2 + 0 \times x_3 + 0 \times x_4 + 10\times x_5 + \epsilon \]

• For bayesian way, compute

\(f({\bf y}|M_k)\) for \(k = 1,2,\ldots, 32\) and select the best model \(\hat M\), such that \[ f({\bf y}|\hat M) \geq f({\bf y}| M_k) \] If \(\hat M = M_{true}\), then assign \(Z_i = 1\), otherwise \(Z_i = 0\).

Repeat the above procedure for \(R=500\) times, and compute \(\sum_{i = 1}^{R} \frac{Z_i}{R} \simeq P(\) selecting true model with Bayes method \()\).

• For Frequentist way, similarly, build linear model for 32 candidate models, compute AIC and BIC. If the smallest AIC or BIC is ture model, then \(W_i = 1\), otherwise \(W_i = 0\). Repeat it 500 times and compute \(\sum_{i = 1}^{R} \frac{W_i}{R}\).

Consider the sample sizes are 20, 50, 100 and 300 so that we can see how sample size impacts the probability of selecting the true model.

• For MCMC, sampling with MCMC method, then compute Eq. (4) for 32 candidates models, select the best model \(\hat M\), such that \[ f({\bf y}|\hat M) \geq f({\bf y}| M_k) \] If \(\hat M = M_{true}\), then assign \(Z_i = 1\), otherwise \(Z_i = 0\)

Repeat the above procedure for \(R=500\) times, and compute \(\sum_{i = 1}^{R} \frac{Z_i}{R} \simeq P(\) selecting true model with MCMC method \()\)

Result