=&\frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \tilde{\eta}_j)) - \frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j)) + \frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j))  - \mathbb{E}_g[H(y; f_j(\cdot; \eta^{\ast}_j))].\nonumber
\end{align}

We then use the weak law of large numbers and the Lipschitz condition (A4) to show that

\begin{align}
&\left|\frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \tilde{\eta}_j)) - \mathbb{E}_g[H(y; f_j(\cdot; \eta^{\ast}_j))]\right|\nonumber\\
=&\left|\frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \tilde{\eta}_j)) - \frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j)) + \frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j))  - \mathbb{E}_g[H(y; f_j(\cdot; \eta^{\ast}_j))]\right|\nonumber\\
\leq&\left|\frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \tilde{\eta}_j)) - \frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j))\right| + \left|\frac{1}{n}\sum_{i=1}^n H(y_i; f_j(\cdot; \eta^{\ast}_j))  - \mathbb{E}_g[H(y; f_j(\cdot; \eta^{\ast}_j))]\right|\textrm{ (tri. in)}\nonumber\\
\leq&\frac{1}{n}\sum_{i=1}^n\left| H(y_i; f_j(\cdot; \tilde{\eta}_j)) - H(y_i; f_j(\cdot; \eta^{\ast}_j))\right| + o_p(1)\textrm{ (tri. in \& WLLN)}\nonumber\\
\leq&\frac{1}{n}\sum_{i=1}^nm_H(y_i)\left|\left| \tilde{\eta}_j - \eta^{\ast}_j\right|\right|_2 + o_p(1) = o_p(1)\nonumber
\end{align}

where $\frac{1}{n}\sum_{i=1}^nm_H(y_i) \stackrel{P}{\longrightarrow} \mathbb{E}\left[m_H(y)\right] < \infty$ by the weak law of large numbers and $\left|\left|\tilde{\eta}_j - \eta^{\ast}_j\right|\right|_2 = O_p(\nicefrac{1}{\sqrt{n}})$ by Proposition <ref>.

As a result

\begin{align}
\frac{1}{n}\log \frac{\mathcal{H}_k(y)}{\mathcal{H}_l(y)} = \mathbb{E}_g[H(y; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(y; f_k(\cdot; \eta^{\ast}_k))] + o_p(1),\nonumber
\end{align}

where $\mathbb{E}_g[H(y; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(y; f_k(\cdot; \eta^{\ast}_k))] < 0$ by assumption. That is, the  ratio accrues evidence in favour of model $l$ at an exponential rate.

Case 2)

Now $\mathbb{E}_g\left[H(y; f_l(\cdot; \eta^{\ast}_{l})\right] = \mathbb{E}_g\left[H(y; f_k(\cdot; \eta^{\ast}_{k})\right]$ and we want to prove that the $\mathcal{H}$-Bayes factor convergence towards the simpler model $k$. The key is to prove that

\begin{align}
\sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) - \sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) &= O_p(1). \nonumber%\\
%\Rightarrow \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) - \log \pi(\tilde{\eta}_k) - \left\lbrace \sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \log \pi(\tilde{\eta}_k)\right\rbrace &= O_p(1)
\end{align}

Note that as model $k$ is nested in model $l$, there exists $\eta^{\ast}_{l\setminus k}$ such that

\begin{equation}
%H(y_i; f_k(\cdot; \tilde{\eta}_k)) = H(y_i; f_l(\cdot; \left(\tilde{\eta}_k, \kappa_l^{\ast}\right))).\nonumber
H(y_i; f_k(\cdot; \eta_k)) = H(y_i; f_l(\cdot; \left(\eta_k, \eta^{\ast}_{l\setminus k}\right)))\quad \forall \eta_k,\nonumber
\end{equation}

where $\eta^{\ast}_{l\setminus k}$ is the value of the hyperparameters of model $l$ that recover model $k$.
Henceforth, denote $\tilde{\eta}^{\prime}_k = \left(\tilde{\eta}_k, \eta^{\ast}_{l\setminus k}\right)$. We can rewrite our objective of interest as

\begin{align}
\sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) - \log \pi_l(\tilde{\eta}^{\prime}_k) - \left\lbrace \sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \log \pi_l(\tilde{\eta}^{\prime}_k)\right\rbrace,\nonumber
\end{align}

provided $\pi_l(\eta_l) > 0$, $\forall \eta_l$ (A1). This facilitates the Taylor expansion of $\tilde{H}_l(\tilde{\eta}^{\prime}_k) := \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}^{\prime}_k)) - \log \pi_l(\tilde{\eta}^{\prime}_k)$ about $\tilde{\eta}_l$, assuming that its first and second derivatives are finite.

\begin{align}
%	 & \tilde{H}_l(\tilde{\eta}^{\prime}_k) = \tilde{H}_l(\tilde{\eta}_l) + (\tilde{\eta}_k - \tilde{\eta}_l)\nabla_{\eta_l} \tilde{H}_l(\tilde{\eta}_l) + \frac{1}{2}(\tilde{\eta}_k - \tilde{\eta}_l)^T\nabla_{\eta_l}^2 \tilde{H}_l(\tilde{\eta}^{\dagger(n)}_l)(\tilde{\eta}_k - \tilde{\eta}_l)\nonumber
& \tilde{H}_l(\tilde{\eta}^{\prime}_k) = \tilde{H}_l(\tilde{\eta}_l) + (\tilde{\eta}_k - \tilde{\eta}_l)\tilde{\Psi}_l(\tilde{\eta}_l) + \frac{1}{2}(\tilde{\eta}_k - \tilde{\eta}_l)^T\nabla_{\eta_l}^2 A_l(\tilde{\eta}^{\dagger}_l)(\tilde{\eta}_k - \tilde{\eta}_l),\nonumber
\end{align}

for some $\eta^{\dagger}_l = \alpha\tilde{\eta}_k + (1-\alpha)\tilde{\eta}_l$ with $\alpha\in [0, 1]^{d_l}$.

Now by definition we have that $\tilde{\Psi}_l(\tilde{\eta}_l) = 0$ for all $n$ and therefore plugging this back into our target equation we get that

\begin{align}
&\sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) - \log \pi_l(\tilde{\eta}^{\prime}_k) - \left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \log \pi_l(\tilde{\eta}_l)\right\rbrace\nonumber\\
=& \log \pi_l(\tilde{\eta}_l) - \log \pi_l(\tilde{\eta}^{\prime}_k) + \frac{1}{2}(\tilde{\eta}_k - \tilde{\eta}_l)^TA_l(\tilde{\eta}^{\dagger}_l)(\tilde{\eta}_k - \tilde{\eta}_l)\nonumber\\
=& \log \pi_l(\tilde{\eta}_l) - \log \pi_l(\tilde{\eta}^{\prime}_k) +  \frac{1}{2}\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l)^T\frac{1}{n}A_l(\tilde{\eta}^{\dagger}_l)\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l).\nonumber
\end{align}

Now, $||\tilde{\eta}_j - \eta_j^{\ast}||_2 = O_p\left(1/\sqrt{n}\right)$, $j = \{k, l\}$ from Proposition <ref>, where in this nested case $\eta_l^{\ast} = \eta_k^{\ast} = \eta^{\ast}$, which further provides that $||\eta^{\dagger}_l - \eta^{\ast}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$ also. From Step 1 we know that $\log \pi_j(\tilde{\eta}_j) \stackrel{P}{\longrightarrow} \log \pi_j(\eta^{\ast}_j)$, $j = k, l$, and therefore $\log \pi_l(\tilde{\eta}_l) - \log \pi_l(\tilde{\eta}^{\prime}_k) \stackrel{P}{\longrightarrow} 0$.

Lastly, Step 2 proved that $\left|\left|\frac{1}{n}A_j(\tilde{\eta}_j) -  A^{\ast}_j(\eta^{\ast}_j)\right|\right|_2 = o_p(1)$

and therefore

\begin{align}
\sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) - \sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) &=\frac{1}{2}\sqrt{n}(O_p(\nicefrac{1}{\sqrt{n}})(A^{\ast}_l(\eta_l^{\ast}) + o_p(1))\sqrt{n}(O_p(\nicefrac{1}{\sqrt{n}})\nonumber\\
&= O_p(1) + O_p(1)o_p(1) = O_p(1)\nonumber
\end{align}

As a result, that log- accrues evidence in favour of simpler model $k$ at a $\log(n)$ rate. That is, the  accrues evidence at polynomial rate.

§.§.§ Asymptotic normality of the

While Theorem <ref> established the consistency and asymptotic normality of the frequentist parameter, the following theorem establishes the asymptotic normality of the  in  (<ref>).

This could be just for the whole $\eta = (\theta, \kappa)$ or given a plug in $\kappa$ and just for $\theta$, would come after the two stage inference

(Asymptotic normality of the  )

Let the regularity conditions of <cit.> hold. There exists a non-singular matrix $K$ such that under the  in  (<ref>) $\pi^H(\eta_k|y)$

\begin{equation}
\sqrt{n}\left(\eta_k - \hat{\eta}_k\right) \stackrel{P}{\longrightarrow} \mathcal{N}\left(0, K^{-1}\right)
\end{equation}

almost surely as $n\rightarrow\infty$ with respect to $x_{1:\infty}\sim g(\cdot)$[$\pi^H(\eta_k|y)$ is here interpreted as a random probability measures and a function of the random observations $y$] and where  $K = K\left(\theta^{\ast}_k\right)$ is given in  (<ref>)

Apply  (4) of <cit.> with $\ell(\theta_k, x) = H(x;f_k(\cdot;\theta_k))$. For the regularity conditions see <cit.> and the supplementary material of <cit.>.

§.§.§ Corollary <ref> (inconsistency for nested models)

Here we investigate whether the model selection criteria provided by <cit.> can provide the same model/loss selection consistency as Theorem <ref>, in particular focusing on nested models. We do so under the same conditions as Theorem  <ref>, where the Laplace approximation of the  were shown to be consistent, but with minor further conditions. Firstly, define

\begin{align}
I_j(\eta_j) :&= \sum_{i=1}^n\left(\nabla_{\eta_j} H(y_i; f_j(\cdot; \eta_j))\right)\left(\nabla_{\eta_j} H(y_i; f_j(\cdot; \eta_j))\right)^T\nonumber\\
I^{\ast}_j(\eta_j) :&= \mathbb{E}_g\left[\left(\nabla_{\eta_j} H(z; f_j(\cdot; \eta_j))\right)\left(\nabla_{\eta_j} H(z; f_j(\cdot; \eta_j))\right)^T\right].\nonumber
\end{align}

We then consider conditions

A5. For model $j = l, k$ there exist functions $m_I(\cdot)$ with $\mathbb{E}\left[m_I(\cdot)\right] < \infty$ so that the following Lipschitz condition holds

\begin{align}
\left|\left|I^{(1)}_j(\eta_j^{(a)}) - I^{(1)}_j(\eta_j^{(b)})\right|\right|_2 \leq m_I(z)\left|\left|\eta_j^{(a)} - \eta_j^{(b)}\right|\right|_2, \quad \forall{}\eta_j^{(a)}, \eta_j^{(b)},\nonumber
\end{align}

where the notation $I^{(1)}_j(\cdot)$ emphasises that only one observation, $z$, is involved.

A6 Models $k$ and $l$ are such that $||\tilde{\eta}_j - \eta^{\ast}_j||_2 > o_p(1/\sqrt{n})$.

A5 provides a standard Lipschitz condition for the matrices $I_j(\eta_j)$, while A6 says that the convergence of $\tilde{\eta}_j$ to $\eta^{\ast}_j$ is not faster than $1/\sqrt{n}$. While Proposition <ref> proved this convergence of $O_p(\nicefrac{1}{\sqrt{n}})$, it is rare to obtain parametric convergence of $o_p(1/\sqrt{n})$, for simplicity here we simply assume this to be the case.

Further, note that while we desired maximisation of our , the is to be minimised.

Assume Conditions A1-A4, A5 and parametric convergence according to A6. Given a sample $y \sim g$, the difference in the proposed by <cit.> of loss $k$ over loss $l$ has the following asymptotic behaviour as $n\rightarrow\infty$.

* When $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))]$, with $k$ being the simpler model

