Bayes Factor with Lindley Paradox and Tow Standard Methods in Model

Xiaoting Nie

Changsha Normal University, Changsha, China.

DOI: 10.4236/ojs.2020.101006                

Abstract

For any statistical analysis, Model selection is necessary and required. In many cases of selection, Bayes factor is one of the important basic elements. For the unilateral hypothesis testing problem, we extend the harmony of frequency and Bayesian evidence to the generalized p-value of unilateral hypothesis testing problem, and study the harmony of generalized P-value and posterior probability of original hypothesis. For the problem of single point hypothesis testing, the posterior probability of the Bayes evidence under the traditional Bayes testing method, that is, the Bayes factor or the single point original hypothesis is established, is analyzed, a phenomenon known as the Lindley paradox, which is at odds with the classical frequency evidence of p-value. At this point, many statisticians have been worked for this from both frequentist and Bayesian perspective. In this paper, I am going to focus on Bayesian approach to model selection, starting from Bayes factors and going within Lindley Paradox, which also briefly talks about partial and fractional Bayes factor. Trying to use a simple way to consider this paradox is the thing what I want to do in the paper. On the other hand, a detailed derivation of BIC and AIC is given in Section 4. The guiding principle of selecting the optimal model is to investigate from two aspects: one is to maximize the likelihood function, the other is to minimize the number of unknown parameters in the model. The larger the likelihood function value, the better the model fitting, but we can not simply measure the model fitting accuracy, which leads to more and more unknown parameters in the model, and the model that becomes more and more complex would have caused an overmatch. Therefore, a good model should be the combination of the fitting accuracy and the number of unknown parameters to optimize the configuration.

Keywords

Bayes Factors, Lindley Paradox, Fractional Bayes Factor, BIC, AIC

1. Introduction

Many statisticians are naturally involved in the question of model selection [1], in case to define the “best model” to fit real data, different approaches have been proposed since last century, many well-known methods such as F-test [2], AIC, BIC [3], Bayesian model averaging [4]. We are focusing on Bayesian approach, as we analyze data from some possible models $M_1\mathrm{,}\cdots \mathrm{,}{M}_k$. We denote $\theta \in \Theta$ as parameter and $\pi \left(\theta \right)$ as prior probability, then for likelihoods $f\left(Data\mathrm{<mo>\; |\; </mo>}\theta \right)$ and prior $g\left(\theta \right)$. The posterior for model $M_k$ with parameter ${\theta }_k$ is proportional to $f_k\left(Data\mathrm{<mo>\; |\; </mo>}{\theta }_k\right)g_k\left({\theta }_k\right){\pi }_k$, we get posterior probability as

$\begin{array}{c}P\left(M_k|Data\right)\propto {\pi }_k{\displaystyle {\int }_{{\Omega }_k}}{f}_k\left(Data|{\theta }_k\right)g_k\left({\theta }_k\right)\text{d}{\theta }_k\\ =\frac{{\pi }_k{\displaystyle {\int }_{{\Omega }_k}}{f}_k\left(Data|{\theta }_k\right)g_k\left({\theta }_k\right)\text{d}{\theta }_k}{{\displaystyle {\sum }_{j=1}^K{\pi }_j{\displaystyle {\int }_{{\Omega }_j}}{f}_j\left(Data|{\theta }_j\right)g_j\left({\theta }_j\right)\text{d}{\theta }_j}}\end{array}$

In a Bayesian analysis, the priors ${\pi }_k$ on each model and $g_k\left({\theta }_k\right)$ on the parameters of model k are proper and subjective. And the Bayesian solutions to do questions are to compute the posterior probability $P\left(M_k|Data\right)$ for each model. For model selection, we would choose the model from Bayesian conclusion as maximizes $P\left(M_k|Data\right)$.

However, Bayes factor has its only limitation, that is Bayes factors itself can only show the difference of how hypothesis model is against a null model [5]. Also, Bayes factor has a close connection with priors, if we change the width of the prior, it will also change the Bayes factor. At this point, we may need to consider about Lindley Paradox.

In Section 2, we give a simple and general explanation of Bayes factor. Following, in Section 3, we will talk about Lindley’s Paradox. And Section 4 can be one of the main parts of the theoretical approach for AIC and BIC, for which we give the derivation. A simple example is given as well to use AIC and BIC.

2. Bayes Factor

Before talking about all things, first we would construct one of the most important variables within Bayesion Methods-Bayes Factor [6].

Suppose we have data D with prior $\theta$ and $M_1\mathrm{,}{M}_2$ as two different models. By Condition Rule, we have:

$P\left(M_1|D\right)=\frac{P\left(D|M_1\right)P\left(M_1\right)}{P(D)}$

Recall for Odds we have $P\left(M_1\right)=1-P\left(M_2\right)$. And for $P\left(D|M_1\right)$ is the marginal likelihood, which $P\left(D\mathrm{<mo>\; |\; </mo>}{M}_1\right)={\displaystyle \int }P\left(D\mathrm{<mo>\; |\; </mo>}\theta \mathrm{,}{M}_1\right)P\left(\theta \mathrm{<mo>\; |\; </mo>}{M}_1\right)\text{d}\theta$. $\theta$ denotes prior. Then, by Bayes’ Rule,