\begin{align}
\SMIC_k - \SMIC_l > o_p(1)
\end{align}

That is that when the models are equally preferable according to Fisher's divergence the will not consistently select the simpler model.

To prove this we first establish that the trace terms, $\textrm{tr}\left(\frac{1}{n}I_j(\tilde{\eta}_j)\left(\frac{1}{n}A_j(\tilde{\eta}_j)\right)^{-1}\right)$, tend to constants. Then we show that under the conditions the difference $\left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace$ is not $o_p(1)$ and therefore for any $n$ there is non-zero probability of $SMIC_k - SMIC_l < 0$ when the model $k$ was sufficient according to Fisher's divergence.

From <cit.>  (34)

\begin{align}
\SMIC_k - \SMIC_l =& \left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace\nonumber \\
& + \textrm{tr}\left(\frac{1}{n}I_k(\tilde{\eta}_k)\left(\frac{1}{n}A_k(\tilde{\eta}_k)\right)^{-1}\right) - \textrm{tr}\left(\frac{1}{n}I_l(\tilde{\eta}_l)\left(\frac{1}{n}A_l(\tilde{\eta}_l)\right)^{-1}\right).\nonumber
\end{align}

Firstly, Step 2 of Theorem <ref> proved that

\begin{align}
||\frac{1}{n}A_j(\tilde{\eta}_j) - A^{\ast}_j(\eta^{\ast})||_2 = o_p(1), \quad j = l, k.\nonumber
\end{align}

Next, we use similar arguments, and A5 to prove that $||\frac{1}{n}I_j(\tilde{\eta}_j) - I^{\ast}_j(\eta^{\ast})||_2 = o_p(1)$.

\begin{align}
\left|\left|\frac{1}{n}I_j(\tilde{\eta}_j) -  I^{\ast}_j(\eta^{\ast}_j)\right|\right|_2 &= \left|\left|\frac{1}{n}I_j(\tilde{\eta}_j) - \frac{1}{n}I_j(\eta^{\ast}_j) + \frac{1}{n}I_j(\eta^{\ast}_j) - I^{\ast}_j(\eta^{\ast}_j)\right|\right|_2\nonumber\\
&\leq \left|\left|\frac{1}{n}I_j(\tilde{\eta}_j) - \frac{1}{n}I_j(\eta^{\ast}_j)\right|\right|_2 + \left|\left|\frac{1}{n}I_j(\eta^{\ast}_j) - I^{\ast}_j(\eta^{\ast}_j)\right|\right|_2\textrm{ (tri. in)}\nonumber\\
&\leq \frac{1}{n}\sum_{i=1}^n\left|\left|I_j(\tilde{\eta}_j) - I_j(\eta^{\ast}_j)\right|\right|_2 + o_p(1)\textrm{ (tri. in \& WLLN.)}\nonumber\\
&\leq \frac{1}{n}\left|\left|\tilde{\eta}_j - \eta^{\ast}_j\right|\right|_2\sum_{i=1}^nm_I(y_i) + o_p(1)= o_p(1),\nonumber
\end{align}

where $\frac{1}{n}\sum_{i=1}^nm_I(y_i) \stackrel{P}{\longrightarrow} \mathbb{E}\left[m_I(y)\right] < \infty$ by the weak law of large numbers and $\left|\left|\tilde{\eta}_j - \eta^{\ast}_j\right|\right|_2 = O_p(\nicefrac{1}{\sqrt{n}})$ as proved in Proposition <ref>. Therefore, for constant $C^{\ast} = \textrm{tr}\left(\frac{1}{n}I_k^{\ast}(\eta^{\ast}_k)\left(\frac{1}{n}A_k^{\ast}(\eta^{\ast}_k)\right)^{-1}\right) - \textrm{tr}\left(\frac{1}{n}I_l^{\ast}(\eta^{\ast}_l)\left(\frac{1}{n}A_l^{\ast}(\eta^{\ast}_l)\right)^{-1}\right)$ we have that

\begin{align}
\SMIC_k - \SMIC_l &= \left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace + C^{\ast} + o_p(1).\nonumber
\end{align}

Lastly, Step 3, Case 2) of Theorem <ref> proved that in the nested case

\begin{align}
%\left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace = O_p(1).\nonumber
\left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace %=& \log \pi_l(\tilde{\eta}_l) - \log \pi_l(\tilde{\eta}^{\prime}_k) +  \frac{1}{2}\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l)^T\frac{1}{n}A_l(\tilde{\eta}^{\dagger}_l)\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l)\nonumber\\
= & o_p(1) + \frac{1}{2}\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l)^T\left(\frac{1}{n}A_l^{\ast}(\eta^{\ast}_l) + o_p(1)\right)\sqrt{n}(\tilde{\eta}_k - \tilde{\eta}_l).\nonumber
\end{align}

As a result unless $||\tilde{\eta}_k - \tilde{\eta}_l||_2 = o_p(1/\sqrt{n})$, which when $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))]$ will not happen by assumption (A6), then

\begin{align}
\left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace > o_p(1) \nonumber
\end{align}

as a result there exists $\epsilon^{\ast}> 0$ and $M^{\ast} > 0$ such that for all $n$

\begin{align}
P\left(\left\lbrace\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \right\rbrace > M^{\ast} \right) > \epsilon. \nonumber
\end{align}

The consequences of this are as follows, in the nested modelling scenario $\{\sum_{i=1}^n H(y_i; f_k(\cdot; \tilde{\eta}_k)) - \sum_{i=1}^n H(y_i; f_l(\cdot; \tilde{\eta}_l)) \} \geq 0$, therefore if $C^{\ast} > 0$, then model $k$ will never be preferred to model $l$. For nested models however it is often the case that $C^{\ast}$ is less than 0. In this case is is possible that $M^{\ast}$ is such that $M^{\ast} + C^{\ast} > 0$ and therefore there is still non-zero probability that model $l$ is chosen over model $k$ even as $n \rightarrow \infty$. Knowing whether this is the case is impossible without knowledge of the data generating density $g$ and therefore a practitioner can not be sure the SMIC will provide consistent model selection. Therefore, we describe it as inconsistent for this task.

Note Corollary <ref> considered the evaluated at the penalised score minimiser $\tilde{\eta}_k$ in order to simplify the extension from Theorem <ref>, however it is straightforward to see how this also holds for the unpenalised $\hat{\eta}_j$ also. A further consequence of this result is that cases where $A_j^{\ast}(\eta_j^{\ast})$ is not a finite matrix, as discussed in Section <ref>, will also affect the performance of the . For the there is no natural method to ameliorate this.

§.§.§ Theorem <ref> ($\mathcal{H}$-Bayes Factor Consistency under non-local priors)

The conditions of Theorem <ref> required that $\pi_j(\theta_j) > 0$ for all $\theta_j$ and therefore these conditions will be violated if a non-local prior <cit.> is placed on any model parameters. Here, we extend Theorem <ref> to allow for these, particularly focusing on nested models

In order to do so we consider the following extended conditions defining the form of the nested models and the non-local prior.

A7 Models $k$ and $l$ are such that $\eta_k = \{\theta_k\}$ with $\theta_k \in \Theta_k = \Theta$ (i.e. no hyperparameters $\kappa_k = \emptyset$) and $\eta_l = \{\theta_l, \kappa_l\}$ with $\kappa_l \neq \emptyset$ and $\theta_l \in \Theta_l = \Theta$. Define $\kappa_{l\setminus k}$ as the hyperparameters of model $l$ that recovers model $k$,

\begin{equation}
H(y_i; f_k(\cdot; \eta_k)) = H(y_i; f_l(\cdot; \left\lbrace\eta_k, \kappa_{l\setminus k}\right\rbrace)), \quad \textrm{for all } y_i.\nonumber
\end{equation}

A8 Models $l$ and $k$ satisfying A7 have priors

i) The prior $\pi_k(\eta_k) = \pi_k(\theta_k)$ for model $k$ satisfies A1 and the prior for model $l$ is $\pi_l(\eta_l) = \pi_l(\theta_l, \kappa_l) = \pi^{\NLP}_l(\kappa_l)\pi_l(\theta_l | \kappa_l)$ where  $\pi_l(\theta_l | \kappa_l)$ satisfy A1 also.

ii) The non-local prior on $\kappa_l$ is $\pi^{\NLP}_l(\kappa_l) = d_l(\left|\left|\kappa_l - \kappa_{l \setminus k}\right|\right|_2)\pi^{\LP}_l(\kappa_l)$,

where $\pi^{\LP}_l(\kappa_l)$ satisfies A1, $d_l(z) > 0$ for all $z \neq 0$, and there is a monotonically decreasing function $c(z): \mathbb{R}\mapsto \mathbb{R}$ such that $\lim_{z\rightarrow 0} c(z) = \infty$, $c(z/a) \leq c(z)/c(a)$, and

$\frac{c(z)}{- \log d_l(z)} = O(1)$ as $z\rightarrow 0$.

A7 defines the form of the nested models we consider. Both models share the same parameter space $\Theta$, but the simpler model $k$ has no hyperparameters. An example of A7 is where $k$ is the Gaussian model and $l$ Tukey's loss and $\kappa_{l\setminus k} = \infty$. A8 specifies that a non-local prior is placed on the additional hyperparameter of model $l$ which penalises their being too close to the value that recovered the simpler model $k$ at a rate as least as fast as $c(\cdot)$. The parameters shared by both models are given local priors.

We define the non-local prior adjusted Laplace approximate Bayes factor as

\begin{equation}
\tilde{B}_{kl}^{\mathcal{H}-\NLP} := \nicefrac{\tilde{\mathcal{H}}_k(y)}{\tilde{\mathcal{H}}^{\LP}_l(y)}\times \nicefrac{1}{d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2)},\label{Equ:LaplaceApproxHScoreNLP}
\end{equation}

where $\tilde{\mathcal{H}}_k(y)$ and $\tilde{\mathcal{H}}^{\LP}_l(y)$ are Laplace approximations to the marginal likelihoods (<ref>) whose respective priors $\pi_k(\eta_k) = \pi_k(\theta_k)$ and $\pi_l(\eta_l) = \pi_l(\theta_l | \kappa_l)\pi_l^{\LP}(\kappa_l)$ satisfy the conditions of A1 according to A8. This allows us to invoke Theorem <ref> for the asymptotic behaviour of $\nicefrac{\tilde{\mathcal{H}}_k(y)}{\tilde{\mathcal{H}}^{\LP}_l(y)}$. Some justification for this objective is provided in Remark <ref> at the end of this section.

The following Theorem proves that (<ref>)
maintains consistent model selection of the  factor as proved in Theorem <ref>, and can improve the rate at which the simpler of two nested models is selected when it is sufficient for minimising Fisher's divergence.

Assume Conditions A1-A4, models satisfying A7 and non-local priors satisfying A8.

* When $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] < 0$ then

\begin{align}
\frac{1}{n}\log \tilde{B}^{\mathcal{H}-\NLP}_{kl} = \mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] + o_p(1).
\end{align}

That is that when the more complex model $l$ decreases Fisher's divergence relative to model $k$, the non-local  factor accrues evidence in favour of model $l$, $\tilde{B}^{\mathcal{H}-\NLP}_{kl}\rightarrow 0$, at an exponential rate.

* When $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))]$, with $k$ being the simpler model, there exists a constant $M_d$ such that with arbitrarily high probability

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} > c(1/\sqrt{n})\left(\frac{d_l - d_k}{2}\frac{\log(n)}{c(1/\sqrt{n})} + o_p(1) + M_d\right).\label{Equ:logBNLP_case2}
\end{align}

That is that when the models are equally preferable according to Fisher's divergence then $\tilde{B}^{\mathcal{H}-\NLP}_{kl}\rightarrow \infty$ at a rate depending on the specification via the function $c(\cdot)$.

To prove this we decompose the log non-local  factor in the local  factor, which was dealt with by Theorem <ref>, and the non-local penalty term. Then, we consider the two cases of the theorem. When the more complicated model is true, we show that the non-local penalty term tends to a positive constant while when the simpler model is sufficient, we show that the properties of the non-local penalty function control the rate of convergence.

Firstly, from (<ref>)

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} = \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2),\label{Equ:logBNLP}
\end{align}

and by Slutsky's Theorem we can investigate both terms separately. The term $\tilde{B}^{(\mathcal{H})}_{kl}$ was constructed so that it satisfied the conditions of Theorem <ref> and therefore we can invoke these results here. Further, $\tilde{\eta}_l = \{\tilde{\theta}_l, \tilde{\kappa}_l\}$ where defined as the score minimisers penalised by the local prior, therefore satisfying the conditions of Proposition <ref> and ensuring that  $||\tilde{\eta}_j - \eta_j^{\ast}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$. Next, we consider the two cases.

Case 1)

From Theorem <ref> we have that

\begin{align}
\frac{1}{n}\log \tilde{B}^{(\mathcal{H})}_{kl} = \mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] + o_p(1).\nonumber
\end{align}

Further, for $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] < 0$ it must be the case that $\kappa^{\ast}_l \neq \kappa_{l\setminus k}$ (by A7). Therefore, $||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \nrightarrow 0$ and by the continuous mapping theorem we have that $|d(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) - d_l(||\kappa^{\ast}_l - \kappa_{l\setminus k}||_2)| = O_p(\nicefrac{1}{\sqrt{n}})$ with  $d_l(||\kappa^{\ast}_l - \kappa_{l\setminus k}||_2) > 0$ as a strictly positive constant (A8ii)). As a result $\frac{1}{n}d(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) = O(1/n))$ and

\begin{align}
\frac{1}{n}\log \tilde{B}^{\mathcal{H}-\NLP}_{kl} = \mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] + o_p(1).\nonumber
\end{align}

Case 2)

Firstly, it is convenient to rewrite

\begin{align}
%\log \tilde{B}_{kl}^{\mathcal{H}-NLP} &= \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(\tilde{\kappa}_l)\nonumber\\
%\Rightarrow\log \tilde{B}_{kl}^{\mathcal{H}-NLP}\times\frac{1}{-\log d_l(\tilde{\kappa}_l)}&= \log \tilde{B}^{(\mathcal{H})}_{kl}\ \times\frac{1}{-\log d_l(\tilde{\kappa}_l)} + 1\nonumber\\
\log \tilde{B}_{kl}^{\mathcal{H}-NLP}\times\frac{c(n)}{-\log d_l(\tilde{\kappa}_l)}&= \log \tilde{B}^{(\mathcal{H})}_{kl}\ \times\frac{c(n)}{-\log d_l(\tilde{\kappa}_l)} + c(n).\nonumber
\end{align}

From Theorem <ref> we have that

\begin{align}
\log \tilde{B}^{(\mathcal{H})}_{kl} = \frac{d_l - d_k}{2}\log(n) + O_p(1),\nonumber
\end{align}

Further, $\mathbb{E}_g[H(y; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(y; f_k(\cdot; \eta^{\ast}_k))] < 0$ requires that $\kappa^{\ast}_l = \kappa_{l\setminus k}$. Now from Proposition <ref> we have that $||\tilde{\kappa}_j - \kappa_j^{\ast}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$ and therefore using A8ii) $\frac{c(n)}{- \log d_l(\tilde{\kappa}_l)} = O_p(1)$. Therefore,

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP}\times O_p(1) &= \left(\frac{d_l - d_k}{2}\log(n) + O_p(1)\right) \times O_p(1) + c(n) \nonumber\\
&= c(n) \left[ 1 + \left(\frac{d_l - d_k}{2} \frac{\log(n)}{c(n)} O_p(1) + O_p(1/c(n))\right)  \right] \nonumber\\
&= c(n) \left[ 1 + \frac{d_l - d_k}{2} \frac{\log(n)}{c(n)} \left( O_p(1) + O_p(1/\log(n)) \right)  \right] \nonumber\\
&= c(n) \left[ 1 + \frac{d_l - d_k}{2} \frac{\log(n)}{c(n)} O_p(1) \right] \nonumber\\
\end{align}

NEW - David

Firstly, it is convenient to rewrite

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} = \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(\tilde{\kappa}_l)\nonumber \\
\color{red}
\log \tilde{B}^{(\mathcal{H})}_{kl} - c(\tilde{\kappa}_l) \frac{\log d_l(\tilde{\kappa}_l)}{c(\tilde{\kappa}_l)} \nonumber
\color{black}
\end{align}

Since $\tilde{\kappa}_l \stackrel{P}{\longrightarrow} 0$ and by A8ii),
with arbitrarily high probability as $n\rightarrow \infty$, there exists a constant $a>0$ such that

\begin{align}
\color{red}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP}=
\log \tilde{B}^{(\mathcal{H})}_{kl} + a c(\tilde{\kappa}_l)
= a c(\tilde{\kappa}_l) + \frac{d_l - d_k}{2}\log(n) + O_p(1)
\end{align}

Since $1/\tilde{\kappa}_l= O_p(\nicefrac{1}{\sqrt{n}})$, there is a $\delta>0$ such that for all fixed $\epsilon > 0$

\begin{align}
P \left( \frac{\sqrt{n}}{\tilde{\kappa}_l} > \delta \right) < \epsilon
\Rightarrow
P \left( c \left(\frac{\sqrt{n}}{\delta} \right) > c(\tilde{\kappa}_l) \right) < \epsilon
\end{align}

for all  $n \geq n_0$, and fixed $n_0$, since $c()$ is a strictly monotone function.

it follows that (proof to make sure)

\begin{align}
\color{red}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP}
= a c(\sqrt{n}) + \frac{d_l - d_k}{2}\log(n) + O_p(1)
\end{align}

NEW - Jack

In this case, $\mathbb{E}_g[H(y; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(y; f_k(\cdot; \eta^{\ast}_k))] = 0$ which by A7 requires that $\kappa^{\ast}_l = \kappa_{l\setminus k}$. Therefore, $||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$, which by the continuous mapping theorem implies $d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) \rightarrow0$. Here we establish the asymptotic rate of this convergence.

By (<ref>), in order to establish (<ref>) we require that

\begin{align}
\forall \epsilon,\quad \exists M_d, N_d > 0 \textrm{ such that } P\left(- \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > M_dc(1/\sqrt{n})\right) > 1 - \epsilon,\quad \forall n > N_d  \nonumber
\end{align}

From A8ii), $\frac{c(z)}{- \log d_l(z)} = O(1)$ as $z\rightarrow 0$ which guarantees the existence of $M_{\delta_{\kappa}}$ and $\delta_{\kappa}>0$ such that

\begin{equation}
\frac{- \log d_l(z)}{c(z)} > M_{\delta_{\kappa}} \textrm{ for all } z \textrm{ such that } |z| < \delta_{\kappa}. \nonumber
\end{equation}

As a result, by the law of total probability

\begin{align}
&P(- \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > M_dc(1/\sqrt{n})) \nonumber\\
=& P(-\log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > M_dc(1/\sqrt{n}) \mid ||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2\leq \delta_{\kappa})P(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \leq \delta_{\kappa}) \nonumber\\
&+ P(-\log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > M_dc(1/\sqrt{n})\mid ||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 > \delta_{\kappa})P(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 > \delta_{\kappa}) \nonumber\\
\geq&  P(c(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > \frac{M_d}{M_{\delta_{\kappa}}}c(1/\sqrt{n}) \mid ||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2\leq \delta_{\kappa})P(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \leq \delta_{\kappa}). \nonumber
\end{align}

Given that $||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$ from Proposition <ref> and $\delta_{\kappa}$ is a constant, the term $P(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \leq \delta_{\kappa})$ is arbitrarily close to 1 as $n\rightarrow\infty$. Regarding the other terms, by the decreasing monotonicity of $c(z)$ (A8ii)) the conditioning events is such that

\begin{align}
\left\lbrace||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \leq \delta_{\kappa}\right\rbrace = \left\lbrace c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2\right) \geq c(\delta_{\kappa})\right\rbrace,\nonumber
\end{align}

and therefore for $\frac{M_d}{M_{\delta_{\kappa}}}c(1/\sqrt{n}) \geq c(\delta_{\kappa})$ we have that