$P\left(M_1\mathrm{<mo>\; |\; </mo>}D\right)=\frac{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_1\right)P\left(M_1\right)}{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_1\right)P\left(M_1\right)+P\left(D\mathrm{<mo>\; |\; </mo>}{M}_2\right)P(M2)}$

$\frac{P\left(M_1\mathrm{<mo>\; |\; </mo>}D\right)}{P\left(M_2\mathrm{<mo>\; |\; </mo>}D\right)}=\frac{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_1\right)}{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_2\right)}\frac{P\left(M_1\right)}{P(M2)}$

where $\frac{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_1\right)}{P\left(D\mathrm{<mo>\; |\; </mo>}{M}_2\right)}$ is defined as Bayes’ factor, and realized it is also the ratio of marginal likelihood. Furthermore, we denote Bayes’ factor as:

$B_{\mathrm{1,2}}\left(y\right)=\frac{P\left(y\mathrm{<mo>\; |\; </mo>}{M}_1\right)}{P\left(y\mathrm{<mo>\; |\; </mo>}{M}_2\right)}$

Bayesian method fits in many models for testing because it can provide a decisiveness of the evidence agree the null model in contrast p-values [7] which are usually just regarded as evidence mearsurement against the alternative [8]. Also, the Bayes factor (Jefferys, 1961) [9] is used in Bayesian hypothesis. Assuming that $p_i\left(D|{\theta }_i\right),i=1,2$ are the likelihoods for D under two competing models $H_1$ and $H_2$, and the parameters are ${\theta }_i$. Meanwhile, let ${\pi }_i\left({\theta }_i\right)$ be their prior distributions [10]. The Bayes factor for $H_2$ against $H_1$ :

$B_{\mathrm{2,1}}=\frac{{\displaystyle \int }{P}_2\left(D\mathrm{<mo>\; |\; </mo>}{\theta }_2\right){\pi }_2\left({\theta }_2\right)\text{d}{\theta }_2}{{\displaystyle \int }{P}_1\left(D\mathrm{<mo>\; |\; </mo>}{\theta }_1\right){\pi }_1\left({\theta }_1\right)\text{d}{\theta }_1}$

Above these, evidence from the data agrees $H_2$, against $H_1$. So Bayes factor can avoid many limitations in p-value testing. The development of Bayes factor in statistical models test can applicate in many areas of research [11].

3. Priors and Lindley Paradox

3.1. Introduction to Lindley Paradox

The Lindley’s Paradox shows how a value (or the number of standard deviations) is used in a Frequent Assumption [12] test results in a completely different inference from Bayesian hypothesis [13].

When we faced with improper priors (like priors can’t integrate to one) in the null hypothesis and model selection, we will find some problems. Such priors can be acceptable, but for other purposes it is also acceptable. So we consider testing the hypotheses:

$H_0\mathrm{<mo>\; :\; </mo>}\theta \in {\Theta }_0\text{vs}{H}_1\mathrm{<mo>\; :\; </mo>}\theta \in {\Theta }_1$

Defining $\theta$ for marginal density, so we can use the following model:

$p\left(\theta \right)=p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_0\right)p\left(H_0\right)+p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_1\right)p(H1)$

Making $p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_0\right)$ and $p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_1\right)$ are proper density functions, the posterior is given by:

$\begin{array}{l}p\left(H_0\mathrm{<mo>\; |\; </mo>}D\right)\\ =\frac{p\left(H_0\right)p\left(D\mathrm{<mo>\; |\; </mo>}{H}_0\right)}{p\left(H_0\right)P\left(D\mathrm{<mo>\; |\; </mo>}{H}_0\right)+p\left(H_1\right)p\left(D\mathrm{<mo>\; |\; </mo>}{H}_1\right)}\\ =\frac{p\left(H_0\right){\displaystyle {\int }_{{\Theta }_0}}p\left(D|\theta \right)p\left(\theta |H_0\right)\text{d}\theta }{p\left(H_0\right){\displaystyle {\int }_{{\Theta }_0}}p\left(D|\theta \right)p\left(\theta |H_0\right)\text{d}\theta +p\left(H_1\right){\displaystyle {\int }_{{\Theta }_1}}p\left(D|\theta \right)p\left(\theta |H_1\right)\text{d}\theta }\end{array}$

Then we can suppose that we use improper priors, making $p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_0\right)\propto z_0$ and $p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_1\right)\propto z_1$. So:

$\begin{array}{c}p\left(H_0|D\right)=\frac{p\left(H_0\right)z_0{\displaystyle {\int }_{{\Theta }_0}}p\left(D|\theta \right)\text{d}\theta }{p\left(H_0\right)z_0{\displaystyle {\int }_{{\Theta }_0}}p\left(D|\theta \right)\text{d}\theta +p\left(H_1\right)z_1{\displaystyle {\int }_{{\Theta }_1}}p\left(D|\theta \right)\text{d}\theta }\\ =\frac{p\left(H_0\right)z_0{s}_0}{p\left(H_0\right)z_0{s}_0+p\left(H_1\right)z_1{s}_1}\end{array}$

Establishing model i that $s_i={\displaystyle {\int }_{{\Theta }_i}}p\left(x|\theta \right)\text{d}\theta$ is the marginal likelihood or the integrated. So we assume that $p\left(H_0\right)=p\left(H_1\right)=\frac{1}{2}$

Then an equation can be obtained:

$p\left(H_0|D\right)=\frac{z_0{s}_0}{z_0{s}_0+z_1{s}_1}=\frac{s_0}{s_0+\left(\frac{z_1}{z_0}\right)s_1}$

So we can use different z that we want to change the posterior arbitrarily. Meanwhile, when using proper and not clear priors might cause similar problems. Because the probability of data in a complex model with a diffuse prior will be very small. So one thing we must know, when we do research in Bayes factor a clearer and simper model is better. It was called the Lindley paradox.

3.2. A Simple Model in Lindley Paradox

Many authors [14] have discussed this so-called paradox [15] in different ways [16]. So I want to find a simple way to consider this problem. The usual point null hypothesis testing problem is to test:

$H_0:\theta =0\text{vs}{H}_1:\theta \ne 0$

In normal model $N\left(x\mathrm{<mo>\; |\; </mo>}\theta \mathrm{,1}\right)$. The prior probability is $p\left(H_0\right)\left(p\left(H_0\right)={\rho }_0\right)$.

Let $\pi \left(\theta \right)=N\left(\theta |0,{\sigma }^2\right)\left(\sigma >0\right)$ be the prior distribution for the unknowm parameter $\theta$ in the model.

The Bayes factor is given by:

$B=\frac{N\left(x|0,1\right)}{{\displaystyle \int }N\left(x|\theta ,1\right)\pi \left(\theta \right)\text{d}\theta }$

In order to consider the paradox, we can formalise it and compare the two following normal models:

$M_0=\left\{N\left(x|0,1\right)={\left(2\pi \right)}^{-\frac{1}{2}}\exp \left(-\frac{1}{2}{x}^2\right)\right\}$

$M_1=\left\{N\left(x|\theta ,1\right)={\left(2\pi \right)}^{-\frac{1}{2}}\exp \left\{-\frac{1}{2}{\left(x-\theta \right)}^2\right\}\right\}$

Consider a physical system where quantity X may be measured and assume. And we need to use the $\sigma$ to define both the priors. The prior of the null hypothesis is ${\rho }_0$ supposing the ${\rho }_0$ can depend on $\sigma$.

Computing the Bayes factor representing the odds of the null hypothesis $H_0$ is:

$B_0=\frac{N\left(x|0,1\right)}{{\displaystyle \int }N\left(x|\theta ,1\right)N\left(\theta |0,{\sigma }^2\right)\text{d}\theta }=\frac{{\text{e}}^{-\frac{1}{2}{x}^2}}{{\text{e}}^{\frac{-\frac{1}{2}{x}^2}{{\sigma }^2+1}}}\sqrt{{\sigma }^2+1}$

In this case, prior probabilities $p\left(H_0\right)$ and $p\left(H_1\right)=1-P\left(H_0\right)$ for two hypotheses can be expressed. Given the result x, in Bayes theory that:

$p\left(H_m\mathrm{<mo>\; |\; </mo>}x\right)p\left(x\right)=p\left(x\mathrm{<mo>\; |\; </mo>}{H}_m\right)p(Hm)$

for $m=0,1$, $p\left(H_m\right)$ is prior probabilities and $p\left(x\mathrm{<mo>\; |\; </mo>}{H}_m\right)$ is the conditional distribution, $p\left(x\right)=p\left(x\mathrm{<mo>\; |\; </mo>}{H}_0\right)P\left(H_0\right)+P\left(x\mathrm{<mo>\; |\; </mo>}{H}_1\right)p\left(H_1\right)$ can outcome the overall distribution. Posterior probability $p\left(H_m\mathrm{<mo>\; |\; </mo>}x\right)$ is in the hypothesis $H_m$. In Bayes theory we can evaluate the posterior probabilities, $p\left(H_0\mathrm{<mo>\; |\; </mo>}x\right)$ is given by:

$\begin{array}{c}p\left(H_0|x\right)=\frac{p\left(x|H_0\right)p\left(H_0\right)}{p\left(x\right)}\\ =\frac{p\left(x|H_0\right)p\left(H_0\right)}{p\left(x|H_0\right)p\left(H_0\right)+p\left(x|H_1\right)p\left(H_1\right)}\\ ={\left[1+\frac{p\left(x|H_1\right)p\left(H_1\right)}{p\left(x|H_0\right)p\left(H_0\right)}\right]}^{-1}\\ ={\left[1+\frac{1-p\left(H_0\right)}{p\left(H_0\right)}\frac{p\left(x|H_1\right)}{p\left(x|H_0\right)}\right]}^{-1}\end{array}$

Then, we can use the mean value in prior distribution with ${\pi }_m\left(\theta \right)\equiv p\left(\theta \mathrm{<mo>\; |\; </mo>}{H}_m\right)$ and make the rest of the prior probability as a normal distribution with variance $\tau$, so:

${\pi }_1\left(\theta \right)={\left(2\pi {\tau }^2\right)}^{-\frac{1}{2}}\exp \left\{\frac{-{\theta }^2}{2{\tau }^2}\right\}$

Evaluating the conditional probabilities:

$p\left(x|H_m\right)={\displaystyle \int }{\pi }_m\left(\theta \right)p\left(x|\theta \right)\text{d}\theta$

We can evaluate $p\left(x\mathrm{<mo>\; |\; </mo>}{H}_0\right)$ and $p\left(x\mathrm{<mo>\; |\; </mo>}{H}_1\right)$, overall:

$p\left(H_0|x\right)={\left[1+\frac{1-{\rho }_0}{{\rho }_0}\frac{{\text{e}}^{-\frac{1}{2}{x}^2/\left({\sigma }^2+1\right)}}{{\text{e}}^{-\frac{1}{2}{x}^2}}\frac{1}{\sqrt{{\sigma }^2+1}}\right]}^{-1}={\left[1+\frac{1-{\rho }_0}{{\rho }_0}\frac{1}{B_0}\right]}^{-1}$

So we have an equation like before, we can talk about the prior $\rho \left(\sigma \right)$. Our approach is to measure the value of alternative assumptions about zero. In Asymptotically Bayesian attribute, if the model is incorrectly specified, the posterior will accumulates in the model. In the case of the Kullback-Leibler divergence, the closest to the real model [17]. As a result, divergence $D_{K}L\left(N\mathrm{<mo>\; |\; </mo>}x\left(\theta \mathrm{,1}\right)\parallel N\left(x\mathrm{<mo>\; |\; </mo>0,1}\right)\right)$ represents the loss. Because we know the prior before. The excepted loss can be given:

${\displaystyle {\int }_{\Theta }}{D}_{K}L\left(N\left(x\theta \mathrm{,1}\right)\parallel N\left(x\mathrm{<mo>\; |\; </mo>0,1}\right)\right)\pi \left(\theta \right)\text{d}\theta ={\displaystyle \int }\frac{1}{2}{\theta }^2\pi \left(\theta \right)\text{d}\theta =\frac{1}{2}{\sigma }^2$

The model prior represent the loss relatied with a probability statement, it also determined self-information loss function. So we have the prior on the alternative model is:

$1-{\rho }_0\left(\sigma \right)\propto {\text{e}}^{{\sigma }^2}$

The prior of the null hypothesis is $\rho \left(\sigma \right)\propto 1$, then we can get:

${\rho }_0\left(\sigma \right)=\frac{1}{1+\exp \left\{\frac{1}{2}{\sigma }^2\right\}}$

Then, this applies to the category of large $\sigma$ and $p\left(H_0\mathrm{<mo>\; |\; </mo>}x\mathrm{,}\sigma \right)$ goes to zero, so $p\left(H_0\right)\to 0$. Therefore, this method is consistent, we do not advocate the choice of big $\sigma$.

4. BIC and AIC

4.1. BIC

4.1.1. Notation (Table 1)

Table 1. Notation 1.

4.1.2. Derivation of BIC

In this section we are going to talk about the basic idea [18] of how BIC (Bayesian information criterion) constructed and given the derivation of BIC [4].

As what we have showed in section one, $B_{\mathrm{1,2}}\left(y\right)=\frac{P\left(y\mathrm{<mo>\; |\; </mo>}{M}_1\right)}{P\left(y\mathrm{<mo>\; |\; </mo>}{M}_2\right)}$ as Bayes factor for two models, then we consider more models $M_i$ which $i\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$

$P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)={\displaystyle \int }f\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_i\right)={\displaystyle \int }\exp \left(\log \left(f\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_i\right)g_i\left({\theta }_i\right)\right)\right)\text{d}{\theta }_i$

$f\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_i\right)g_i\left({\theta }_i\right)$ where ${\theta }_i$ is the vector of parameters in the model $M_i$, L is the likelihood function and $g_i\left({\theta }_i\right)$ is the p.d.f. of the distribution of parameters ${\theta }_i$

Denoting ${\tilde{\theta }}_i$ as the posterior mode, then we use Taylor expansion, let $Q\left({\theta }_i\right)=\log \left(f\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_i\right)g_i\left({\theta }_i\right)\right)$, $\begin{array}{c}Q\left({\theta }_i\right)=\log \left(f\left(y|{\theta }_i\right)g_i\left({\theta }_i\right)\right)\\ \approx \log \left(f\left(y|{\tilde{\theta }}_i\right)\right)g_i\left({\tilde{\theta }}_i\right)+\left({\theta }_i-{\tilde{\theta }}_i\right){{\nabla }_{{\theta }_i}Q|}_{{\tilde{\theta }}_i}+\frac{1}{2}{\left({\theta }_i-{\tilde{\theta }}_i\right)}^{\text{T}}{H}_{{\theta }_i}\left({\theta }_i-{\tilde{\theta }}_i\right)\end{array}$.

where $H_{{\theta }_i}$ is a $\left|{\theta }_i\right|\left|{\theta }_i\right|$ matrix such that $H_{mn}={\frac{{\partial }^{2}Q}{\partial {\theta }_m\partial {\theta }_n}|}_{{\tilde{\theta }}_i}$, where $\left|{\theta }_i\right|=d_i=dimension\left({\theta }_i\right)$. since Q attains its maximum, the Hessian matrix $H_{{\theta }_i}$ is negative definite. Let us denote ${\bar{H}}_{{\theta }_i}=-H_{{\theta }_i}$, and then approximate $P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)$ :