\begin{align}
P\left(c(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > \frac{M_d}{M_{\delta_{\kappa}}}c(1/\sqrt{n}) \mid c(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > c(\delta_{\kappa})\right) \geq P\left(c(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > \frac{M_d}{M_{\delta_{\kappa}}}c(1/\sqrt{n})\right).\nonumber
\end{align}

Since $||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 = O_p(\nicefrac{1}{\sqrt{n}})$, by definition for every $\epsilon' > 0$ there exists $M_{\epsilon'} >0$ and $N_{\epsilon'} > 0$ such that

\begin{align}
P(\sqrt{n}||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2  < M_{\epsilon'}) > 1 - \epsilon',\quad \forall n > N_{\epsilon}. \nonumber
\end{align}

Further, since $c(z)$ is strictly decreasing and $c(z/a) \leq c(z)/c(a)$ by Assumption A8ii), we have that

\begin{align}
P(c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \right) / c(1/\sqrt{n}) > c(M_{\epsilon'}))
\geq P(c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 /(1/\sqrt{n})\right) > c(M_{\epsilon'})) > 1 - \epsilon'\quad  \forall n > N_{\epsilon'}.\nonumber
%\Rightarrow& P(c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \right) > c(\delta_{\epsilon})c(1/\sqrt{n})) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber\\
%\Rightarrow&  P(c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \right) / c(1/\sqrt{n}) > c(M_{\epsilon})) > \epsilon\quad  \forall n > N_{\epsilon}\nonumber
\end{align}

Consider the particular choice $M_d = c(M_{\epsilon'})M_{\delta_{\kappa}}$ so that $\frac{M_d}{M_{\delta_{\kappa}}} = c(M_{\epsilon'})$ giving that

\begin{align}
P(c\left(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2 \right) > \frac{M_d}{M_{\delta_{\kappa}}} c(1/\sqrt{n})) > 1 - \epsilon'\quad  \forall n > N_{\epsilon'}\nonumber
\end{align}

Therefore, taking $(1 - \epsilon') > (1 - \epsilon)$

\begin{align}
P(- \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > M_dc(1/\sqrt{n})) > (1 - \epsilon')P(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2) > (1-\epsilon)\nonumber
\end{align}

for large enough $n$. Further, invoking Theorem <ref> we have that with arbitrarily high probability there exists constant $M_d$ such that for sufficiently large $n$

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-NLP} &= \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2)\nonumber\\
&> \log \tilde{B}^{(\mathcal{H})}_{kl} + M_d \times c(1/\sqrt{n})\nonumber\\
&= \frac{d_l - d_k}{2}\log(n) + O_p(1) + M_d \times c(1/\sqrt{n})\nonumber\\
&= c(1/\sqrt{n})\left(\frac{d_l - d_k}{2}\frac{\log(n)}{c(1/\sqrt{n})} + \frac{O_p(1)}{c(1/\sqrt{n})} + M_d\right).\nonumber
\end{align}

It is convenient to rewrite

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} &= \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)\nonumber \\
&= \log \tilde{B}^{(\mathcal{H})}_{kl} - c(|\tilde{\kappa}_l - \kappa_{l\setminus k}|) \frac{\log d_l(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)}{c(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)} \nonumber
\end{align}

From A8ii), $\frac{c(\kappa^{\prime})}{- \log d_l(\kappa^{\prime})} = O(1)$ as $\kappa^{\prime}\rightarrow 0$ guarantees the existence of $M_{\delta_{\kappa}}$ and $\delta_{\kappa}>0$ such that

\begin{equation}
\frac{- \log d_l(\kappa^{\prime})}{c(\kappa^{\prime})} > M_{\delta_{\kappa}} \textrm{ for all } \kappa_l \textrm{ s.t. } |\kappa^{\prime}| < \delta_{\kappa}. \nonumber
\end{equation}

Now, from Proposition <ref> $||\tilde{\kappa}_l - \kappa_{l\setminus k}|| = O_p(\nicefrac{1}{\sqrt{n}})$ and therefore for all $\epsilon$ there exists $N_{\epsilon}$ and $\delta_{\epsilon}$ such that

\begin{align}
P(\sqrt{n}|\tilde{\kappa}_l - \kappa_{l\setminus k}| < \delta_{\epsilon}) > \epsilon,\quad \forall n > N_{\epsilon} \nonumber
\end{align}

As a result, with arbitrarily high probability
$|\tilde{\kappa}_l - \kappa_{l\setminus k}| < \frac{\delta_{\epsilon}}{\sqrt{n}} < \delta_{\epsilon}$ and therefore that

\begin{align}
\frac{- \log d_l(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)}{c(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)} > M_{\delta_{\epsilon}}\label{Equ:d_over_c}
\end{align}

Next, we use Proposition <ref> again to investigate $c(|\tilde{\kappa}_l - \kappa_{l\setminus k}|)$. By the monotonically decreasing nature of $c(\kappa^{\prime})$ (A8ii)) we have that

\begin{align}
P(\sqrt{n}|\tilde{\kappa}_l - \kappa_{l\setminus k}|\geq \delta_{\epsilon}) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber\\
\Rightarrow P(c\left(\sqrt{n}|\tilde{\kappa}_l - \kappa_{l\setminus k}|\right) \leq c(\delta_{\epsilon})) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber
\end{align}

and further by the fact that $c(\kappa^{\prime}/a) = c(\kappa^{\prime})/c(a)$

\begin{align}
P(c\left(|\tilde{\kappa}_l - \kappa_{l\setminus k}|\right)/c(1/\sqrt{n}) \leq c(\delta_{\epsilon})) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber\\
%\Rightarrow P(c\left(|\tilde{\kappa}_l - \kappa_{l\setminus k}|\right) \leq c(\delta_{\epsilon})c(1/\sqrt{n})) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber\\
\Rightarrow P(c\left(|\tilde{\kappa}_l - \kappa_{l\setminus k}|\right) > c(\delta_{\epsilon})c(1/\sqrt{n})) > \epsilon\quad  \forall n > N_{\epsilon}\nonumber
\end{align}

Therefore with arbitrarily high probability we have that $c\left(|\tilde{\kappa}_l - \kappa_{l\setminus k}|\right) > c(\delta_{\epsilon})c(1/\sqrt{n}))$ where $c(1/\sqrt{n}) \rightarrow\infty$ as $n\rightarrow\infty$ by A7(ii). Combining this with (<ref>) Am I okay to bound these two probabilistic bounds in this way. I think they key to this is that we have a non-probabilistic bound, that is activated by the probabilistic bound, there exists a constant $M_c$ such that with arbitrarily high probability

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-NLP} &> \log \tilde{B}^{(\mathcal{H})}_{kl} + M_c \times c(1/\sqrt{n}).\nonumber
\end{align}

Lastly, from Theorem <ref> we have that

\begin{align}
\log \tilde{B}^{(\mathcal{H})}_{kl} = \frac{d_l - d_k}{2}\log(n) + O_p(1),\nonumber
\end{align}

and as a result

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-NLP} &> c(1/\sqrt{n})\left(\frac{d_l - d_k}{2}\frac{\log(n)}{c(1/\sqrt{n})} + \frac{O_p(1)}{c(1/\sqrt{n})} + M_c\right)\nonumber
\end{align}

We follow Theorem <ref> with the special case corollary considering the inverse-gamma non-local prior applied to $\nu_2 = \frac{1}{\kappa_2^2}$ where $\kappa_2$ was the cut-off parameter of Tukey's loss (Section <ref>). First, we generalise that prior set-up from Section <ref> as the following condition.

A9 Under the model setup of A7

i) Let $\kappa_{l} = \{\nu\}\in \mathbb{R}$ be a univariate hyperparameter with $\kappa_{l\setminus k} = \nu_0$ and $\nu \geq \nu_0$.

ii) $\pi_l^{\LP}(\nu) = 2\mbox{I}_{\nu \geq \nu_0}\mathcal{N}(\nu - \nu_0; 0, s_0^2)$, a half-Gaussian distribution with scale $s_0$.

iii) $\pi_l^{\NLP}(\nu) = \mathcal{IG}(\nu - \nu_0; a_0, b_0)$, an inverse-gamma distribution with shape $a_0$ and scale $b_0$.

We note that the corresponding local prior to an inverse-gamma non-local prior would normally be a gamma distribution. However, this does not have finite second derivative as $\nu - \nu_0\rightarrow 0$ and would therefore violate the conditions of Theorem <ref>.

Assume Conditions A1-A4, models satisfying A7 and non-local priors satisfying A8(i) and A9.

* When $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] < 0$ then

\begin{align}
\frac{1}{n}\log \tilde{B}^{\mathcal{H}-\NLP}_{kl} = \mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] - \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))] + o_p(1).
\end{align}

That is that when the more complex model $l$ decreases Fisher's divergence relative to model $k$, the non-local  factor accrues evidence in favour of model $l$, $\tilde{B}^{\mathcal{H}-\NLP}_{kl}\rightarrow 0$, at an exponential rate.

* When $\mathbb{E}_g[H(z; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(z; f_k(\cdot; \eta^{\ast}_k))]$, with $k$ being the simpler model, there exists a constant $M_d$ such that with arbitrarily high probability

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} > \sqrt{n}(o_p(1) + M_d).
\end{align}

That is that when the models are equally preferable according to Fisher's divergence then $\log \tilde{B}^{\mathcal{H}-\NLP}_{kl}\rightarrow \infty$ at a rate at least as fast as $\sqrt{n}$ as $n \rightarrow\infty$.

this proof doesn't actually use anything from Theorem <ref>, it was useful to prove from first principles, however in the paper we can just put the simplified version below

To prove this we decompose the log non-local  factor in the local  factor, which is dealt with by Theorem <ref> and the non-local penalty term provided by the inverse-gamma prior. Then, we consider the two cases of the theorem. When the more complicated model is true, the we show the non-local penalty term tends to a positive constant while, when the simpler model is sufficient, we show that the properties of the inverse-gamma non-local prior control the rate of convergence.

Firstly, the non-local penalty function $d_l(|\nu - \nu_0|) = d_l(\nu - \nu_0) = \frac{\pi^{\NLP}_l(\nu)}{\pi^{\LP}_l(\nu)}$ given by the prior specification of A9 is

\begin{align}
d_l(|\nu - \nu_0|) &= \frac{\pi_l^{\NLP}(\nu)}{\pi_l^{\LP}(\nu)}\nonumber\\
%&= \frac{\frac{b_0^{a_0}}{\Gamma(a_0)}\left(\frac{1}{\nu_2}\right)^{a_0 + 1}\exp\left(-\frac{b_0}{\nu_2}\right)}{2\frac{1}{\sqrt{2\pi}\sqrt{s_0^2}}\exp\left(- \frac{\nu_2^2}{s_0^2}\right)}\nonumber\\
&\propto \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)}\left(\frac{1}{\nu - \nu_0}\right)^{a_0 + 1}\exp\left(\frac{(\nu - \nu_0)^2}{2s_0^2} - \frac{b_0}{\nu - \nu_0}\right)\nonumber
\end{align}

and therefore

\begin{align}
- \log d_l(|\nu - \nu_0|) &= - \log \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)} + (a_0 + 1)\log \left(\nu - \nu_0\right) - \frac{\left(\nu - \nu_0\right)^2}{2s_0^2} + \frac{b_0}{\left(\nu - \nu_0\right)}.\nonumber
\end{align}

Then there exists $\delta_{\nu}, M_{\delta_{\nu}} >0$ such that for all $\nu$ s.t. $|\nu - \nu_0| < \delta_{\nu}$

\begin{align}
- \log d_l(|\nu - \nu_0|) > M_{\delta_{\nu}}\frac{1}{\left|\nu - \nu_0\right|}. \label{Equ:d_bound}
\end{align}

Now we consider the two cases.

Case 1)

Follows directly from the proof of Theorem <ref>, noting that $d_l(z) > 0$ for all $0 < z < \infty$.

Case 2)

In this case $\mathbb{E}_g[H(y; f_l(\cdot; \eta^{\ast}_l))] = \mathbb{E}_g[H(y; f_k(\cdot; \eta^{\ast}_k))] = 0$ which by A7 and A9 requires that $\nu^{\ast}_l = \nu_0$. Now from Proposition <ref> we have that $|\tilde{\nu} - \nu_{0}| = O_p(\nicefrac{1}{\sqrt{n}})$, which by the continuous mapping theorem implies $- \log d_l(|\tilde{\nu} - \nu_{0}|) \rightarrow 0$. Here we establish the asymptotic rate of this convergence.

Firstly, from (<ref>)

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} = \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2)\nonumber
\end{align}

in order to establish the result we require that

\begin{align}
\forall \epsilon > 0 ,\quad \exists M_d, N_d > 0 \textrm{ such that } P( -\log d_l(|\tilde{\nu} - \nu_{0}|) > M_d\sqrt{n}) > 1 - \epsilon,\quad \forall n > N_d.  \nonumber
\end{align}

Now from (<ref>), there exists $\delta_{\nu} >0 $ and $M_{\delta_{\nu}} > 0$ such that

\begin{align}
&P(-\log d_l(|\tilde{\nu} - \nu_{0}|) > M_d\sqrt{n}) \nonumber\\
=& P(-\log d_l(|\tilde{\nu} - \nu_{0}|) > M_d\sqrt{n} \mid |\tilde{\nu} - \nu_{0}| \leq \delta_{\nu})P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu}) \nonumber\\
&+ P(-\log d_l(|\tilde{\nu} - \nu_{0}|) > M_d\sqrt{n} \mid |\tilde{\nu} - \nu_{0}| > \delta_{\nu})P(|\tilde{\nu} - \nu_{0}| > \delta_{\nu}) \nonumber\\
>& P(M_{\delta_{\nu}}\frac{1}{|\tilde{\nu} - \nu_{0}|} > M_d\sqrt{n} \mid |\tilde{\nu} - \nu_{0}| \leq \delta_{\nu})P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu}) \nonumber\\
%    =& P(\frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n} \mid |\tilde{\nu} - \nu_{0}| \leq \delta_{\nu})P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu}) \nonumber\\
=& P(\frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n} \mid \frac{1}{|\tilde{\nu} - \nu_{0}|} \geq \frac{1}{\delta_{\nu}})P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu}) \nonumber
\end{align}

Given that $\left|\tilde{\nu} - \nu_0\right| = O_p(\nicefrac{1}{\sqrt{n}})$ from Proposition <ref> and $\delta_\nu$ is a constant, the term $P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu})$ is arbitrarily close to 1 as $n \rightarrow \infty$.
Regarding the other term, if $\frac{1}{\delta_{\nu}} \leq \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n}$
then we can say that

\begin{align}
P(\frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n} \mid \frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{1}{\delta_{\nu}}) > P(\frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n}).\nonumber
\end{align}

Since $\left|\tilde{\nu} - \nu_0\right| = O_p(\nicefrac{1}{\sqrt{n}})$, by definition for every $\epsilon'$ there exists $M_{\epsilon'}>0$ and $N_{\epsilon'}>0$ such that

\begin{align}
%P(\sqrt{n}\left|\tilde{\nu} - \nu_0\right|\geq M_{\epsilon}) \leq \epsilon\quad  \forall n > N_{\epsilon}\nonumber\\
%\Rightarrow P\left(\frac{1}{\sqrt{n}}\frac{1}{\left|\tilde{\nu} - \nu_0\right|} \leq \frac{1}{M_{\epsilon}}\right) \leq \epsilon\quad  \forall n > N_{\epsilon}\\
P\left(\frac{1}{\left|\tilde{\nu} - \nu_0\right|} > \frac{\sqrt{n}}{M_{\epsilon'}}\right) > 1 - \epsilon' \quad  \forall n > N_{\epsilon'}.\nonumber
\end{align}

Consider the particular choice  $M_d = \frac{M_{\delta_{\nu}}}{M_{\epsilon'}}$, so that

$M_{\epsilon'} = \frac{M_{\delta_{\nu}}}{M_d}$, giving that

\begin{align}
P(\frac{1}{|\tilde{\nu} - \nu_{0}|} > \frac{M_d}{M_{\delta_{\nu}}}\sqrt{n}) > 1 - \epsilon' \quad  \forall n > N_{\epsilon'}.\nonumber
\end{align}

Therefore, taking $ 1 - \epsilon' > 1 - \epsilon$

\begin{align}
P(-\log d_l(|\tilde{\nu} - \nu_{0}|) > M_d\sqrt{n}) > (1 - \epsilon')P(|\tilde{\nu} - \nu_{0}| \leq \delta_{\nu}) > (1 - \epsilon) \nonumber%, \quad \forall n > \max\{N_{\epsilon}, \left(\frac{M_{\delta_{\nu}}}{\delta_{\nu}M_d}\right)^2\}\nonumber.
\end{align}

for large enough $n$.

Further invoking Theorem <ref> we have that with arbitrarily high probability there exists $M_d$ such that for sufficiently large $n$

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-NLP} &= \log \tilde{B}^{(\mathcal{H})}_{kl} - \log d_l(|\tilde{\nu} - \nu_{0}|)\nonumber\\
&> \frac{d_l - d_k}{2}\log(n) + O_p(1) + M_d\sqrt{n}\nonumber\\
&= \sqrt{n}\left(\frac{d_l - d_k}{2}\frac{\log(n)}{\sqrt{n}} + \frac{O_p(1)}{\sqrt{n}} + M_d\right)\nonumber\\
&= \sqrt{n}(o_p(1) + M_d)\nonumber
\end{align}

Simplified proof directly calling Theorem <ref>

To prove this we derive the non-local penalty function associated with the prior setup of A9 and calculate its bounding function $c(\cdot)$ associated with A8ii). This allows us to invoke Theorem <ref> to provide the desired asymptotic behaviour.

Firstly, from A9, $\nu \geq \nu_0$ and therefore $d_l(|\nu - \nu_0|) = d_l(\nu - \nu_0)$. Then, the non-local penalty function $d_l(\nu - \nu_0) = \frac{\pi^{\NLP}_l(\nu)}{\pi^{\LP}_l(\nu)}$ given by the prior specification of A9 is

\begin{align}
d_l(\nu - \nu_0) &= \frac{\pi_l^{\NLP}(\nu)}{\pi_l^{\LP}(\nu)}\nonumber\\
%&= \frac{\frac{b_0^{a_0}}{\Gamma(a_0)}\left(\frac{1}{\nu_2}\right)^{a_0 + 1}\exp\left(-\frac{b_0}{\nu_2}\right)}{2\frac{1}{\sqrt{2\pi}\sqrt{s_0^2}}\exp\left(- \frac{\nu_2^2}{s_0^2}\right)}\nonumber\\
&\propto \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)}\left(\frac{1}{\nu - \nu_0}\right)^{a_0 + 1}\exp\left(\frac{(\nu - \nu_0)^2}{2s_0^2} - \frac{b_0}{\nu - \nu_0}\right),\nonumber
\end{align}

and therefore

\begin{align}
- \log d_l(\nu - \nu_0) &= - \log \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)} + (a_0 + 1)\log \left(\nu - \nu_0\right) - \frac{\left(\nu - \nu_0\right)^2}{2s_0^2} + \frac{b_0}{\left(\nu - \nu_0\right)}.\nonumber
\end{align}

Therefore, there exists $\delta_{\nu}>0$ and $M_{\delta_{\nu}} >0$ such that for all $\nu$ s.t. $|\nu - \nu_0| < \delta_{\nu}$

\begin{align}
- \log d_l(\nu - \nu_0) > M_{\delta_{\nu}}\frac{1}{\nu - \nu_0}, \nonumber%\label{Equ:d_bound}
\end{align}

and as a result, $c(z) = \frac{1}{z}$. Now we consider the two cases.

Case 1)

Follows directly from the proof of Theorem <ref>, noting that $d_l(z) > 0$ for all $0 < z < \infty$.

Case 2)

Using Theorem <ref> with $c(z) = \frac{1}{z}$ provides that with arbitrarily high probability there exists $M_d$ such that for sufficiently large $n$

\begin{align}
\log \tilde{B}_{kl}^{\mathcal{H}-\NLP} &> c(1/\sqrt{n})\left(\frac{d_l - d_k}{2}\frac{\log(n)}{c(1/\sqrt{n})} + \frac{O_p(1)}{c(1/\sqrt{n})} + M_d\right)\nonumber\\
&= \sqrt{n}\left(\frac{d_l - d_k}{2}\frac{\log(n)}{\sqrt{n}} + \frac{O_p(1)}{\sqrt{n}} + M_d\right)\nonumber\\
&= \sqrt{n}(o_p(1) + M_d),\nonumber
\end{align}

as required.

Lastly, we provide some justification for the objective function (<ref>) to conduct model selection using non-local priors.
Using the arguments of <cit.>, the integrated  (<ref>) associated with the non-local prior specification A8 can be rewritten as

\begin{align}
\mathcal{H}^{\NLP}_l(y) &= \int\pi^{\NLP}_l(\kappa_l)\pi_l(\theta_l|\kappa_l) \exp \left\{ -\sum_{i=1}^nH(y_i; f_l(\cdot;\theta_l, \kappa_l)) \right\} d\theta_ld\kappa_l\nonumber\\
&= \mathcal{H}^{\LP}_l(y)\int d_l(||\kappa_l - \kappa_{l\setminus k}||_2) \frac{\pi^{\LP}_l(\kappa_l)\pi_l(\theta_l| \kappa_l) \exp \left\{ -\sum_{i=1}^nH(y_i; f_l(\cdot;\theta_l, \kappa_l)) \right\} }{\mathcal{H}^{\LP}_l(y)}d\theta_ld\kappa_l\nonumber\\
&= \mathcal{H}^{\LP}_l(y)\int d_l(||\kappa_l - \kappa_{l\setminus k}||_2) \pi^{\LP}_l(\theta_l, \kappa_l|y)d\theta_ld\kappa_l\nonumber
\end{align}

\begin{align}
\pi^{\LP}_l(\theta_l, \kappa_l|y) &\propto \pi_l^{\LP}(\kappa_l)\pi_l(\theta_l| \kappa_l) \exp \left\{ -\sum_{i=1}^nH(y_i; f_l(\cdot;\theta_j, \kappa_j)) \right\},\nonumber\\
\mathcal{H}^{\LP}_l(y) &= \int\pi_l^{\LP}(\kappa_l)\pi_l(\theta_l| \kappa_l) \exp \left\{ -\sum_{i=1}^nH(y_i; f_l(\cdot;\theta_l, \kappa_l)) \right\} d\theta_ld\kappa_l,\nonumber
\end{align}

As a result we can write the  factor of nested models satisfying A7 and non-local priors according to A8

\begin{align}
B_{kl}^{\mathcal{H}-\NLP} = \frac{\mathcal{H}_k(y)}{\mathcal{H}^{\NLP}_l(y)} = \frac{\mathcal{H}_k(y)}{\mathcal{H}^{\LP}_l(y)}\frac{1}{\int d_l(||\kappa_l - \kappa_{l\setminus k}||_2) \pi^{\LP}_l(\theta_l, \kappa_l|y)d\theta_ld\kappa_l}.\nonumber
\end{align}