$P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)\approx {\displaystyle \int }\exp \left\{{Q|}_{{\tilde{\theta }}_i}+\left({\theta }_i-{\tilde{\theta }}_i\right){{\nabla }_{{\theta }_i}Q|}_{{\tilde{\theta }}_i}+\frac{1}{2}{\left({\theta }_i-{\tilde{\theta }}_i\right)}^{\text{T}}{H}_{{\theta }_i}\left({\theta }_i-{\tilde{\theta }}_i\right)\right\}\text{d}{\theta }_i$

Then, by higher dimension normal distribution,

$\because {\displaystyle \int }\frac{1}{{\left(2\pi \right)}^{\frac{d_i}{2}}{\left|{\tilde{H}}_{{\theta }_i}^{-1}\right|}^{\frac{1}{2}}}\exp \left(-\frac{1}{2}{\left({\theta }_i-{\hat{\theta }}_i\right)}^{{\rm T}}{\stackrel{\to }{H}}_{{\theta }_i}\left({\theta }_i-{\hat{\theta }}_i\right)\right)\text{d}{\theta }_i=1$

$\Rightarrow {\displaystyle \int }\exp \left(-\frac{1}{2}{\left({\theta }_i-{\hat{\theta }}_i\right)}^{{\rm T}}{\stackrel{\to }{H}}_{{\theta }_i}\left({\theta }_i-{\hat{\theta }}_i\right)\right)\text{d}{\theta }_i={\left(2\pi \right)}^{\frac{d_i}{2}}{\left|{\tilde{H}}_{{\theta }_i}^{-1}\right|}^{\frac{1}{2}}$

$\Rightarrow P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)=f\left(y\mathrm{<mo>\; |\; </mo>}{\tilde{\theta }}_i\right)g_i\left({\tilde{\theta }}_i\right){\left(2\pi \right)}^{\frac{d_i}{2}}{\left|{\tilde{H}}_{{\theta }_i}^{-1}\right|}^{\frac{1}{2}}$

$\Rightarrow \log P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)=\log f\left(y\mathrm{<mo>\; |\; </mo>}{\tilde{\theta }}_i\right)+\log g_i\left({\tilde{\theta }}_i\right)+\frac{d_i}{2}\log \left(2\pi \right)+\frac{1}{2}\log \left|{\tilde{H}}_{{\theta }_i}^{-1}\right|$

Furthermore, let us think about Weak Law of Large Numbers. For y is given data, $f\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_i\right)$ is the likelihood $L\left({\theta }_i\mathrm{<mo>\; |\; </mo>}y\right)$ and L attains its maximum at the maximum likelihood estimate ${\theta }_i-{\hat{\theta }}_i$.