Now, the terms involved in the construction of $\mathcal{H}_k(y)$ and $\mathcal{H}^{\LP}_l(y)$ satisfy the conditions of Theorem <ref> (by A8i)), and therefore the behaviour of the Laplace approximations (<ref>) to the  factor,  $\nicefrac{\tilde{\mathcal{H}}_k(y)}{\tilde{\mathcal{H}}^{\LP}_l(y)}$, is detailed in Theorem <ref>. We further consider a somewhat looser, but convenient, approximation of $\int d_l(||\kappa_l - \kappa_{l\setminus k}||_2) \pi^{\LP}_l(\theta_l, \kappa_l|y)d\theta_ld\kappa_l$ as $d_l(||\tilde{\kappa}_l - \kappa_{l\setminus k}||_2)$ where $\tilde{\kappa}_l \in \tilde{\eta}_l$ are the parameters maximising the  $\pi^{\LP}_l(\theta_l, \kappa_l|y)$. As a result we consider (<ref>) to approximate the  factor under non-local priors

We highlight that calculating the approximation $\tilde{\mathcal{H}}^{\NLP}_l(y)$ uses the maximum a posteriori estimates and evaluates the observed Hessian matrices according to the  under the local prior, before evaluating the prior density as it appears in the Laplace approximation (<ref>) at the non-local prior, as shown in  (<ref>). This is convenient as it allows for a direct extension of Theorem <ref>.

§.§ Derivation of scores

Section <ref> provides the score for the Gaussian log-likelihood and for Tukey's loss.

§.§.§ The score of the Gaussian Model and Tukey's Loss

Here we provide the derivations of the score applied to the Gaussian model (<ref>) and Tukey's loss (<ref>). Firstly, for the Gaussian model $f_1(y_i; x_i, \theta_1, \kappa_1) = \mathcal{N}(y_i;x_i^T \beta,\sigma^2)$ we have that

\begin{align}
\log f_1(y_i; x_i, \theta_1, \kappa_1) &= -\frac{1}{2}\log (2\pi) - \frac{1}{2}\log \sigma^2 -\frac{(y_i - x_i^T \beta)^2}{2\sigma^2}\nonumber\\
\frac{\partial }{\partial y_i}\log f_1(y_i;x_i, \theta_1, \kappa_1) &= -\frac{(y_i - x_i^T \beta)}{\sigma^2}\nonumber\\
\frac{\partial^2 }{\partial y_i^2}\log f_1(y_i;x_i, \theta_1, \kappa_1) &= -\frac{1}{\sigma^2},\nonumber
\end{align}

which results in score

\begin{align}
H_{1}(y_i; f(\cdot;x_i, \theta_1)) &= -\frac{2}{\sigma^2} + \frac{(y_i - x_i^T \beta)^2}{\sigma^4}.\nonumber
\end{align}

For Tukey's loss, $\log f_2(y_i; x_i, \theta_2, \kappa_2) = -\ell_2(y_i; x_i, \theta_2, \kappa_2)$ given in  (<ref>) and

\begin{align}
\frac{\partial}{\partial y_i}\ell_2(y_i; x_i, \theta_2, \kappa_2) &= \begin{cases}
\frac{(y_i - x_i^T \beta)}{\sigma^2}-\frac{2(y_i - x_i^T \beta)^3}{\kappa_2^2\sigma^4}+\frac{(y_i - x_i^T \beta)^5}{\kappa_2^4\sigma^6} &\textrm{ if } |y_i - x_i^T \beta|\leq \kappa_2\sigma\\
0 &\textrm{ otherwise}
\end{cases}\nonumber\\
\frac{\partial^2}{\partial y_i^2}\ell_2(y_i; x_i, \theta_2, \kappa_2) &= \begin{cases}
\frac{1}{\sigma^2} - \frac{6(y_i - x_i^T \beta)^2}{\kappa_2^2\sigma^4} + \frac{5(y_i - x_i^T \beta)^4}{\kappa_2^4\sigma^6} &\textrm{ if } |y_i - x_i^T \beta|\leq \kappa_2\sigma\\
0 &\textrm{ otherwise}
\end{cases},\nonumber
\end{align}

which results in score

\begin{align}
H_{2}(y_i; f(\cdot;x_i, \theta_2, \kappa_2)) &=  \mbox{I}(|y_i - x_i^T \beta|\leq \kappa_2\sigma)\left\lbrace \left(\frac{(y_i - x_i^T \beta)}{\sigma^2}-\frac{2(y_i - x_i^T \beta)^3}{\kappa_2^2\sigma^4}+\frac{(y_i - x_i^T \beta)^5}{\kappa_2^4\sigma^6}\right)^2 \right.\nonumber\\
&\qquad\qquad\qquad\left. - 2\left(\frac{1}{\sigma^2}-\frac{6(y_i - x_i^T \beta)^2}{\kappa_2^2\sigma^4}+\frac{5(y_i - x_i^T \beta)^4}{\kappa_2^4\sigma^6}\right)\right\rbrace.\nonumber%\label{Equ:HscoreTukeysLoss}
\end{align}

§.§.§ The score of kernel density estimation

We derive the score for the kernel density estimation examples implemented in Section <ref>. Combining the kernel density estimate improper density (<ref>) with the power $w$  (<ref>) results in in-sample

\begin{align}
\hat{g}_{w, h}(y_i) =& \left(\frac{1}{n\sqrt{2\pi} h}\sum_{j=1}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)^{w},\nonumber
\end{align}

or equivalently defining the kernel density loss function to be the log-density

\begin{align}
\log \hat{g}_{w, h}(y_i) =& w\log\left(\frac{1}{n\sqrt{2\pi} h}\sum_{j=1}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)\nonumber\\
%=& w\log\left(\frac{1}{n\sqrt{2\pi\sigma^2}}exp\left(0\right) + \frac{1}{n\sqrt{2\pi\sigma^2}}\sum_{j\neq i}^n\exp\left(-\frac{(x_i - x_j)^2}{2\sigma^2}\right)\right)\\
=& w\log\left(\frac{1}{n\sqrt{2\pi} h} + \frac{1}{n\sqrt{2\pi} h}\sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right).\nonumber
\end{align}

The score is then composed of the first an second derivatives of the log-density, given by

\begin{align}
\frac{\partial}{\partial y_i} \log \hat{g}_{w,h}(y_i) =& w\frac{-\sum_{j\neq i}^n\frac{(y_i - y_j)}{n\sqrt{2\pi}h^{3}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)}{\left(\frac{1}{n\sqrt{2\pi} h} + \frac{1}{n\sqrt{2\pi h}}\sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)}\nonumber\\
\frac{\partial^2}{\partial y_i^2} \log \hat{g}_{w, h}(y_i) =& w\frac{\sum_{j\neq i}^n\frac{(y_i - y_j)^2}{n\sqrt{2\pi}h^{5}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right) - \sum_{j\neq i}^n\frac{1}{n\sqrt{2\pi}h^{3}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)}{\left(\frac{1}{n\sqrt{2\pi} h} + \frac{1}{n\sqrt{2\pi} h}\sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)}\nonumber\\
&\quad - w\frac{\left(\sum_{j\neq i}^n\frac{(y_i - y_j)}{n\sqrt{2\pi}h^{3}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)^2}{\left(\frac{1}{n\sqrt{2\pi} h} + \frac{1}{n\sqrt{2\pi} h}\sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)^2}.\nonumber
\end{align}

As a result, the score of the kernel density estimate is

\begin{align}
H(y_i, \log \hat{g}_{w,h}(\cdot))
&= 2w\frac{\sum_{j\neq i}^n\frac{(y_i - y_j)^2}{h^{4}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right) - \sum_{j\neq i}^n\frac{1}{h^{2}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)}{\left(1 + \sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)}\nonumber\\
&\quad -(2w - w^2)\frac{\left(\sum_{j\neq i}^n\frac{(y_i - y_j)}{h^{2}}\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)^2}{\left(1 + \sum_{j\neq i}^n\exp\left(-\frac{(y_i - y_j)^2}{2h^2}\right)\right)^2}.
\end{align}

§.§ The Tsallis score

The robust regression example introduced in Section <ref> focused on Tukey's loss as this provided a traditional alternative to the Gaussian likelihood that had one parameter, was nested in the Gaussian loss and corresponded to a improper model. A loss function satisfying similar properties is the Tsallis score <cit.> (often referred to as the density power divergence <cit.> or $\beta$-divergence) which applied to proper probability model $f(y_i; \theta)$ is given by

\begin{align}
\ell_{\beta}(y_i, \theta) := - \frac{1}{\beta} f(y_i; \theta)^{\beta} + \frac{1}{\beta + 1}\int  f(z; \theta)^{\beta+1}dz.
\end{align}

As $\beta \rightarrow 0$ then $\ell_{\beta}(y_i, \theta) \rightarrow -\log f(y_i, \theta)$ and therefore $f(y_i; \theta)$ is nested within $\exp\left\{-\ell_{\beta}(y_i, \theta)\right\}$. Generally, $\exp\left\{-\ell_{\beta}(y_i, \theta)\right\}$ is also improper for $\beta > 0$. For $y_i$ such that $f(y_i;\theta)\rightarrow 0$, it can be seen that $\exp\left\{-\ell_{\beta}(y_i, \theta)\right\} \rightarrow \exp\left\{-\frac{1}{\beta + 1}\int  f(z; \theta)^{\beta+1}dz\right\}$ which is a constant in $y_i$ and strictly greater than 0.

Figure <ref> plots $\ell_{\beta}(y_i, \theta)$ where $f(y_i; \theta) = \mathcal{N}(y_i; \mu, \sigma^2)$ and its corresponding pseudo density. Comparing this with the Tukey's loss in Figure <ref>, we can see that the Tsallis score behaves very similarly. In fact since the first version of this paper appeared, <cit.> provided methodology analogous to that developed in this paper to set $\beta$ in a data-driven manner.

Squared-error loss and Tsallis score applied to a Gaussian model (left) and corresponding (improper) densities (right). The improper densities for Tsallis score are scaled to match the mode of the Gaussian density.

§.§ Implementation Details

This sections provides a thorough description of the experiments conducted in Section <ref> and <ref> of the main paper.
Section <ref> outlines our implementation of the  in Stan.
Section <ref> discusses a restriction placed on Tukey's loss parameters to prevent degenerate parameter estimates based on including a single observation in the fit.
Given that Tukey's loss is not twice continuously differentiable, Section <ref> discusses a differentiable approximation to Tukey's loss that allows using second order sampling and optimisation methods, as well as satisfying the conditions of Theorem <ref>.
Finally, Section <ref> discusses how to set a non-local prior on Tukey's loss cut-off hyper-parameter $\kappa_2$.

Code to reproduce the examples of Sections <ref> and <ref> can be found in the repository

§.§.§ Computation

To calculate the Laplace approximations of the  we used the second order LBFGS optimisation routine implemented in the probabilistic programming language Stan <cit.>.
Further, 's implementation of the No U-Turn Sampler <cit.> can further be used to sample from the General Bayesian posteriors.
The code associated with all of our experiments is available as part of the supplementary material.

§.§.§ The breakdown of Tukey's loss

As discussed in Section <ref>, we define a restriction of the parameter space of Tukey's loss $(\beta,\sigma,\kappa_2)$ that prevents degenerate solutions where $\kappa_2$ is set to exclude all observations but one. Intuitively, the idea is to avoid $\kappa_2$ being too small. More specifically, we motive the restriction via a related notion of the breakdown point.
The breakdown of an estimator is the number of observations that can be arbitrarily perturbed without causing arbitrary changes to the estimator. The squared loss can be seen to have breakdown 0, while clearly for Tukey's-loss any observation outside the threshold can be perturbed arbitrarily without changing the estimator.
There exists abundant literature arguing that

the breakdown can be no greater than $1/2$, which intuitively means that an estimator should depend on more than half of the data.
See <cit.> and <cit.> for discussion and thorough investigation of Tukey's loss breakdown in the context of S-estimation.
The authors showed that for Tukey's breakdown to be less than $1/2$ the condition in (<ref>) must hold.

Therefore, we restrict the parameter space to $(\beta, \sigma, \kappa_2)$ satisfying (<ref>).

We note that (<ref>) applies to S-estimation rather than to our s however, we found that (<ref>) provides a simple rule that showed good performance in our examples, hence we recommend it as a default.

§.§.§ A differentiable approximation to the indicator

In order to satisfy the conditions of Theorem <ref> as well as facilitate second order optimisation (or sampling) methods, it is necessary to have a log-density that is twice continuously differentiable. While Tukey's loss itself is defined to ensure it that this is the case, the score applied to Tukey's loss, which is a function of its first and second derivatives is not. This is because the indicator $\mbox{I}(|y_i - x_i^T \beta|\leq \kappa_2\sigma)$
is not differentiable at $|y_i - x_i^T \beta| = \kappa_2\sigma$. We ensured the required conditions by the approximating

\begin{align}
|x| &\approx \sqrt{x^2 + \frac{1}{k_1}} \textrm{ for large } k_1\nonumber\\
%\mathbb{I}_{x\geq 0} &\approx \sigma_{k_2}(x) = \frac{1}{1 + \exp( - k_2x)} \textrm{ for large } k_2\nonumber
\mbox{I}(x\geq 0) &\approx \sigma_{k_2}(x) = \frac{1}{1 + \exp( - k_2x)} \textrm{ for large } k_2\nonumber
\end{align}

in (<ref>).
For our results we set $k_1 = k_2 = 100$.

§.§.§ Non-Local Prior Specification for Tukey's cut-off hyper-parameter

To place a non-local prior on $\kappa_2$, it is convenient to consider the reparametrisation  $\nu_2 = \frac{1}{\kappa_2^2}$ (as discussed in Section <ref>, this is also needed to ensure that the Hessian of Tukey's loss is finite as $\kappa_2\rightarrow\infty$ ($\nu_2\rightarrow 0$)).
Given that the Gaussian model is recovered at $\nu_2 = 0$, any prior setting positive density at $\nu_2=0$ satisfies the definition of being a local prior. In our examples we considered a half-Gaussian prior $\pi_2^{\LP}(\nu_2) \propto \mathcal{N}(\nu_2; 0, s_0^2) I(\nu_2\geq 0)$ where $I()$ is the indicator function and $s_0 = 1$.
A possible alternative that also assigns non-zero density at $\nu_2=0$ is to consider a gamma prior with shape parameter $<1$, but the associated log-prior does not have a finite Hessian as $\nu_2\rightarrow 0$ and violates Condition A3 of Theorem <ref>.

As a non-local prior we considered an inverse-gamma distribution $\pi_2^{\NLP}(\nu_2) = \mathcal{IG}(\nu_2; a_0, b_0)$, which has zero density at $\kappa_2=0$. The values of the prior hyper-parameters $(a_0, b_0)$ can have a considerable affect on the finite sample performance of the  model selection, and therefore setting these in a principled manner is important. Fortunately, $\kappa_2$ has a natural interpretation that allows restricting attention to a range of values that would be reasonable in most applications.
Recall from Section <ref> that the role of $\kappa_2$ is to exclude observations that are more than $\kappa_2$ `standard deviations' $\sigma$ away from the mean $\mu$.
First, if the data are Gaussian then there is $>0.99$ probability of an observation being within $3\sigma$'s of the mean. We thus set $(a_0,b_0)$ such that there is high prior probability that $\kappa_2 \leq 3$.
Second, given that the Gaussian has a substantial amount of central data within 1$\sigma$ from the mean, we set high prior probability that $\kappa_2 > 1$, to ensure that this central component is captured.

Given these considerations, we set default parameter values $(a_0, b_0) = (4.35, 1.56)$, which ensures a prior probability $P(\kappa_2 \in (1, 3)) = 0.95$. Figure <ref> plots these local and non-local priors for $\nu_2$ and the corresponding priors under the original parametrisation, $\kappa_2$.

Left: Density of the local half-standard Gaussian and non-local $ \mathcal{IG}(\nu_2; 4.35, 1.56)$ priors for $\nu_2$. Right: The corresponding prior density for $\kappa_2 = \nicefrac{1}{\sqrt{\nu_2}}$. The non-local prior sets $P(\kappa_2 \in (1,3))=0.95$.

Decomposition of the Inverse-Gamma Non-Local Prior Following Corollary <ref> the rate at which the  selects the simpler of nested models when supported by the data under the non-local prior specification depends on the non-local penalty function $d_l(\kappa_l) = \frac{\pi^{\NLP}_l(\kappa_l)}{\pi^{\LP}_l(\kappa_l)}$ given in  (<ref>). Given $\pi_2^{\NLP}(\nu_2) = \mathcal{IG}(\nu_2; a_0, b_0)$

and $\pi_2^{\LP}(\nu_2) = 1_{\nu_2\geq 0}\mathcal{N}(\nu_2; 0, s_0^2)$[
the corresponding local prior to an Inverse-Gamma prior would normally be a Gamma distribution, however this does not have finite second derivative as $\nu_2\rightarrow 0$ and would therefore violate the conditions of Theorem <ref>.]

the non-local penalty function is

\begin{align}
d_2(\nu_2) &= \frac{\pi_2^{\NLP}(\nu_2)}{\pi_2^{\LP}(\nu_2)}\nonumber\\
%&= \frac{\frac{b_0^{a_0}}{\Gamma(a_0)}\left(\frac{1}{\nu_2}\right)^{a_0 + 1}\exp\left(-\frac{b_0}{\nu_2}\right)}{2\frac{1}{\sqrt{2\pi}\sqrt{s_0^2}}\exp\left(- \frac{\nu_2^2}{s_0^2}\right)}\nonumber\\
&= \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)}\left(\frac{1}{\nu_2}\right)^{a_0 + 1}\exp\left(\frac{\nu_2^2}{2s_0^2} - \frac{b_0}{\nu_2}\right)\nonumber
\end{align}

and therefore

\begin{align}
- \log d_2(\nu_2) &= - \log \frac{b_0^{a_0}\sqrt{2\pi}\sqrt{s_0^2}}{2\Gamma(a_0)} + (a_0 + 1)\log \nu_2 - \frac{\nu_2^2}{2s_0^2} + \frac{b_0}{\nu_2}\nonumber
\end{align}

goes to $\inifty$ like $\frac{b_0}{\nu_2}$ as $\nu_2 \rightarrow 0$. Further for $\tilde{\nu_2} = O_p(\nicefrac{1}{\sqrt{n}})$ then we can say that

\begin{align}
\frac{b_0\sqrt{n}}{- \log d_l(\tilde{\kappa}_l)} = O_p(1)\nonumber
\end{align}

§.§.§ Asymptotic MSE

In Figure <ref> we compare the  for fixed $\kappa_2$ with an asymptotic estimate of the root mean squared error for the location parameter $\hat{\beta}(\kappa_2)$ of $\hat{\theta}_2(\kappa_2)$ minimising Tukey's loss $\ell_2(y; \theta_2, \kappa_2)$.
<cit.> propose a form of mean squared error in situations where the data generating distribution $g(y)$ has a $\epsilon$-contamination form

\begin{equation}
g(y) = (1 - \epsilon)f(y; \theta_0) + \epsilon h(y),\nonumber
\end{equation}

like that considered in Section <ref>. Here $(1-\epsilon)\times 100\%$ of the data comes from the parametric model with generating parameter $\theta_0$ and $\epsilon \times 100\%$ comes from contaminating distribution $h(\cdot)$.

In this scenario, an asymptotic approximation to the MSE of $\hat{\theta}_2(\kappa_2)$ estimated using $n$ observations from $g(y)$ is given by

\begin{align}
\mathbb{E}_g\left[\left(\hat{\theta}_2(\kappa_2) - \theta_0\right)^T\left(\hat{\theta}_2(\kappa_2) - \theta_0\right)\right] = \left(\tilde{\theta}^*_{\kappa_2} - \theta_0\right)^T\left(\tilde{\theta}^*_{\kappa_2} - \theta_0\right) + \frac{1}{n}\textrm{tr}\left\{J^{-1}_{\kappa_2}(\tilde{\theta}^*_{\kappa_2})K_{\kappa_2}(\tilde{\theta}^*_{\kappa_2})J^{-1}_{\kappa_2}(\tilde{\theta}^*_{\kappa_2})\right\},\label{Equ:AsympMSE}
\end{align}

where $\textrm{tr}\left\{\cdot\right\}$ is the trace operation and

\begin{align}
\tilde{\theta}^*_{\kappa_2} &:= \argmin_{\theta_2}\mathbb{E}_g\left[\ell_2(y; \theta_2, \kappa_2)\right]\nonumber\\
K_{\kappa_2}(\theta_2) &:= \mathbb{E}_g\left[\nabla_{\theta_2}\ell_2(y_1; \theta_2, \kappa_2)\nabla_{\theta_2}\ell_2(y_1; \theta_2, \kappa_2)^T\right]\nonumber\\
J_{\kappa_2}(\theta_2) &:= \mathbb{E}_g\left[\nabla^2_{\theta_2}\ell_2(y_1; \theta_2, \kappa_2)\right]\nonumber
\end{align}

In (<ref>), the first term is the summed squared bias of the minimiser of Tukey's loss in expectation relative to the data generating parameter $\theta_0$. The second term is the asymptotic sandwich covariance matrix of $\hat{\theta}_2(\kappa_2)$ based on $n$ observations.

In Figure <ref> we consider only the MSE of the location component $\beta$ of $\theta_2$ and therefore $\beta_0 = 0$ is the mean of the uncontaminated data. Although (<ref>) has a very different form to the  in (<ref>) of the main paper, both are optimised for similar parameter values.

§.§.§ Finite Gaussian Mixture Model Implementation

Section <ref> implemented Gaussian Mixture Models using the marginal likelihood to select the number of components. Here we provide details of the implementation used. We estimated Gaussian mixture models of the form

\begin{align}
f(y; m, \mu, \sigma, J) = \sum_{j=1}^J m_j N(y; \mu_j, \sigma_j)\nonumber
\end{align}

using Normal-Inverse-Gamma-Dirichlet priors

\begin{align}
(m_1, \ldots, m_J) &\sim \textup{Dir}(\alpha_1, \ldots, \alpha_j)\nonumber\\
\sigma_j^2 &\sim \mathcal{IG}\left(\frac{\nu_0}{2}, \frac{S_0}{2}\right), \quad j = 1,\ldots, J\nonumber\\
\mu_j | \sigma_j &\sim \mathcal{N}(0, \sqrt{\kappa}\sigma_j), \quad j = 1,\ldots, J\nonumber
\end{align}

with $\alpha_j = 1$ $j = 1,\ldots, J$, $\nu_0 = 5$, $S_0 = 0.2$ and $\kappa = 5.68$ following the recommendations of <cit.>. We first estimated the number of components as

\begin{align}
\hat{J} = \argmax_{1\leq J \leq10} \int f(y; \mu, \sigma, J)\pi(\mu |\sigma^2)\pi(\sigma^2)\pi(m){}d\mu{} d\sigma^2 dm\nonumber
\end{align}

where the integral is estimated using the posterior probability that one cluster is empty under each possible $J$ <cit.> and is implemented in the mombf <cit.> in . Then given $\hat{J}$, we found MAP estimates for $\{m_j\}_{j=1}^{\hat{J}}$, $\{\mu_j\}_{j=1}^{\hat{J}}$ and $\{\sigma_j\}_{j=1}^{\hat{J}}$.

§.§.§ The prior on the KDE bandwidth

When using the  to conduct Bayesian Kernel Density Estimation in Section <ref> we considered prior distributions of the form $\pi(h^2, w) = \mathcal{IG}(h^2; a_0, b_0)\textup{Exp}(w; \lambda_0)$. The inverse-gamma shape parameter was set to $a_0 = 2$ resulting in $\pi(h^2)$ with infinite variance. Prior parameters $b_0$ and $\lambda_0$ were elicited in order to provide regularisation towards the values that perform well when the data is Gaussian. The elicitation procedure is outlined below.

Suppose that data were truly Gaussian and consider simulating

$z \in \mathbb{R}^n$ from a standard Gaussian, i.e. $z_i \sim N(0, 1)$ for $i=1,\ldots,n$. This is consistent with the recommendations of the dirichletprocess package to  standardised the data.
Then the mean integrated squared error (MISE) of $\hat{g}_{h, w}$ with bandwidth $h$ and tempering parameter $w$ from the generating distribution $N(0, 1)$ is

\begin{equation}
\textup{MISE}[\hat{g}_{h, w}] := \mathbb{E}_{z}\left[ \int [\hat{g}_{h, w}(y^{\prime})  - N(y^{\prime}; 0, 1) ]^2 dy^{\prime}\right].\nonumber
\end{equation}

where the expectation is over $(z_1, \ldots, z_n) \overset{iid}{\sim} N(0, 1)$.
We then elicit prior parameters $b_0$ and $\lambda_0$, given $a_0 = 2$ such that

\begin{equation}
b_0, \lambda_0 := \argmin_{b, \lambda}\mathbb{E}_{\pi(h^2, w)} \textup{MISE}[\hat{g}_{h, w}].\nonumber
\end{equation}

This results in a default $\pi(h^2,w)$ centered on $(h,w)$ values such that, if the data were truly Gaussian, the estimated density would be accurate, subject to the prior being minimally informative (e.g. $h^2$ has infinite prior variance).

In the case without tempering when $w = 1$, we consider only prior $\pi(h^2) = \mathcal{IG}\left(h^2; 2, b_0\right)$ and $b_0 = 0.061$ was elicited. In the case where we also estimated tempering parameter $w$, we elicited $b_0 = 0.024$ and $\lambda_0 = 0.725$.

§.§ Supplementary results

§.§.§ Variable selection for and dataset

To illustrate our methodology for robust regression in Section <ref>, we selected between the Gaussian model and Tukey's loss when regressing the and gene expressions on a subset of the variables available in the full data sets. The procedures for which variables were selected was outlined in Sections <ref> and <ref>. To ensure that our results are reproducible, below we indicate the selected covariates and the supplementary material contains code for these variable pre-screening steps. We also provide supplementary figures providing a visual inspection of the Gaussian fit to the residuals of the and examples.

For the analysis we focused on 7 of the 10172 genes available in the data set that appear in the ‘1 pathway’ according
to the KEGGREST package in R (Tenenbaum, 2016). These were the VIT, PDE4B, ATP8B1, MAGEA11, PDE6C, PDE9A and SEPTIN4 genes.

For the analysis we selected the 15 genes with the 5 highest loadings in the first 3 principal components of the original 57 predictors. This procedure selected the following genes C15orf52, BRAT1, CYP26C1, SLC35B4, GRLF1, RXRA, RAB3GAP2, NOTCH2NL, SDC4, TTC22, PTCH2, ECH1, CSF2RA, TP53AIP1, and RRP1B.

Supplementary figures:

We provide Figure <ref> to visually inspect the Gaussian approximation to the fitted standardised residuals of the and datasets. The left panels overlay the fitted Gaussian models to the histogram of the residuals, while the right are Q-Q Normal residual plots. The top panels show that the data is well-approximated by a Gaussian model, whereas the bottom panels show that the data set exhibits heavier tails. These graphical inspections support the  findings from Section <ref> and Figure <ref>.

Top: data, where the integrated  selected the Gaussian model. Bottom: data, where the integrated  selected Tukey's loss. Left: Gaussian approximations to the residuals. Right: Q-Q normal plot of the fitted residuals according to the Gaussian model.

§.§.§ Further Gaussian Mixture Examples

We provide the results of the  bandwidth estimation for 's fixing $w = 1$, in addition to the results presented in Section <ref>.

We also provide two further Gaussian mixture data sets from <cit.>

* Asymmetric: 2-components,  $\mu_1 = 0$, $\mu_2 = 1.5$, $\sigma_1 = 1$, $\sigma_2 = 1/3$, $m_1 = 0.75$ and $m_2 = 0.25$.

* Kurtotic: 2-components,  $\mu_1 = \mu_2 = 0$, $\sigma_1 = 1$ $\sigma_2 = 0.1$, $m_1 = 2/3$ and $m_2 = 1/3$.

Figure <ref> plots the approximation to the underlying $g(y)$ of the methods under consideration in Section <ref> as the well as the  with $w = 1$, while Table <ref> estimates the Fisher's divergence of the estimates (posterior mean where appropriate) to $g(y)$. For the datasets considered in Section <ref> we see that generally learning the tempering parameter $w$ improved the density estimation relative to $w = 1$, however in most cases $w = 1$ remains competitive with the other methods considered. For both the extra examples the finite mixture model provides the best estimate according the Fisher's divergence. For the Kurtotic data set, the  's perform comparably to the and better than standard the implementations. However, for the Asymmetric data set the  's under-smooth the data and perform considerably worse than the other methods.

Fisher's divergence between the density estimates and the data-generating Gaussian mixtures in four simulation scenarios. The best performing method in each is highlighted in bold face.

(ucv)         Finite Mixture Model     $\mathcal{H}$-KDE ($w = 1$)     $\mathcal{H}$-KDE ($w \neq 1$)

Bimodal     1.03     0.37     0.13     0.10     0.26     0.09

Claw     13.77     6.09     15.25     2.17     3.37     2.51

Trimodal     0.26     0.12     0.32     0.33     0.18     0.18

Skewed     21.34     16.15     17.61     6.12     10.52     9.51

Asymmetric     0.25     3.76     0.25     0.07     48.84     2.33

Kurtotic     3.14     7.68     2.32     0.47     2.56     2.44

Density estimation for Gaussian mixture data. Histograms of observed data, standardised to 0 mean and unit variance, and estimated density by
's density function with the bandwidth rule of <cit.>, unbiased cross validated bandwidth estimation,
, Bayesian mixtures with marginal likelihood selection of components, the  estimate with no tempering ($w = 1$),
and with tempering  ($w \neq 1$). For Bayesian methods the density associated with MAP estimates is plotted.

§.§.§ The  inference for multimodal datasets

We provide two extended experiments to show that Fisher's divergence/scores' inability to learn the mixture weights for an assumed mixture model when the components are well-separated <cit.> is not an issue in our Tukey's loss or kernel density estimation examples, as discussed in Section <ref>.

Robust Regression with Tukey's loss: We extended the experiments conducted in Section <ref>, estimating Tukey's loss to outlier contaminated data. Here we considered $g(y) = 0.9\mathcal{N}(y; 0, 1) + 0.1\mathcal{N}(y; \mathbf{10}, 3)$ with the outlying mean twice as far away from the uncontaminated mean as it was in Section <ref>. The resulting $g(y)$ has an area of almost 0 density in between its two modes. The left hand side of Figure <ref> plots the $\mathcal{H}_2(y, \kappa_2)$ compared with an asymptotic approximation to the RMSE of $\hat{\beta}(\kappa_2)$ <cit.> (see Section <ref>) and the right hand side of Figure <ref> compares the fitted Gaussian and Tukey's loss densities to a histogram of the generated data

For this example $\mathcal{H}_2(y, \kappa_2)$ is maximised at  $\kappa_2 = 7.5$ and similarly to the results in Section <ref>, we see that the resulting approximation of Tukey's loss to the data captures its central component and excludes outliers. We further see that values of $\kappa_2$ that achieved high 's appear to also correspond to values of $\kappa_2$ whose Tukey's loss minimising estimates $\hat{\beta}(\kappa_2)$ have small RMSE. The separation of the modes appears not to have adverse effects on the 's performance.

Left: The integrated  $\mathcal{H}_2(y; \kappa_2)$ (black)
and asymptotic approximation to the RMSE of $\hat{\beta}(\kappa_2)$ (grey) for varying values of $\kappa_2$.
Right: Histogram of $n = 500$ observations from $g(y) = 0.9\mathcal{N}(y;0, 1) + 0.1\mathcal{N}(y;\mathbf{10}, 3)$ and fitted densities of Tukey's loss with $\hat{\kappa}_2 = 7.5$ and the Gaussian model (the height of Tukey's loss was set to match the mode of $g(y)$).

Kernel density estimation We extend the experiments in Section <ref>, estimating the bandwidth parameter for kernel density estimation. Here we consider dataset, Bimodal2, generated from

\begin{align}
%g(y)= \sum_{j=1}^J m_j N(y; \mu_j, \sigma_j)
g(y)= 0.5\mathcal{N}\left(y; -2, \frac{1}{\sqrt{6}}\right) + 0.5 \mathcal{N}\left(y; 2, \frac{1}{\sqrt{6}}\right)
\nonumber
\end{align}

Importantly, data generated form Bimodal2 has two disjoint modes with areas of almost 0 density in between. Figure <ref> plots a histogram of the generated data and the density estimated by the methods considered in Section <ref>. Table <ref> presents the Fisher's divergence from each density estimate to the true generating $g(y)$. Similarly to Section <ref> we see that the  KDE performs comparably with standard methods when the tempering parameter $w$ is estimated alongside the bandwidth.

Density estimation for Gaussian mixture data. Histograms of observed data, standardised to 0 mean and unit variance, and estimated density by 's density function with the bandwidth rule of <cit.>, unbiased cross validated bandwidth estimation, , Bayesian mixtures with marginal likelihood selection of components, the  estimate with no tempering ($w = 1$), and with tempering  ($w \neq 1$). For Bayesian methods the density associated with MAP estimates is plotted.

Fisher’s divergence between the density estimates and the data-generating Gaussian mixtures in two simulation scenarios.  The best performing method in each is highlighted in bold face.

(ucv)         Finite Mixture Model     $\mathcal{H}$-KDE ($w = 1$)     $\mathcal{H}$-KDE ($w \neq 1$)

Bimodal2     8.36     2.40     0.62     0.23     29.27     0.24