We set $g_i\left({\theta }_i\right)=(\begin{array}{cc}\mathrm{1,}& \theta \in \left[{\tilde{\theta }}_i-\frac{1}{2}\mathrm{,}{\tilde{\theta }}_i+\frac{1}{2}\right]\\ \mathrm{0,}& \text{else}\end{array}$, then each element in the matrix, ${\tilde{H}}_{{\theta }_i}$, can be expressed as:

${\tilde{H}}_{mn}=-{\frac{{\partial }^2\log L\left({\theta }_i|y\right)}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}$

Then, for ${\tilde{H}}_{{\theta }_i}$ as a Fisher information matrix that,

$\begin{array}{c}{\tilde{H}}_{mn}=-{\frac{{\partial }^2\log \left({\displaystyle {\prod }_{j=1}^{n}L\left({\theta }_i|y_j\right)}\right)}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\\ =-{\frac{{\partial }^2{\displaystyle {\sum }_{j=1}^n\log L\left({\theta }_i|y_j\right)}}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\\ =-{\frac{{\partial }^2\left(\frac{1}{n}{\displaystyle {\sum }_{j=1}^{n}n\log L\left({\theta }_i|y_j\right)}\right)}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\end{array}$

In this case, for the data $y_1\mathrm{,}\cdots \mathrm{,}{y}_n$ is IID, and n is large, we would apply Weak Law of Large number here, as random variable $X_j=n\log L\left({\theta }_i|y_j\right)$ we have $\frac{1}{n}{\displaystyle {\sum }_{j=1}^n}n\log L\left({\theta }_i\mathrm{<mo>\; |\; </mo>}{y}_j\right)\stackrel{P}{\to }E\left(n\log L\left({\theta }_i\mathrm{<mo>\; |\; </mo>}{y}_j\right)\right)$, Moreover, for Fisher information matrix:

$\begin{array}{c}{\tilde{H}}_{mn}=-{\frac{{\partial }^{2}E\left[n\log L\left({\theta }_i|y_j\right)\right]}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\\ =-{n\frac{{\partial }^{2}E\left[\log L\left({\theta }_i|y_j\right)\right]}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\\ =-{n\frac{{\partial }^{2}E\left[\log L\left({\theta }_i|y_1\right)\right]}{\partial {\theta }_m\partial {\theta }_n}|}_{{\theta }_i={\hat{\theta }}_i}\\ =n{I}_{mn}\end{array}$

$\Rightarrow \left|{\tilde{H}}_{{\theta }_i}\right|=n^{\left|{\theta }_i\right|}|I_{{\theta }_i}$

For which $I_{{\theta }_i}$ is the Fisher information matrix for a single data point $y_1$, and after substituting we final get for BIC:

$2\log P\left(y\mathrm{<mo>\; |\; </mo>}{M}_i\right)=2\log L\left({\hat{\theta }}_i\mathrm{<mo>\; |\; </mo>}y\right)+2\log g_i\left({\tilde{\theta }}_i\right)+\left|{\theta }_i\right|\log \left(2\pi \right)-\left|{\theta }_i\right|\log n-\log \left|I_{{\theta }_i}\right|$

4.2. AIC

4.2.1. Notation (Table 2)

Table 2. Notation 2.

4.2.2. Derivation of AIC

We can measure the quality of ${\hat{p}}_j\left(y\right)$ (as an estimate of p) by the Kullback-Leibler distance [19] :

$\begin{array}{c}K\left(p\mathrm{,}{\hat{p}}_j\right)={\displaystyle \int }p\left(y\right)\log \left(\frac{p\left(y\right)}{{\hat{p}}_j\left(y\right)}\right)\text{d}y\\ ={\displaystyle \int }p\left(y\right)\log p\left(y\right)\text{d}y-{\displaystyle \int }p\left(y\right)\log {\hat{p}}_j\left(y\right)\text{d}y\end{array}$

So, we want to minimize $K\left(p\mathrm{,}{\hat{p}}_j\right)$ over j, which is the same as maximizing

$K_j={\displaystyle \int }p\left(y\right)\log p\left(y\mathrm{<mo>\; |\; </mo>}{\hat{\theta }}_j\right)\text{d}y$

For calculating $K_j$, we can use Monte Carlo method to do an estimate

${\bar{K}}_j=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\log p\left(Y_i|{\hat{\theta }}_j\right)=\frac{{\mathcal{l}}_j\left({\hat{\theta }}_j\right)}{n}$

However, this estimate is very biased because the data are being used twice: first to get the MLE and second to estimate the integral by Monte Carlo method, and the bias is approximayely $\frac{d_j}{n}$. That means we should prove [20]

${\bar{K}}_j-\frac{d}{n}\approx K_j$

Choose ${\theta }_{j_0}$, s.t. $p\left(y|{\theta }_{j_0}\right)={\max }_{{\theta }_j\in {\Theta }_j}p\left(y|{\theta }_j\right)$, and let

$s\left(y\mathrm{,}{\theta }_j\right)=\frac{\partial \log p\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_j\right)}{\partial {\theta }_j}\mathrm{,}H\left(y\mathrm{,}{\theta }_j\right)=\frac{{\partial }^2\log p\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_j\right)}{\partial {\theta }_j^2}$

So, $s\left(y\mathrm{,}{\theta }_j\right)$ is the Jacobi martix of $\log p\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_j\right)$, and $H\left(y\mathrm{,}{\theta }_j\right)$ is the Hessian martix of $\log p\left(y\mathrm{<mo>\; |\; </mo>}{\theta }_j\right)$.

$\begin{array}{l}\Rightarrow K_j\approx {\displaystyle \int }p\left(y\right)(\log p\left(y|{\theta }_{j_0}\right)+{\left(\hat{\theta }-{\theta }_{j_0}\right)}^{\text{T}}s\left(y,{\theta }_{j_0}\right)\begin{array}{c}\\ \end{array}\\ +\frac{1}{2}{\left(\hat{\theta }-{\theta }_{j_0}\right)}^{\text{T}}H\left(y,{\theta }_{j_0}\right)\left(\hat{\theta }-{\theta }_{j_0}\right))\text{d}y\\ =K_0-\frac{1}{2n}{Z}_n^{\text{T}}J{Z}_n\end{array}$

where

$K_0={\displaystyle \int }p\left(y\right)\log p\left(y|{\theta }_{j_0}\right)\text{d}y,Z_n=\sqrt{n}\left({\hat{\theta }}_j-{\theta }_{j_0}\right),J=-E\left[H\left(y,{\theta }_{j_0}\right)\right].$

$\begin{array}{l}\Rightarrow {\bar{K}}_j\approx \frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(\mathcal{l}\left(Y_i,{\theta }_{j_0}\right)+{\left(\hat{\theta }-{\theta }_{j_0}\right)}^{\text{T}}s\left(Y_i,{\theta }_{j_0}\right)+\frac{1}{2}{\left(\hat{\theta }-{\theta }_{j_0}\right)}^{\text{T}}H\left(Y_i,{\theta }_{j_0}\right)\left(\hat{\theta }-{\theta }_{j_0}\right)\right)\\ =K_0+A_n+{\left(\hat{\theta }-{\theta }_0\right)}^{\text{T}}{S}_n-\frac{1}{2n}{Z}_n^{\text{T}}{J}_n{Z}_n\\ =K_0+A_n+\frac{Z_n^{\text{T}}{S}_n}{\sqrt{n}}-\frac{1}{2n}{Z}_n^{\text{T}}J{Z}_n\end{array}$

where

$A_n=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left(\mathcal{l}\left(Y_i,{\theta }_0\right)-K_0\right),S_n=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}s\left(Y_i,{\theta }_0\right)$

and

$J_n=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}H\left(Y_i,{\theta }_0\right)\stackrel{P}{\to }J$

${\bar{K}}_j-K_j\approx A_n+\frac{\sqrt{n}{Z}_n^{\text{T}}{S}_n}{n}\approx A_n+\frac{Z_n^{\text{T}}J{Z}_n}{n}$

From the knowledge of asymptotic distribution, we have three claims [21] :

Claim 4.1 $Z_n=\sqrt{n}\left(\hat{\theta }-{\theta }_{j_0}\right)\to N\left(0,J^{-1}V{J}^{-1}\right)$, where $V=Var\left(s\left(Y,{\theta }_{j_0}\right)\right)=J^{-1}$.

Claim 4.2 $\sqrt{n}{S}_n=\frac{\sqrt{n}}{n}{\displaystyle {\sum }_{i=1}^n}s\left(Y_i,{\theta }_{j_0}\right)\to N\left(0,V\right)$

Claim 4.3 Let $ϵ$ be a random vector with mean $\mu$ and covariance $\Sigma$, and $Q={ϵ}^{\text{T}}Aϵ$, then,

$E\left(Q\right)=trace\left(A\Sigma \right)+{\mu }^{\text{T}}A\mu$

So, with these calims above,

$\begin{array}{l}\Rightarrow E\left(\bar{K}-K\right)\approx E\left(A_n\right)+E\left(\frac{Z_n^{\text{T}}J{Z}_n}{n}\right)=0+\frac{E\left(Z_n^{\text{T}}J{Z}_n\right)}{n}\\ =\frac{trace\left(J{J}^{-1}V{J}^{-1}\right)}{n}+0^{\text{T}}J0\\ =\frac{trace\left(V{J}^{-1}\right)}{n}=\frac{trace\left(I\right)}{n}=\frac{d_j}{n}\end{array}$

$\Rightarrow {\hat{K}}_j=\frac{{\mathcal{l}}_j\left({\hat{\theta }}_j\right)}{n}-\frac{d_j}{n}={\bar{K}}_j-\frac{d_j}{n}$

So, we define

$AIC\left(j\right)=2n{\hat{K}}_j=2{\mathcal{l}}_j\left({\hat{\theta }}_j\right)-2d_j$

4.3. Example of Simple Model

Let us consider again with the example in section 3, if we take data $Y_1\mathrm{,}\cdots \mathrm{,}{Y}_n\sim N\left(\theta \mathrm{,1}\right)$, and compare it with two models, such that, $M_0\mathrm{<mo>\; :\; </mo>}N\left(\mathrm{0,1}\right)$ and $M_1\mathrm{<mo>\; :\; </mo>}N\left(\theta \mathrm{,1}\right)$. Then take the same hypothesis as in section 3.2, we test:

$H_0\mathrm{<mo>\; :\; </mo>}\theta =0\text{vs}{H}_1\mathrm{<mo>\; :\; </mo>}\theta \ne 0$

By standard normal distribution we have,

$Z=\frac{\bar{Y}-0}{\sqrt{Var\left(\bar{Y}\right)}}=\sqrt{n}\bar{Y}$

In case to avoid Type I error in our test, for $\alpha =0.05$, by Z table, we would reject $H_0$ if $\left|Z\right|>\frac{z_{\alpha }}{2}\approx 1.96$ (we take $\left|Z\right|>2$ ). Which implies if $\left|\bar{Y}\right|>\frac{2}{\sqrt{n}}$, we reflect $H_0$.

Case 1: BIC

For what we have showed in section 4.1, we proved that $BIC=2\log L\left({\hat{\theta }}_i\mathrm{<mo>\; |\; </mo>}y\right)+2\log g_i\left({\tilde{\theta }}_i\right)+\left|{\theta }_i\right|\log \left(2\pi \right)-\left|{\theta }_i\right|\log n-\log \left|I_{{\theta }_i}\right|$. However, in case to make comparison with two models, we could get away some unnecessary part, we take $BIC=\log L\left({\hat{\theta }}_i|y\right)-\frac{\left|{\theta }_i\right|}{2}\log n$. Thus,

For $H_0$,

$BIC=\log L\left(0\right)-\frac{0}{2}\log n=-\frac{n{\bar{Y}}^2}{2}-\frac{n{S}^2}{2}$

and $H_1$,

$BIC=\log L\left(\hat{\theta }\right)-\frac{1}{2}\log n=-\frac{n{S}^2}{2}-\frac{1}{2}\log n$

where $S^2={\displaystyle {\sum }_i}{\left(Y_i-\bar{Y}\right)}^2$. If we want to choose $M_1$ as a better model, then we would make $-\frac{n{\bar{Y}}^2}{2}-\frac{n{S}^2}{2}<-\frac{n{S}^2}{2}-\frac{1}{2}\log n$, in other words, $\left|\bar{Y}\right|>\sqrt{\frac{\log n}{n}}$. And BIC is an estimate of a function of the posterior probability of a model under Bayesian setup.

Case 2: AIC

And from section 4.2, for $AIC=2{\mathcal{l}}_j\left({\hat{\theta }}_j\right)-2d_j$, for which as what we have defined above $S^2={\displaystyle {\sum }_i}{\left(Y_i-\bar{Y}\right)}^2$ that $\mathcal{l}\left(\theta \right)=-\frac{n{\left(\bar{Y}-\theta \right)}^2}{2}-\frac{n{S}^2}{2}$. Further deduce $AIC={\mathcal{l}}_S-\left|S\right|$ Thus,

For $H_0$,

$AIC=\mathcal{l}\left(0\right)-0=-\frac{n{\bar{Y}}^2}{2}-\frac{n{S}^2}{2}$

and $H_1$,

$AIC=\mathcal{l}\left(\theta \right)-1=-\frac{n{S}^2}{2}-1$

If we want to choose $M_1$ as a better model at this point, we would take $-\frac{n{\bar{Y}}^2}{2}-\frac{n{S}^2}{2}<\frac{n{S}^2}{2}-1$, implies $\left|\bar{Y}\right|>\frac{\sqrt{2}}{\sqrt{n}}$. Which AIC is estimate a constant plus the relative distance between unknow likelihood function.

5. Conclusion

The question of how to choose a best model and what is a best model, it is hard to define. More precise, the controversy has existed for a long time, and no doubt it will continue longer. In this paper, we have discussed Bayes factor in hypothesis. It is obviously that Bayes factor is increasingly used in many fields of statistic research. For Bayes factor standard methods, AIC and BIC, we would consider to use for model selection. However, we also should notice that for all methods they all have their own limitation, such as the sensitivity of priors in Lindley’s paradox. Even both frequentist and Bayesian statisticians have came up with different new ideas, it is still hard to be implemented or understand by all other. Moreover, from statistic point, the method also needs to be general enough to apply. Such as for Lindley’s paradox, the partial Bayes factor in case to avoid the sensitive of priors, it takes the minimal training sample from data set to get prior and then apply with rest of the data. Partial Bayes factor at some point did deduce the influence of sensitivity of prior, but how to find the minimal training sample could also be a hard problem. Same as fractional Beyes factor, even it proves the method of choosing data for partial Bayes facto, it still has many limitations we need consider.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Vehtari, A., Gelman, A. and Gabry, J. (2017) Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC. Statistics and Computing, 27, 1413-1432.
[2]	Hartman, B.M. and Groendyke, C. (2013) Model Selection and Averaging in Financial Risk Management. North American Actuarial Journal, 17, 216-228. https://doi.org/10.1080/10920277.2013.824374
[3]	Bayarri, M.J., Berger, J.O., Forte, A. and Garca-Donato, G. (2012) Criteria for Bayesian Model Choice with Application to Variable Selection. The Annals of Statistics, 40, 1550-1577. https://doi.org/10.1214/12-AOS1013
[4]	Wasserman, L. (2000) Bayesian Model Selection and Model Averaging. Journal of Mathematical Psychology, 44, 92-107. https://doi.org/10.1006/jmps.1999.1278
[5]	Ardila, F., et al. (2008) Root Polytopes and Growth Series of Root Lattices. SIAM Journal on Discrete Mathematics, 25, 1-17.
[6]	Kass, R.E. and Raftery, A.E. (1999) Bayes Factors. Journal of the American Statistical Association, 90, 773-795.
[7]	Bayarri, M.J. and Berger, J.O. (2004) The Interplay of Bayesian and Frequentist Analysis. Statistical Science, 19, 58-80.
[8]	Spezzaferri, F., Verdinelli, I. and Zeppieri, M. (2007) Bayes Factors for Goodness of Fit Testing. Journal of Statistical Planning and Inference, 137, 43-56. https://doi.org/10.1016/j.jspi.2005.09.002
[9]	Jeffery, H. (1961) Theory of Probability.
[10]	Bhat, H.S. and Kumar, N. (2010) On the Derivation of the Bayesian Information Criterion.
[11]	Kadane, J.B. and Lazar, N.A. (2004) Methods and Criteria for Model Selection. Journal of the American Statistical Association, 99, 279-290.
[12]	Merhav, N. and Feder, M. (1998). Universal Prediction. IEEE Transactions on Information Theory, 44, 2124-2147. https://doi.org/10.1109/18.720534
[13]	Lindley, D.V. (1957) A Statistical Paradox. Biometrika, 44, 187-192. https://doi.org/10.1093/biomet/44.1-2.187
[14]	Bartlett, M.S. (1957) A Comment on D. V. Lindley’s Statistical Paradox. Biometrika, 44, 533-534. https://doi.org/10.2307/2332888
[15]	Robert, C.P. (1993) A Note on Jeffreys-Lindley Paradox. Statistica Sinica, 3, 601-608.
[16]	Villa, C. and Walker, S. (2015) On the Mathematics of the Jeffreys-Lindley Paradox. Communications in Statistics—Theory and Methods, 46, 12290-12298.
[17]	Cousins, R.D. (2014) The Jeffreys-Lindley Paradox Discovery Criteria in High Energy Physics. Synthese, 194, 395-432. https://doi.org/10.1007/s11229-015-0687-3
[18]	Shively, T. and Walker, S.G. (2013) On the Equivalence between Bayesian and Classical Hypothesis Testing.
[19]	Vardeman, S.B. (1987) Testing a Point Null Hypothesis: The Irreconcilability of p-Values and Evidence Comment. Journal of the American Statistical Association, 82, 130-131. http://www.jstor.org/stable/2289136 https://doi.org/10.2307/2289136
[20]	Penny, W.D. (2012) Comparing Dynamic Causal Models Using AIC, BIC and Free Energy. Neurolmage, 59, 319-330. https://doi.org/10.1016/j.neuroimage.2011.07.039
[21]	Zellner, A. and Siow, A. (1980) Posterior Odds Ratios for Selected Regression Hypotheses. Trabajos de Estadistica Y de Investigacion Operativa, 31, 585-603. https://doi.org/10.1007/BF02888369