Classification has been a major task for building intelligent systems as it enables decision-making under uncertainty. Classifier design aims at building models from training data for representing feature-label distributions–either explicitly or implicitly. In many scientific or clinical settings, training data are typically limited, which makes designing accurate classifiers and evaluating their classification error extremely challenging.
While transfer learning (TL) can alleviate this issue by incorporating data from relevant source domains to improve learning in a different target domain, it has received little attention for performance assessment, notably in error estimation. In this paper, we fill this gap by investigating knowledge transferability in the context of classification error estimation within a Bayesian paradigm. We introduce a novel class of Bayesian minimum mean-square error (MMSE) estimators for optimal Bayesian transfer learning (OBTL), which enables rigorous evaluation of classification error under uncertainty in a small-sample setting. Using Monte Carlo importance sampling, we employ the proposed estimator to evaluate the classification accuracy of a broad family of classifiers that span diverse learning capabilities. Experimental results based on both synthetic data as well as real-world RNA sequencing (RNA-seq) data show that our proposed OBTL error estimation scheme clearly outperforms standard error estimators, especially in a small-sample setting, by tapping into the data from other relevant domains.

§ INTRODUCTION

Transfer learning provides  promising means to repurpose the data and/or scientific knowledge available in other relevant domains for new applications in a given domain. The ability to transfer relevant data/knowledge across different domains practically enables learning effective models in target domains with limited data. Classifier design can take advantage of transfer learning (TL) to address small-sample challenges we often face in various scientific applications. However, rigorous error estimators that can leverage such transferred data/knowledge for better estimation of classification error have been missing to date, which makes the design framework epistemologically incomplete [1].
Generally, the scientific validity of any predictive model is assessed by the ability to generalize outside the observed training sample. However, the available sample is often too small in many scientific applications (e.g., bio-marker discovery) to hold out sufficient data just for testing purpose, which makes the reuse of training data for both classifier design and error estimation inevitable.

While various error estimation schemes exist to date, their accuracy and reliability in a small-sample setting are often questioned [2]. For instance, in [3] many classification studies of cancer gene expression data have been listed where the performance was assessed by cross-validation (CV) based on small-size training datasets. Analyses in [4] have shown that CV error estimators derived based on small-size samples show large variance, which explains the controversy across many biological studies that relied on data-driven CV [5]. Model-based error estimation also faces practical challenges as non-informative modeling assumptions may mislead the error estimators in case of model mismatch.

The ability for accurate error estimation based on small samples is also critical in other contexts, an example being continual learning [10], where a series of labeled datasets are sequentially fed to the learner as in realistic learning scenarios.
In recent years, continual learning regained attention as a promising strategy for avoiding “catastrophic forgetting” that may arise when the training data are split for a series of small learning operations called tasks [7]. Such a continual learning setting is getting prevalent these days, where retaining the observed training data is either undesirable (confidentiality) or intractable (high-throughput systems), and developing reliable task-specific error estimators is indispensable. For instance, an intuitive approach to continual learning from a Bayesian perspective is to leverage the posterior of the current task to update the prior of the next task [6]. However, analysis in [8] has shown that evaluation approaches for this prior-focused set-up suffer from severe bias in realistic scenarios, particularly for finely partitioned data. Recent work in [9] provided a solution for test data scarcity by reusing the same test set in the context of a continuously evolving classification problem. To avoid overfitting the test data, the authors employed a reusable holdout mechanism based on the AUC (area under the receiver operating characteristic curve) metric. Nevertheless, this approach remains contingent on the availability of an independent test set.
For these reasons, there is a pressing need to develop novel error estimators that can effectively overcome data-scarcity limitations. For assessing different classification models in the context of small-size training datasets, having an accurate error estimator with TL capabilities that can take advantage of relevant datasets in other domains would be highly beneficial. Such an estimator would be readily applicable to continual learning as cross-task datasets can be seen as related source-target samples. In the remainder of this introductory section, we provide a brief review of the standard error estimation techniques along with prevalent transfer learning scenarios.

§.§ Relevant prior studies

§.§.§ Error Estimation

For unknown feature-label distributions, the classification error of a given classifier is typically estimated by leveraging a large sample collected from the true distribution. However, limiting factors, such as the excessive cost of large-scale data acquisition, make it often infeasible to collect and hold out large test sets. Consequently, the available small-size sample may have to be used for both training and assessing the classifier, and researchers have striven to devise practical methods for accurate error estimation.

Existing error estimation schemes can be broadly categorized into parametric and non-parametric methods. Non-parametric estimators compute the error rate by counting the misclassified points, where widely used estimators include the resubstitution, cross-validation (CV), and bootstrap estimators. Resubstitution assesses the error committed on the training data, as a result of which it has an optimistic bias and frequently raises overfitting issues [11].
In $k$-fold CV [12, 13], the data is randomly split into $k$-folds and the classification error is estimated by the average error of the $k$ distinct classifiers, each trained on $k$-1 folds and evaluated on the left-out fold. Leave-one-out CV (LOO) is a non-randomized special case of CV where the number of folds $k$ is equal to the sample size. Unlike resubstitution, CV and LOO tend to have small bias but they suffer from large variance [4].
The bootstrap estimator [14, 15] takes a different randomized approach. Instead of partitioning the data into $k$-folds, bootstrap generates $k$ sets of equal size by sampling (with replacement) from the original set.

In this sampling approach, bootstrap sets are formed on average by 63.2% of unique data points from the original training data. This randomization, while reducing the variance, results in a pessimistic bias as it excludes parts of the data from model training. To reconcile this deviation, the 0.632-bootstrap uses a weighted average of the standard bootstrap and the optimistically biased resubstitution with 63.2% of the approximation being a bootstrap estimate.
Bolstered error estimation [16] employs density kernels centered at validation data points to smoothly approximate the true error. Instead of counting the misclassified points, bolstered estimator quantifies the overflow of bolstering kernels through the classifier decision boundaries. This technique, when associated with optimal kernel variances, has shown to result in a significant improvement over resubstitution and CV in terms of variance and bias.

Parametric methods include the popular plug-in estimator that naively estimates the true error from an empirical model. One major limitation for this technique is the strong dependence on the estimated parameters, which may fail drastically due to poor parameter estimation.
The Bayesian minimum mean square error estimator (BEE) proposed in [3, 17] significantly enhances the robustness by computing the expected true error with respect to the posterior of the model parameters. The BEE has shown notable improvements over standard estimators as it effectively handles the uncertainty about the feature-label distributions [3, 17].

§.§.§ Transfer Learning

Transfer Learning (TL) was originally proposed to provide remedies for pitfalls caused by training-data scarcity in a target domain by utilizing available data from different yet relevant source domains [18]. Based on the properties of source and target domains, two scenarios of TL are available. The first one, commonly known as “homogeneous TL”, occurs when the source and target domains share the same feature space. The second scenario is called “heterogeneous TL” and is considered when differences exist between domains in terms of feature space or data dimensionality.

Practically, TL is applied across domains that share identical category structures, where the sample and label spaces remain unchanged while only the probability distributions differ across domains. This specific case of TL is also known as domain adaptation (DA) [19, 20]. In order to maximize the similarity between feature-label distributions in the source and target domains, the vast majority of existent DA methods proceed by adapting three domain aspects: the data, the model parameters, or a cross-domain latent space.
For data adaptation, methods such as instance re-weighting [21] and transfer adaptive boosting [22] are mostly suited for homogeneous TL as they act by assigning optimized weights to source and target samples to build an augmented dataset. In [23], a Bayesian framework was proposed for TL, in which the source and target domains are related through the joint prior density of the model parameters.
For parameter adaptation, popular algorithms redesign classifier models built in source domains to make them effective in target domains. This is enabled by leveraging the available target data to update the classifier parameters [24, 25, 26].
Another approach for TL is to project the source and target data into an intermediate invariant space that unifies the two domains [27]. This approach showed promising results in heterogeneous TL [28], taking advantage of deep learning models that can map complex feature spaces in the respective domains into an invariant latent domain [29, 30].

§.§ Main contributions of this paper

In this paper, we propose a novel transfer learning framework for robust estimation of classification error based on a rigorous Bayesian paradigm. To the best of our knowledge, this study is the first work on TL-based BEE, which can significantly enhance our understanding of transferability across domains in the context of error estimation. Building on the Bayesian transfer learning framework proposed in [23], we introduce a novel TL-based BEE estimator that can enhance the error estimation accuracy in the target domain by utilizing the data available in a relevant source domain. We present a rigorous study of error estimation in the context of Bayesian TL and show that our proposed TL-based BEE effectively represents and exploits the relatedness (or dependency) between different domains to improve error estimates in a challenging small-sample setting.
For applicability of the proposed TL-based BEE estimator in real-world problems for arbitrary classifiers, we introduce an efficient and robust importance sampling set-up with control variates where the importance density and the control variates function are carefully defined to reduce the variance of the estimator while keeping the overall sampling process computationally feasible and scalable. For this purpose, we utilize Laplace approximations for fast evaluation of matrix-variate confluent and Gauss hypergeometric functions.
The performance of the TL-based BEE estimator is extensively evaluated using both synthetic datasets as well as real-world biological datasets. A wide range of classifiers with different levels of learning capabilities are considered, in order to demonstrate the general applicability of our TL-based BEE estimation scheme. We also show the outstanding performance of the proposed estimator with respect to standard error estimation techniques that are commonly used.

The remainder of this paper is organized as follows. In Sec. <ref> we present the Bayesian transfer learning validation framework. Section <ref> defines the novel TL-based Bayesian MMSE estimator and derives its expression for arbitrary classifiers. In Sec. <ref>, details of the robust importance sampling set-up is provided with the definition of the control variates. Section <ref> elaborates the experimental set-up for both synthetic and real-world datasets. The performance analysis results are presented in Sec. <ref> and we conclude the paper in Sec. <ref> with potential future research directions.
Appendix A reviews the binary classification setting. Appendix B states the definition of the standard Bayesian MMSE estimator. Appendix C recalls the definition of matrix-variate hypergeometric functions and provides the Laplace approximations of confluent and Gauss hypergeometric functions of matrix argument. Appendix D states some useful results for importance sampling and provides the derivation of the control variates function. Finally, Appendix E provides additional results for linear classifiers.

§ BAYESIAN TRANSFER LEARNING FRAMEWORK FOR BINARY CLASSIFICATION

We consider a binary classification problem in the context of supervised transfer learning where there are two common classes in each domain.
Let $\mathcal{D}_{s}$ and $\mathcal{D}_{t}$ be two labeled datasets from the source and target domains with sizes $N_{s}$ and $N_{t}$, respectively.
We are interested in the scenario where $N_t \ll N_{s}$.
Let $\mathcal{D}_{s}^{y}=\left\{ \mathbf{x}_{s,1}^{y}, \mathbf{x}_{s,2}^{y}, \cdots, \mathbf{x}_{s,n_{s}}^{y}\right \}$, $y \in \left \{0,1 \right \}$, where $n_{s}^{y}$ denotes the size of source data in class $y$.
Likewise, let $\mathcal{D}_{t}^{y}=\left\{ \mathbf{x}_{t,1}^{y}, \mathbf{x}_{t,2}^{y}, \cdots, \mathbf{x}_{t,n_{t}}^{y}\right \}$, $y \in \left \{0,1 \right \}$, where $n_{t}^{y}$ denotes the size of target data in class $y$.

We consider a $d$-dimensional homogeneous transfer learning scenario where $\mathcal{D}_{s}$ and $\mathcal{D}_{t}$ are normally distributed and separately sampled from the source and target domains, respectively.

\begin{equation}
\mathbf{x}_{z}^{y}  \sim \mathcal{N}\left ( \mathbf{\mu}_{z}^{y},\left( \mathbf{\Lambda}_{z}^{y}\right)^{-1}\right ),~y \in \left\{0,1\right\},
\end{equation}

where $z\in\left\{s,t\right\}$, $\mathbf{\mu}_{z}^{y}$ is a $\left(d\times 1\right)$ mean vector in domain $z$ for class $y$, and $\mathbf{\Lambda}_{z}^{y}$ is a $\left(d \times d\right)$ matrix that denotes the precision matrix (inverse of covariance) in domain $z$ for label $y$.
An augmented feature vector $\mathbf{x}^{y}=\begin{bmatrix}\mathbf{x}_{t}^{y}\\\mathbf{x}_{s}^{y}\end{bmatrix}$ is a joint sample point from two related source and target domains given by

\begin{equation}
\mathbf{x}^{y}\sim \mathcal{N}\left ( \mathbf{\mu}^{y},\left( \mathbf{\Lambda}^{y}\right)^{-1}\right ),~y \in \left\{0,1\right\},
\end{equation}

\begin{equation}
\mathbf{\mu}^{y}=\begin{bmatrix}\mathbf{\mu}_{t}^{y}\\\mathbf{\mu}_{s}^{y}\end{bmatrix}, \mathbf{\Lambda}^{y}=\begin{bmatrix}
\mathbf{\Lambda}_{t}^{y} & \mathbf{\Lambda}_{ts}^{y}\\
{\mathbf{\Lambda}_{ts}^{y}}^{T} & \mathbf{\Lambda}_{s}^{y}
\end{bmatrix}
\end{equation}

where $X^{T}$ denotes the transpose of matrix $X$.
This sampling is enabled through a joint prior distribution for $\mathbf{\Lambda}_{s}^{y}$ and $\mathbf{\Lambda}_{t}^{y}$ that marginalizes out the off-diagonal block matrix $\mathbf{\Lambda}_{ts}^{y}$.
Using a Gaussian-Wishart distribution as the joint prior for mean and precision matrices, the joint model factorizes as

\begin{equation}\label{joint_prior}
p\left(\mathbf{\mu}_{s}^{y}, \mathbf{\mu}_{t}^{y}, \mathbf{\Lambda}_{s}^{y}, \mathbf{\Lambda}_{t}^{y}\right)=p\left(\mathbf{\mu}_{s}^{y}, \mathbf{\mu}_{t}^{y}|\mathbf{\Lambda}_{s}^{y}, \mathbf{\Lambda}_{t}^{y}\right)~p\left(\mathbf{\Lambda}_{s}^{y}, \mathbf{\Lambda}_{t}^{y}\right).
\end{equation}

For conditionally independent mean vectors given the covariances, the joint prior in (<ref>) further expands to

\begin{equation}
p\left(\mathbf{\mu}_{s}^{y}, \mathbf{\mu}_{t}^{y}, \mathbf{\Lambda}_{s}^{y}, \mathbf{\Lambda}_{t}^{y}\right)=p\left(\mathbf{\mu}_{s}^{y}|\mathbf{\Lambda}_{s}^{y}\right)~p\left(\mathbf{\mu}_{t}^{y}|\mathbf{\Lambda}_{t}^{y}\right)~p\left(\mathbf{\Lambda}_{s}^{y}, \mathbf{\Lambda}_{t}^{y}\right).
\end{equation}

The block diagonal precision matrices $\mathbf{\Lambda}_{z}^{y}$ for $z \in \left \{t,s \right \}$ are obtained after sampling $\mathbf{\Lambda}^{y}$ from a predefined joint Wishart distribution as defined in [23] such that
$\mathbf{\Lambda}^{y} \sim W_{2d}\left(\mathbf{M}^{y},\nu^{y} \right)$, where $\nu^{y}$ is a hyperparameter for the degrees of freedom that satisfies $\nu^{y} \geq 2d$ and $\mathbf{M}^{y}$ is a $\left(2d \times 2d\right)$ positive definite scale matrix of the form $\mathbf{M}^{y}=\begin{pmatrix}
\mathbf{M}_{t}^{y} & \mathbf{M}_{ts}^{y}\\
{\mathbf{M}_{ts}^{y}}^{T} & \mathbf{M}_{s}^{y}
\end{pmatrix}$.
$\mathbf{M}_{t}^{y}$ and $\mathbf{M}_{s}^{y}$ are also positive definite scale matrices and $\mathbf{M}_{ts}$ denotes the off-diagonal component that models the interaction between source and target domains.
Given $\mathbf{\Lambda}_{z}^{y}$, and assuming normally distributed mean vectors we get

\begin{equation}
\mathbf{\mu}_{z}^{y}\sim \mathcal{N}\left ( \mathbf{m}_{z}^{y},\left( \kappa_{z}^{y}~\mathbf{\Lambda}_{z}^{y}\right)^{-1}\right ),~z \in \left\{s,t\right\} \mbox{~and~} y \in \left\{0,1\right\},
\end{equation}

where $\mathbf{m}_{z}^{y}$ is the $\left(d \times 1\right)$ mean vector of the mean parameter $\mathbf{\mu}_{z}^{y}$ and $\kappa_{z}^{y}$ is a positive scalar hyperparameter.

The joint prior distribution $p\left(\mathbf{\Lambda}_{s}^{y},\mathbf{\Lambda}_{t}^{y}\right)$ as derived in [23] acts like a bridge through which the useful knowledge transfers from the source to the target domain, making the posterior of the target parameters of the underlying feature-label distribution distributed more narrowly  around the true values.

§ BAYESIAN MMSE ESTIMATION VIA TRANSFER LEARNING

In this section we propose a novel class of Bayesian MMSE error estimators for transfer learning where the observed sample is a mixture of source and target data. The basic  classification setting and a brief review of the standard BEE estimator are presented in the Appendices A and B.

Rooted in signal estimation, the BEE has been motivated by optimal filtering for functions of random variables [3]. For a function of two random variables $g\left(X,Y\right)$, the optimal estimator $\hat{g}\left(Y\right)$ of a filter $g\left(Y\right)$ after observing only $Y$ in the mean-square sense is given by

\begin{equation}
\hat{g}\left(Y\right)=E_{X}\left[g\left(X,Y\right)|Y\right].
\end{equation}

Replacing $X$ with the parameter vector $\theta$ of the feature-label distribution and $Y$ by the sample $S_{n}$ (of size $n$), leads to the standard BEE that has been introduced in [3] as

\begin{equation}
\hat{\varepsilon}\left (S_{n} \right )=E_{\mathbf{\theta}}\left[\varepsilon_{n}\left( \mathbf{\theta},S_{n}~|~S_{n}\right) \right].
\end{equation}

In TL, the sample $S_{n}$ is a mixture of source and target data such that $S_{n}=\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n}$ with $n=N_{s}+N_{t}$ and the classifier $\psi_{n}$ is designed either on $\mathcal{D}_{t}$, $\mathcal{D}_{s}$, or $\mathcal{D}_{s}\cup \mathcal{D}_{t}$. This requires a close attention as the TL-based BEE is valid only for fixed classifiers given the sample. This assumption carries limitations. For instance, classifiers that are only fixed given $\mathcal{D}_{t}$ but not $\mathcal{D}_{s}$ are not deterministic for every set of parameters estimated based on $\mathcal{D}_{s}\cup \mathcal{D}_{t}$. We introduce in this paper the novel TL-based BEE for TL defined as

\begin{equation}\label{BEE_TL_long}
\hat{\varepsilon}\left (\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n} \right )=E_{\mathbf{\theta}}\left[\varepsilon_{n}\left( \mathbf{\theta},\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n}\right)~|~\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n} \right],
\end{equation}

where $\theta=\left[\theta_{t}, \theta_{s}\right]$ denotes the parameter vector of the joint model formed by the target parameters $\theta_{t}$ and source parameters $\theta_{s}$. For a fixed classifier given $\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n}$, this estimator is optimal on average in the mean-square sense and unbiased when averaged over all parameters and samples.
For classification in the target domain, the posterior density $\pi^{*}\left(\theta\right)$ reduces to the posterior of the target parameters after observing the target and source data and takes the form

\begin{equation}
\pi^{*}\left(\theta_{t}\right)=\pi^{*}\left(\theta_{t}~|~\mathcal{D}_{s}, \mathcal{D}_{t}\right),
\end{equation}

where $\pi^{*}\left(\theta_{t}~|~\mathcal{D}_{s}, \mathcal{D}_{t}\right)$ is obtained by marginalizing out the source domain parameters.
Ultimately, the BEE for TL takes the form

\begin{align}
\begin{split}
\hat{\varepsilon}\left (\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n} \right )&=E_{\mathbf{\theta}_{t}}\left[\varepsilon_{n}\left( \mathbf{\theta}_{t},\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n}\right)~|~\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n} \right]\\
&=E_{\pi^{*}\left(\mathbf{\theta}_{t}\right)}\left[\varepsilon_{n}\left( \mathbf{\theta}_{t},\left(\mathcal{D}_{s}\cup \mathcal{D}_{t}\right)_{n}\right)\right].
\end{split}
\end{align}

For the sake of simplicity we write

\begin{equation}
\hat{\varepsilon}=E_{\pi^{*}}\left[\varepsilon_{n}\right],
\end{equation}

where $\pi^{*}=\pi^{*}\left(\mathbf{\theta}_{t}~|~\mathcal{D}_{t},\mathcal{D}_{s}\right)$ denotes the posterior of the target parameters after observing the hybrid sample $\mathcal{D}_{t} \cup \mathcal{D}_{s}$.

§.§ Bayesian transfer learning for error estimation

The advantage of the mathematical formulation that underlies the proposed TL-based BEE (and also the original TL Bayesian framework in [23]) is that it articulates a unified Bayesian inference model that assumes a specified prior distribution governing the parameter vector $\theta_{t}$ and acting like a bridge to help update $\pi^{*}\left(\theta_{t}\right)$ after observing $\mathcal{D}_{t}$ and $\mathcal{D}_{s}$. From this standpoint, the derivation of the TL-based BEE for TL depends on determining $\pi^{*}\left(\theta_{t}\right)$.
To determine the TL-based BEE in the context of the presented Bayesian transfer learning framework we evoke the following theorem.

Theorem 1 [23]: Given the target $\mathcal{D}_{t}$ and source $\mathcal{D}_{s}$ data, the posterior distribution of target mean $\mu_{t}^{y}$ and the target precision matrix $\mathbf{\Lambda}_{t}^{y}$ for the classes $y\in\left\{0,1\right\}$ has Gaussian-hypergeometric function distribution given by

\begin{align}
\begin{split}
\end{split}
\end{align}

where $A^{y}$ is a constant of proportionality given by

\begin{align}
\begin{split}
\left(A^{y}\right)^{-1}&=\left(\frac{2\pi}{\kappa_{t,n}^{y}}\right)^{\frac{d}{2}}2^{\frac{d\left(\nu^{y}+n_{t}^{y}\right)}{2}}\Gamma_{d}\left(\frac{\left(\nu^{y}+n_{t}^{y}\right)}{2}\right)\left|\mathbf{T}_{t}^{y}\right|^{\frac{\left(\nu^{y}+n_{t}^{y}\right)}{2}}\\
&\times~\pFq{2}{1}{\frac{\nu^{y}+n_{s}^{y}}{2}, \frac{\nu^{y}+n_{t}^{y}}{2}}{\frac{\nu^{y}}{2}}{\mathbf{T}_{s}^{y}~\mathbf{F}^{y}~\mathbf{T}_{t}^{y}~{\mathbf{F}^{y}}^{T}},
\end{split}
\end{align}

\begin{align}
\begin{split}
\kappa_{t,n}^{y}&=\kappa_{t}^{y}+n_{t}^{y},\\
\mathbf{m}_{t,n}^{y}&=\frac{\kappa_{t}^{y}~\mathbf{m}_{t}^{y}+n_{t}^{y}~\overline{\mathbf{x}}_{t}^{y}}{\kappa_{t}^{y}+n_{t}^{y}},\\
\left(\mathbf{T}_{t}^{y}\right)^{-1}&=\left(\mathbf{M}_{t}^{y}\right)^{-1}+{\mathbf{F}^{y}}^{T}~\mathbf{C}^{y}~\mathbf{F}^{y}+\mathbf{S}_{t}^{y}\\
\left(\mathbf{T}_{s}^{y}\right)^{-1}&=\left(\mathbf{C}^{y}\right)^{-1}+\mathbf{S}_{s}^{y}\\
\end{split}
\end{align}

with sample means and covariances for $z\in\left\{s,t\right\}$ given by

\begin{align*}
\begin{split}
\overline{\mathbf{x}}_{z}^{y}&=\frac{1}{n_{z}^{y}}\sum_{i=1}^{n_{z}^{y}}\mathbf{x}_{z,i}^{y}\\
\mathbf{S}_{z}^{y}&=\sum_{i=1}^{n_{z}^{y}}\left(\mathbf{x}_{z,i}^{y}-\overline{\mathbf{x}}_{z}^{y}\right)\left(\mathbf{x}_{z,i}^{y}-\overline{\mathbf{x}}_{z}^{y}\right)^{T}.
\end{split}
\end{align*}

$\pFq{1}{1}{a}{b}{\mathbf{X}}$ and $\pFq{2}{1}{a,b}{c}{\mathbf{X}}$ are respectively the confluent and Gauss matrix-variate hypergeometric functions reviewed in Appendix C.

Now, using Theorem 1 and assuming that $c$, $\theta_{t}^{0}$, and $\theta_{t}^{1}$ are independent prior to observing $\mathcal{D}_{t}$ and $\mathcal{D}_{s}$, the BEE for TL is given by

\begin{equation}
\hat{\varepsilon}=E_{\pi^{*}}\left[c \right]~E_{\pi^{*}}\left[\varepsilon_{n}^{0} \right]+\left(1-E_{\pi^{*}}\left[c \right]\right)~E_{\pi^{*}}\left[\varepsilon_{n}^{1} \right],
\end{equation}

\begin{equation}\label{BEE_final}
\end{equation}

with $\Theta_{t}^{y}$ being the parameter space that contains all possible values for $\theta_{t}^{y}$.

§.§ Computing TL-based BEE for arbitrary classifiers

Computing the TL-based BEE for an arbitrary classifier $\psi_{n}$ involves the evaluation of the integral in (<ref>).
Even when we have an analytic expression for the true error of the studied classifier, the closed-form expression for the TL-based BEE cannot be easily derived due to the complex expression of the target posterior in the presence of the matrix-variate hypergeometric functions. With non-linear classifiers, this becomes practically impossible as no closed-form expression exists for the true error itself. The standard way to approximate the true error in this case is to consider the test error. For a specified parameter $\theta_{t}$, a large test set is generated from $f_{\theta_{t}}\left(\mathbf{x},y\right)$ and the performance of $\psi_{n}$ is evaluated on that test set. This requires sampling from $\pi^{*}\left(\theta_{t}^{y}\right)$ so that the integral in (<ref>) can be approximated by a finite sum. Suppose we have $N$ posterior sample points $\theta_{t,i}^{y}\sim \pi^{*}\left(\theta_{t}^{y}\right),~i=1\cdots N$. Then the approximation is given by

\begin{equation}
\end{equation}

Because of the generalized confluent and Gauss hypergeometric functions in the expression of $\pi^{*}$, sampling directly from the posterior is very laborious and the computational cost of applying Markov Chain Monte Carlo (MCMC) methods is exorbitant as the execution may take several weeks even on high performance computing clusters. To address this issue, we propose in the next section an efficient self-normalized importance sampling set-up with control variates that provides accurate estimates for the TL-based BEE and significantly reduces the computation time to make the proposed TL-based BEE feasible.

§ SELF-NORMALIZED IMPORTANCE SAMPLING WITH CONTROL VARIATES

§.§ Importance sampling

Importance sampling (IS) is a variance reduction technique that provides a remedy to sampling from complex distributions [33]. To estimate $E_{\pi^{*}}\left[\varepsilon_{n}^{y}\right]$, IS makes a multiplicative adjustment to $\varepsilon_{n}^{y}$ to compensate for sampling from an alternative importance distribution $\Phi^{*}$ instead of $\pi^{*}$. If $\Phi^{*}$ is a positive probability density function on $\Theta_{t}^{y}$, we can write

\begin{align}\label{importance_expectation}
\begin{split}
\end{split}
\end{align}

Achieving an accurate IS estimation is contingent on selecting an appropriate importance density that is nearly proportional to $\varepsilon_{n}^{y}\left(\theta_{t}^{y}\right)~\pi^{*}\left(\theta_{t}^{y}\right)$. By analogy to [34, 35], a plausible and cogent candidate for $\Phi^{*}$ emanates as the posterior of target parameters upon observation of target-only data. Obviously, both distributions are tracking the same model parameters in the target domain upon observation of data. To determine $\Phi^{*}\left(\theta_{t}^{y}\right)= p\left(\mu_{t}^{y},\mathbf{\Lambda}_{t}^{y}~|~\mathcal{D}_{t}^{y}\right)$ we require the following lemma:

Lemma 1 [36]: If $\mathcal{D}=\left\{\mathbf{x}_{1},\cdots,\mathbf{x}_{n}\right\}$ where $\mathbf{x}_{i}$ is a $d\times 1$ vector and $\mathbf{x}_{i}\sim\mathcal{N}\left(\mu, \left(\mathbf{\Lambda}\right)^{-1}\right)$, for $i=1,\cdots,n$, and $\left(\mu, \mathbf{\Lambda}\right)$ has a Gaussian-Wishart prior, such that $\mu|\mathbf{\Lambda}\sim \mathcal{N}\left(\mathbf{m},\left(\kappa\mathbf{\Lambda}\right)^{-1}\right)$ and $\mathbf{\Lambda\sim W_{d}\left(\mathbf{M},\nu\right)}$, then the posterior of $\left(\mu,\mathbf{\Lambda}\right)$ upon observing $\mathcal{D}$ is also a Gaussian-Wishart distribution such that

\begin{align}
\begin{split}
\mu|\mathbf{\Lambda},\mathcal{D}&\sim\mathcal{N}\left(\mathbf{m}_{n},\left(\kappa_{n}\mathbf{\Lambda}\right)^{-1}\right);\,\,%\\
\mathbf{\Lambda}|\mathcal{D}\sim W_{d}\left(\mathbf{M}_{n},\nu_{n}\right),
\end{split}
\end{align}

\begin{align}
\begin{split}
\kappa_{n}&=\kappa+n, \quad%\\
\nu_{n}=\nu+n,\\
\mathbf{m}_{n}&=\frac{\kappa \mathbf{m}+n\overline{\mathbf{x}}}{\kappa+n}, \mbox{ and} \\
\mathbf{M}_{n}^{-1}&= \mathbf{M}^{-1}+\mathbf{S}+\frac{\kappa n}{\kappa+n}\left(\mathbf{m}-\overline{\mathbf{x}}\right)\left(\mathbf{m}-\overline{\mathbf{x}}\right)^{T},
\end{split}
\end{align}

depending on the sample mean and covariance matrix

\begin{align}
\begin{split}
\overline{\mathbf{x}}&=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}, \quad%\\
\mathbf{S}=\sum_{i=1}^{n}\left(\mathbf{x}_{i}-\overline{\mathbf{x}}\right)\left(\mathbf{x}_{i}-\overline{\mathbf{x}}\right)^{T}.
\end{split}
\end{align}

Using Lemma 1 we now get the expression of the importance density $\Phi^{*}$ given by

\begin{align}
\begin{split}
\end{split}
\end{align}

\begin{align}\label{empirical_estimates}
\begin{split}
\kappa_{t,n}^{y}&=\kappa_{t}^{y}+n_{t}^{y},\\
\mathbf{m}_{t,n}^{y}&=\frac{\kappa_{t}^{y}~\mathbf{m}_{t}^{y}+n_{t}^{y}~\overline{\mathbf{x}}_{t}^{y}}{\kappa_{t}^{y}+n_{t}^{y}},\\
\left(\mathbf{M}_{t,n}^{y}\right)^{-1}&=\left(\mathbf{M}_{t}^{y}\right)^{-1}+\mathbf{S}_{t}^{y}\\
\end{split}
\end{align}

with sample mean and covariance given by

\begin{align*}
\begin{split}
\overline{\mathbf{x}}_{t}^{y}&=\frac{1}{n_{t}^{y}}\sum_{i=1}^{n_{t}^{y}}\mathbf{x}_{t,i}^{y}\\
\mathbf{S}_{t}^{y}&=\sum_{i=1}^{n_{t}^{y}}\left(\mathbf{x}_{t,i}^{y}-\overline{\mathbf{x}}_{t}^{y}\right)\left(\mathbf{x}_{t,i}^{y}-\overline{\mathbf{x}}_{t}^{y}\right)^{T}.
\end{split}
\end{align*}

After simplifications, the expression of the TL-based BEE in (<ref>) takes the form

\begin{equation}
\end{equation}

where $\theta_{t}^{y}=\left(\mu_{t}^{y}, \mathbf{\Lambda_{t}^{y}}\right)$ and $\mathcal{L}\left(\theta_{t}^{y}\right)$ is the likelihood ratio given by

\begin{align}
\begin{split}
&\mathcal{L}\left(\mu_{t}^{y}, \mathbf{\Lambda_{t}^{y}}\right)=\mbox{etr}\left(-\frac{1}{2}\left[\left(\mathbf{T}_{t}^{y}\right)^{-1}-\left(\mathbf{M}_{t,n}^{y}\right)^{-1} \right ]\mathbf{\Lambda}_{t}^{y}\right)\\
&\times~\left|\frac{\mathbf{M}_{t,n}^{y}}{\mathbf{T}_{t}^{y}}\right|^{\frac{\left(\nu^{y}+n_{t}^{y}\right)}{2}}\frac{\pFq{1}{1}{\frac{\nu^{y}+n_{s}^{y}}{2}}{\frac{\nu^{y}}{2}}{\frac{1}{2}~\mathbf{F}^{y}~\mathbf{\Lambda}_{t}^{y}~{\mathbf{F}^{y}}^{T}~\mathbf{T}_{s}^{y}}}{\pFq{2}{1}{\frac{\nu^{y}+n_{s}^{y}}{2}, \frac{\nu^{y}+n_{t}^{y}}{2}}{\frac{\nu^{y}}{2}}{\mathbf{T}_{s}^{y}~\mathbf{F}^{y}~\mathbf{T}_{t}^{y}~{\mathbf{F}^{y}}^{T}}}.
\end{split}
\end{align}

Although the likelihood ratio has a simplified expression, computing the hypergeometric functions involves the computation of series of zonal polynomials, which is computationally expensive and not scalable to high dimensions. To mitigate this limitation, we use the Laplace
approximations of these functions (see Appendix C). To rectify possible disproportionalities in likelihood ratios due to approximations, we consider the self-normalized importance sampling estimate given by

\begin{equation}
\widehat{E}_{\Phi^{*}}\left[\varepsilon_{n}^{y}\right]\approx\frac{\sum_{i=1}^{N}\varepsilon_{n}^{y}\left(\theta_{t,i}^{y}\right)~\mathcal{L}\left(\theta_{t,i}^{y}\right)}{\sum_{i=1}^{N}~\mathcal{L}\left(\theta_{t,i}^{y}\right)}
\end{equation}

with $\theta_{t,i}^{y}\sim \Phi^{*}\left(\theta_{t}^{y}\right),~i=1\cdots N$.

§.§ Control variates

For more stable and efficient estimates, we further combine importance sampling with control variates.
Using control variates in conjunction with importance sampling is a variance reduction technique, in particular when a significant portion of a model for estimating the expectation can be solved explicitly.
In our case, a useful control variates function (CVF) $\mathcal{V}\left(\theta_{t}^{y}\right)$ satisfies

\begin{align}
\begin{split}
\end{split}
\end{align}

where $\delta$ is a constant. Under such circumstances, a more stable estimate for the TL-based BEE can be derived as

\begin{equation}
\widetilde{E}_{\Phi^{*}}\left[\varepsilon_{n}^{y}\right]\approx\frac{\sum_{i=1}^{N}\varepsilon_{n}^{y}\left(\theta_{t,i}^{y}\right)~\mathcal{L}\left(\theta_{t,i}^{y}\right)}{\sum_{i=1}^{N}~\mathcal{L}\left(\theta_{t,i}^{y}\right)}-\frac{1}{N}\sum_{i=1}^{N}\frac{\beta\mathcal{V}\left(\theta_{t,i}^{y}\right)}{\Phi^{*}\left(\theta_{t,i}^{y}\right)}+\beta\delta,
\end{equation}

where $\theta_{t,i}^{y}\sim \Phi^{*}\left(\theta_{t}^{y}\right),~i=1\cdots N$ and $\beta$ is a weighting coefficient tuned to reduce the variance of the estimate. The optimal value of $\beta$ is given by

\begin{equation}
\beta_{opt}=\frac{\mbox{cov}\left[\zeta_{n}^{y}\left(\theta_{t}^{y}\right),\mathcal{V}\left(\theta_{t}^{y}\right)\right]}{\mbox{var}\left[\mathcal{V}\left(\theta_{t}^{y}\right)\right]},
\end{equation}

\begin{equation}
\zeta_{n}^{y}\left(\theta_{t}^{y}\right)=\frac{\varepsilon_{n}^{y}\left(\theta_{t}^{y}\right)~\mathcal{L}\left(\theta_{t}^{y}\right)}{\frac{1}{N}\sum_{i=1}^{N}~\mathcal{L}\left(\theta_{t,i}^{y}\right)}
\end{equation}

and $\mbox{cov}\left[\cdot, \cdot\right]$ and $\mbox{var}\left[\cdot\right]$ denote covariance and variance, respectively (see Appendix D for more details). In practice, it is not likely that we know $\beta_{opt}$ beforehand, but it is estimated from the Monte Carlo sample. It turns out that $\widetilde{E}_{\Phi^{*}}$ has lower variance than $\widehat{E}_{\Phi^{*}}$ by a factor of $\left(1-\mbox{corr}\left[\zeta_{n}^{y}\left(\theta_{t}^{y}\right),\mathcal{V}\left(\theta_{t}^{y}\right)\right]\right)$, where $\mbox{corr}\left[\mathbf{a}, \mathbf{b}\right]$ denotes the correlation coefficient between $\mathbf{a}$ and $\mathbf{b}$ and given by

\begin{equation}
\mbox{corr}\left[\mathbf{a}, \mathbf{b}\right]=\frac{\mbox{cov}\left[\mathbf{a}, \mathbf{b}\right]}{\sqrt{\mbox{var}\left[\mathbf{a}\right]}~\sqrt{\mbox{var}\left[\mathbf{b}\right]}}.
\end{equation}

To select an appropriate CVF we need to consider two criteria. First, its expectation with respect to $\Phi^{*}$ should have an exact evaluation. Second, it has to be correlated with the estimated error. A favorable candidate is the analytic true error of linear classifiers. We consider in this study a CVF given by the true error of an LDA classifier defined by
$g_{N_{t}}\left(\mathbf{x}\right)=\mathbf{a}_{N_{t}}^{T}\mathbf{x}+b_{N_{t}}$ where $\mathbf{a}_{N_{t}}=\mathbf{S}_{t}^{-1}\left(\overline{\mathbf{x}}_{t}^{1}-\overline{\mathbf{x}}_{t}^{0}\right)$, $b_{N_{t}}=-\frac{1}{2}\mathbf{a}^{T}\left(\overline{\mathbf{x}}_{t}^{1}+\overline{\mathbf{x}}_{t}^{0}\right)+\ln{\frac{n_{t}^{1}}{n_{t}^{0}}}$ and the pooled covariance $\mathbf{S}_{t}$ is given by

\begin{equation}
\mathbf{S}_{t}=\frac{\left(n_{t}^{0}-1\right)\mathbf{S}_{t}^{0}+\left(n_{t}^{1}-1\right)\mathbf{S}_{t}^{1}}{N_{t}-2}.
\end{equation}

$\overline{\mathbf{x}}_{t}^{y}$ and $\mathbf{S}_{t}^{y}$ are the empirical estimates utilized in (<ref>). Thus, the CVF is given by

\begin{equation}
\mathcal{V}\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)=\mathbf{\Phi}\left(\frac{\left(-1\right)^{y}g_{N_{t}}\left(\mu_{t}^{y}\right)}{\sqrt{\mathbf{a}_{N_{t}}^{T}~\left(\mathbf{\Lambda}_{t}^{y}\right)^{-1}~\mathbf{a}_{N_{t}}}}\right)
\end{equation}

with $\mathbf{\Phi}$ denoting the standard normal Gaussian cumulative distribution function (CDF).
Now it remains only to determine $E_{\Phi^{*}}\left[\mathcal{V}\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)\right]$ in closed-form to fully define the estimation set-up. We can show after simplifications and using results from [17] that

\begin{align}
\begin{split}
&E_{\Phi^{*}}\left[\mathcal{V}\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)\right]=\frac{1}{2}\\
\end{split}
\end{align}

where $\mbox{sgn}\left(\cdot\right)$ is the sign function,

\begin{equation}
\end{equation}

and $\mathcal{I}\left(\cdot;\cdot,\cdot\right)$ denotes the regularized incomplete beta function given by

\begin{equation}
\mathcal{I}\left(x;a,b\right)=\frac{\Gamma\left(a+b\right)}{\Gamma\left(a\right)\Gamma\left(b\right)}\int_{0}^{x}t^{a-1}\left(1-t\right)^{b-1}dt,
\end{equation}

with $\Gamma\left(\cdot\right)$ being the regular univariate gamma function. Details for simplifying $E_{\Phi^{*}}\left[\mathcal{V}\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)\right]$ are covered in Appendix D.

The complete specification of the CVF concludes our importance sampling set-up. We enumerate some advantages of the proposed set-up over direct sampling methods. First, the importance density $\Phi^{*}$ is much simpler than the nominal density $\pi^{*}$ that involves matrix-variate hypergeometric functions. Second, our set-up successfully combines two variance reduction techniques that enable accurate estimation. Last, and most importantly, the independence of the generated Monte Carlo samples w.r.t source data permits the reuse of the sampled parameters with various source datasets for fixed models. This reusability significantly reduces the computational cost of sampling from $\Phi^{*}$ and makes the utilization of advanced MCMC methods amenable as the whole process could be accelerated by a factor of 10$\sim$20, that also scales up with the dimensionality and the number of used source datasets. For efficient sampling from $\Phi^{*}$, we use Hamiltonian Monte Carlo (HMC), proven to have a superior performance to standard MCMC samplers [38]. For this purpose, we utilize the STAN software that offers a full Bayesian statistical inference framework with HMC [38].

§ EXPERIMENTS AND DATASETS

To evaluate the performance of the proposed error estimator, we consider the mean-square error (MSE) as a performance measure to understand the joint behaviour of the classification error $\varepsilon_{n}$ and its estimate $\hat{\varepsilon}$. For the random vector $\left(\varepsilon_{n}, \hat{\varepsilon}\right)$, the MSE is defined as

\begin{equation}
\mbox{MSE}\left(\hat{\varepsilon}\right)=\mbox{E}\left[\left|\hat{\varepsilon}-\varepsilon_{n}\right|^{2}\right].
\end{equation}

In what follows,  we present the details of the experimental set-up for demonstrating the performance of the proposed TL-based BEE based on three different types of classifiers applied to both synthetic data as well as real biological datasets.

§.§ Synthetic datasets

To simulate and verify the extent of knowledge transferability across domains, we consider a wide range of joint prior densities that model the different levels of relatedness between the source and target domains. The proposed set-up is as follows. We consider a binary classification problem in the context of homogeneous TL with dimensions 2, 3, and 5. In the simulated datasets, the number of source data points per class varies between 10 and 500 and between 5 and 50 for target datasets. This mimics realistic settings of small-size sample conditions (especially, in the target domain) as reported in the literature [3]. We set up the data distributions as follows. $\nu=\nu^{y}=d+20$, $\kappa_{t}=\kappa_{t}^{y}=100$, $\kappa_{s}=\kappa_{s}^{y}=100$, $\mathbf{m}_{t}^{0}=\mathbf{0}_{d}$, $\mathbf{m}_{t}^{1}=\vartheta \times\mathbf{1}_{d}$, $\mathbf{m}_{s}^{0}=\mathbf{m}_{t}^{0}+10\times\mathbf{1}_{d}$, $\mathbf{m}_{s}^{1}=\mathbf{m}_{t}^{1}+10\times\mathbf{1}_{d}$, where $\vartheta$ is an adjustable scalar used to control the Bayes error in the target domain, and $\mathbf{0}_{d}$ and $\mathbf{1}_{d}$ are $d\times 1$ all-zero and all-one vectors, respectively. For the scale matrices of Wishart distributions we set $\mathbf{M}_{t}^{y}=k_{t}\mathbf{I}_{d}$, $\mathbf{M}_{s}^{y}=k_{s}\mathbf{I}_{d}$, and $\mathbf{M}_{ts}^{y}=k_{ts}\mathbf{I}_{d}$ where $\mathbf{I}_{d}$ is the identity matrix of rank $d$.

To ensure that the joint scale matrix $\mathbf{M}^{y}=\begin{pmatrix}
\mathbf{M}_{t}^{y} & \mathbf{M}_{ts}^{y}\\
{\mathbf{M}_{ts}^{y}}^{T} & \mathbf{M}_{s}^{y}
\end{pmatrix}$ is positive definite $\forall y\in\left\{0,1\right\}$, we set $k_{ts}=\alpha \sqrt{k_{t}k_{s}}$ with $k_{t}>0$, $k_{s}>0$, and $\left|\alpha\right|<1$. As in [23], the value of $\left|\alpha\right|$ controls the amount of relatedness between the source and target domains.

To control the level of relatedness by adjusting only $\left|\alpha\right|$ without involving other confounding factors, we set $k_{t}=k_{s}=1$.
To sample from the joint prior, we first sample from a non-singular Wishart distribution $W_{2d}\left(\mathbf{M}^{y}, \nu\right)$ to get a block partitioned sample of the form $\mathbf{\Lambda}^{y}=\begin{pmatrix}
\mathbf{\Lambda}_{t}^{y} & \mathbf{\Lambda}_{ts}^{y}\\
{\mathbf{\Lambda}_{ts}^{y}}^{T} & \mathbf{\Lambda}_{s}^{y}
\end{pmatrix}$ from which we extract $\left(\mathbf{\Lambda}_{t}^{y}, \mathbf{\Lambda}_{s}^{y}\right)$. Afterwards, we sample $\mu_{z}^{y} \sim \mathcal{N}\left(\mathbf{m}_{z}^{y}, \left(\kappa_{z}^{y}\mathbf{\Lambda}_{z}^{y}\right)^{-1}\right)$ for $z\in \left\{s,t\right\}$ and $y\in \left\{0,1\right\}$. We use in our simulations two types of datasets. Training datasets that contain samples from both domains and testing datasets that contain only samples from the target domain. In all the simulations we consider testing datasets of 1,000 data points per class and we assume equal prior probabilities for the classes.

§.§ RNA sequencing (RNA-seq) datasets

To evaluate the performance of the TL-based BEE on real-world data, we consider classifying patients diagnosed with Schizophrenia using transcriptomic profiles collected from psychiatric disorder studies [39]. Based on two RNA-seq datasets listed in Table <ref>, we selected the transcriptomic profiles of 3 genes, based on a stringent feature selection procedure comprising the analysis of differential gene expression, clustering of gene-gene interactions, and statistical testing for multivariate normality. More specifically, we focus on analyzing the astrocyte-related cluster of differentiation 4 (CD4), found to be significantly up-regulated in subjects with Schizophrenia [39]. We select the top three hub genes that collectively satisfy the Royston's multivariate normality test applied to the full datasets for both classes at a significance level of 99%. The identified genes satisfying all the aforementioned criteria include SOX9, AHCYL1, and CLDN10 with an average module centrality of $0.86$ measured by genes' module membership (kME) [39].
In addition to normalization and quality control performed in [39], the selected features in both datasets have been further standardized to zero means and unit variances across both classes as in [24, 23].

2*Disease     2c|Number of samples     1c|2*Brain region     1c|2*Dataset
Case     Control     1c|     1c|
2*Schizophrenia     53     53     1c|Frontal cortex     1c|syn4590909 [39]
262     293     1c|DLPFC     1c|syn2759792 [40]
TOTAL     315     346

Independent datasets sampled from two different brain tissues.

We consider the dataset syn2759792, sampled from the brain dorsolateral prefrontal cortex area (DLPFC), as a target dataset and syn4590909, sampled from the frontal cortex region (FC), as a source dataset. Among 555 postmortem brain samples in syn2759792, we randomly draw 5 samples per class as training data and we use the remaining samples to evaluate the classification error.
This process is repeated 10,000 times to estimate the average MSE deviation of the TL-based BEE from the true error.
To tune the model hyperparameters, we assume shared values for case and control samples in source and target domains. We set $\nu=10d=30$, $n_{t}=5$, and we consider source datasets of different sizes ($n_{s}\in\left\{10,~30,~50\right\}$) to conduct a greedy search to identify the optimal relatedness coefficient $\left|\alpha\right|$ across domains. At each iteration, we randomly permute the source samples for statistical significance. The remaining parameters are set as follows: $\kappa_{t}=n_{t}$, $\kappa_{s}=n_{s}$, and $k_{t}=k_{s}=\frac{1}{\nu}$ such that the mean of the Wishart precision matrices will be equal to the identity matrix which matches the normal standardization. For mean vectors $\mathbf{m}_{t}$ and $\mathbf{m}_{s}$, we pool all case and control samples in each domain and consider their means, respectively.

§.§ Classifier design

For comprehensive evaluation of our novel TL-based error estimator, we design and perform a set of experiments. The proposed TL-based estimator is applied to a collection of classifiers with different levels of learning capacities and tested under various scenarios.
To separate error estimation from classifier design, we start by analyzing the performance of the TL-based BEE estimator for fixed classifiers that do not depend on training data. This set-up distinctly reveals the major characteristics of the TL-based BEE, excluding any confounding factors that may stem from classifier design and the performance of the resulting classifier.

Next, we also conduct a comparative study of the TL-based BEE performance with respect to other widely-used error estimators, which include resubstitution, CV, LOO, and the 0.632-bootstrap estimators. As these popular data-driven estimators involve  classifier design on the training data, we will also consider a TL-based classifier designed on target and source data that operates in the target domain for comparison. For this, we employ the optimal Bayesian transfer learning (OBTL) classifier introduced in [23], which shares the same Bayesian framework on which our TL-based BEE is developed. In what follows, we recall the definition of each classifier considered in our evaluations and also present the details of the evaluation experiments performed in this study.

In the first set of experiments, we employ a fixed quadratic classifier assuming we know beforehand the true target parameters. For normally distributed data, this quadratic classifier corresponds also to the Bayes classifier that is optimal for the given feature-label distributions. Using quadratic discriminant analysis (QDA), we define

$\Psi_{QDA}\left(\mathbf{x}\right)=\mathbf{x}^{T}\mathbf{A}\mathbf{x}+\mathbf{b}^{T}\mathbf{x}+c$, where

\begin{align}
\begin{split}
&\mathbf{A}=-\frac{1}{2}\left(\mathbf{\Lambda_{t}^{1}}-\mathbf{\Lambda_{t}^{0}}\right), \quad%\\
\mathbf{b}=\mathbf{\Lambda_{t}^{1}}\mu_{t}^{1}-\mathbf{\Lambda_{t}^{0}}\mu_{t}^{0},\\
\end{split}
\end{align}

The error estimation problem turns out to be an estimation of the Bayes error that coincides here with the true error of the designed QDA. Obviously, this classifier is independent from any observed sample as it is fixed assuming known true model parameters. Without loss of generality, we apply the TL-based BEE using labeled observations from a compound dataset compiled from target and source domains.

In the second set of experiments we investigate the behaviour of the TL-based BEE within the class of sub-optimal classifiers. To this end, we consider a linear classifier derived through linear discriminant analysis (LDA) and we define
$\Psi_{LDA}\left(\mathbf{x}\right)=\mathbf{a}^{T}\mathbf{x}+b$ where $\mathbf{a}=\mathbf{S}_{t}^{-1}\left(\mu_{t}^{1}-\mu_{t}^{0}\right)$, $b=-\frac{1}{2}~\mathbf{a}^{T}\left(\mu_{t}^{1}+\mu_{t}^{0}\right)$, and the average covariance $\mathbf{S}_{t}$ is given by

\begin{equation}
\mathbf{S}_{t}=\frac{\left(\mathbf{\Lambda}_{t}^{0}\right)^{-1} +\left(\mathbf{\Lambda}_{t}^{1}\right)^{-1} }{2}.
\end{equation}

Our goal is then to approximate the true error of this sub-optimal classifier using transfer learning.

Next, we evaluate the performance of the TL-based BEE for the OBTL classifier that can take advantage of both source and target domain data. The OBTL classifier is defined by

\begin{equation}
\Psi_{OBTL}\left(\mathbf{x}\right) = \operatorname*{arg~max}_{y \in \left\{0,1\right\}} \mathcal{O}_{OBTL}\left(\mathbf{x}|y\right),
\end{equation}

where the objective function $\mathcal{O}_{OBTL}\left(\mathbf{x}|y\right)$ denotes the effective class-conditional density $p\left(\mathbf{x}|y\right)$ given by the following theorem:

Theorem 2 [23]:
The effective class-conditional density, denoted by $p\left(\mathbf{x}|y\right) = \mathcal{O}_{OBTL}\left(\mathbf{x}|y\right)$, in the target domain is given by

\begin{align}
\begin{split}
&\times \Gamma_{d}^{-1}\left(\frac{\nu^{y}+n_{t}^{y}}{2}\right)\left|\mathbf{T}_{x}^{y}\right|^{\frac{\nu^{y}+n_{t}^{y}+1}{2}}\left|\mathbf{T}_{t}^{y}\right|^{-\frac{\nu^{y}+n_{t}^{y}}{2}}\\
&\times \pFq{2}{1}{\frac{\nu^{y}+n_{s}^{y}}{2}, \frac{\nu^{y}+n_{t}^{y}+1}{2}}{\frac{\nu^{y}}{2}}{\mathbf{T}_{s}^{y}\mathbf{F}^{y}\mathbf{T}_{x}^{y}{\mathbf{F}^{y}}^{T}}\\
&\times \pFq{2}{1}{\frac{\nu^{y}+n_{s}^{y}}{2}, \frac{\nu^{y}+n_{t}^{y}}{2}}{\frac{\nu^{y}}{2}}{\mathbf{T}_{s}^{y}\mathbf{F}^{y}\mathbf{T}_{t}^{y}{\mathbf{F}^{y}}^{T}}^{-1},
\end{split}
\end{align}

\begin{align}
\begin{split}
\kappa_{x}^{y} &= \kappa_{t,n}^{y}+1 = \kappa_{t}^{y}+n_{t}^{y}+1,\\
\left(\mathbf{T}_{x}^{y}\right)^{-1} &= \left(\mathbf{T}_{t}^{y}\right)^{-1} + \frac{\kappa_{t,n}^{y}}{\kappa_{t,n}^{y}+1}\left(\mathbf{m}_{t,n}^{y}-\mathbf{x}\right)\left(\mathbf{m}_{t,n}^{y}-\mathbf{x}\right)^{T}.
\end{split}
\end{align}

§.§ Simulation set-up

Figure <ref> provides a combined illustration of the simulation set-up for all three classifiers.
For rigorous evaluation of the performance of the proposed TL-based BEE, we primarily focus our experiments on assessing the impact of using different types and amounts of source data. This is enabled by the joint prior imposed over the model parameters and controlled by the relatedness coefficient $\left|\alpha\right|$ that dictates the extent of interaction between the features in the two domains. For this purpose, we repeatedly conduct experiments following the flow chart in Fig. <ref> with different relatedness values ($\left|\alpha\right| = \left[0.1, 0.3, 0.5, 0.7, 0.9, 0.95\right]$) where $\left|\alpha\right|=0.1$ corresponds to the lowest relatedness bewteen the two domains and $\left|\alpha\right|=0.95$ reflects the highest relatedness within the range of studied values.

startstop = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, align=center, draw=black, fill=green!10]
figure = [rectangle, rounded corners, minimum width=1cm, minimum height=1cm, align=center, draw=black, fill=blue!10]
io = [trapezium, trapezium left angle=70, trapezium right angle=110, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=blue!30]
process = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=orange!30]
decision = [diamond, rounded corners, minimum width=0cm, minimum height=0cm, text centered, draw=black, fill=yellow!10]
parameter = [ellipse, minimum width=1cm, minimum height=1cm, text centered, draw=black, fill=red!10]
data = [cloud, minimum width=1cm, minimum height=1cm, text centered, draw=black, fill=blue!10]
arrow = [thick,->,>=stealth]
line = [draw, -latex']
empty = [ellipse, minimum width=0cm, minimum height=0cm, text centered, draw=white, fill=white!0]

[node distance=1.5 cm, scale=0.8]
(L_0_N_0) [startstop]

$d$, $Y$, $\mathbf{m}_{z}$, $\kappa_{z}$, $\nu$, $\mathbf{M}$, $\left|\alpha\right|$

$N_{p}$, $N_{d}$, $N_{t}$, $N_{s}$, $\vartheta$, $\tau$

(L_1_N_0) [parameter, below of=L_0_N_0] $\Lambda^{y}$;

(L_2_N_0) [parameter, below of=L_1_N_0, xshift=3.5cm]$\mu_{s}^{y}$;
(L_2_N_1) [parameter, left of=L_2_N_0, xshift=-0.5cm]$\Lambda_{s}^{y}$;
(L_2_N_2) [parameter, left of=L_2_N_1, xshift=-1.5cm]$\Lambda_{t}^{y}$;
(L_2_N_3) [parameter, left of=L_2_N_2, xshift=-0.5cm]$\mu_{t}^{y}$;

[on background layer]
(B_L_2_N_0)[fit=(L_2_N_0)(L_2_N_1), draw, dotted, thick, blue, fill=black!3];

[on background layer]
(B_L_2_N_1) [fit=(L_2_N_2)(L_2_N_3), draw, dotted, thick, blue, fill=black!3];

(L_3_N_0)[data, below of=L_2_N_0]$\mathcal{D}_{s}^{y}$;
(L_3_N_1)[data, below of=L_2_N_1]$\mathcal{D}_{t}^{y}$;
(L_3_N_2)[data, below of=B_L_2_N_1]$\mathcal{D}_{t}^{y}$;

[on background layer]
(B_L_3_N_0) [fit=(L_3_N_0)(L_3_N_1), draw, dotted, thick, blue, fill=black!3];

(L_4_N_0)[process, align=center, below of=B_L_3_N_0]Derive
$\Psi_{OBTL}\left ( \mu_{t}, \Lambda_{t}\right )$;
(L_4_N_1)[process, align=center, below of=L_3_N_2, xshift=0.5cm]Derive
$\Psi_{LDA}\left ( \mu_{t}, \Lambda_{t}\right )$;
(L_4_N_2)[process, align=center, left of=L_4_N_1, xshift=-2.5cm]Derive
$\Psi_{QDA}\left ( \mu_{t}, \Lambda_{t}\right )$;

(L_5_N_0)[process, align=center, below of=L_4_N_0]Estimate
(L_5_N_1)[process, align=center, below of=L_4_N_1]Estimate
(L_5_N_2)[process, align=center, below of=L_4_N_2]Estimate

[on background layer]
(B_L_4_N_0) [fit=(L_4_N_0)(L_5_N_0), draw, dotted, thick, blue, fill=black!3];

[on background layer]
(B_L_4_N_1) [fit=(L_4_N_1)(L_5_N_1), draw, dotted, thick, blue, fill=black!3];

[on background layer]
(B_L_4_N_2) [fit=(L_4_N_2)(L_5_N_2), draw, dotted, thick, blue, fill=black!3];

(L_6_N_0)[process, align=center, below of=L_5_N_0]Draw
$\left( \mu_{t}^{y}, \Lambda_{t}^{y}\right )_{i} \sim \pi^{*} \left (\mu_{t}^{y}, \Lambda_{t}^{y} | \mathcal{D}_{t}^{y},\mathcal{D}_{s}^{y}\right )$;

(L_7_N_0)[process, below of=L_6_N_0]
$BEE\left ( \varepsilon_{\Psi}^{y} \right ) =E_{\pi^{*}}\left [ \varepsilon_{\Psi}^{y}  \right ]$

(L_8_N_0)[process, below of=L_7_N_0]$SE\left ( BEE\right )=\left ( BEE\left ( \varepsilon_{\Psi}\right )-\varepsilon_{\Psi}\right )^{2}$;

(L_9_N_0)[decision, below of=L_5_N_2, yshift=-1cm]$\varepsilon_{Bayes}=\tau$?;

(L_10_N_0)[decision, below of=L_8_N_0, yshift=-1cm]$N_{d}$ repetitions?;

(L_11_N_0)[process, below of=L_10_N_0, yshift=-1cm]MSE = $\frac{1}{N_{d}}\sum_{1}^{N_{d}}\left ( BEE\left ( \varepsilon_{\Psi}\right )-\varepsilon_{\Psi}\right )_{1..N_{d}}^{2}$;

(L_12_N_0)[process, below of=L_11_N_0, align=center]Repeat for $N_{s}$
$\left [ 10, 50:50:500\right ]$;

(L_13_N_0)[process, below of=L_12_N_0, align=center]Repeat for $N_{t}$
$\left [ 5, 10:10:50\right ]$;

(L_14_N_0)[decision, below of=L_13_N_0, yshift=-1cm]$N_{p}$ repetitions?;

(L_14_N_1)[process, left of=L_14_N_0, xshift=-6cm]$MSE_{Avg}$ = $\frac{1}{N_{p}}\sum_{1}^{N_{p}}\left ( MSE\left ( \varepsilon_{\Psi}\right )\right )_{1..N_{p}}$;

(L_15_N_0)[startstop, left of=L_11_N_0, align=center, xshift=-6cm]Evaluate MSE
For each $n_{t}\in N_{t}$

Evaluate $MSE_{Avg}\left ( BEE \left( \varepsilon_{\Psi}\right )\right )=f\left (n_{s}\in N_{s} \right )$

(Display $n_{t}~\&~\varepsilon_{Bayes}$)

For each $n_{s}\in N_{s}$

Evaluate $MSE_{Avg}\left ( BEE \left( \varepsilon_{\Psi}\right )\right )=f\left (n_{t}\in N_{t} \right )$

(Display $n_{s}~\&~\varepsilon_{Bayes}$);

(L_16_N_0)[figure, above of= L_15_N_0, xshift=2cm, yshift=2.5cm]
[scale=0.6, transform shape]
[->] (0.7, 0) – (3, 0) node[below] $n_s$;
[->] (0.7, 0) – (0.7, 2) node[above] $MSE$;
[scale=1, domain=0.8:2.5, smooth, variable=, red, very thick] plot (, 1/(*));
[step=0.3,gray, very thin] (0.8,0) grid (2.7,1.7);

(L_16_N_1)[figure, above of= L_15_N_0, xshift=-2cm, yshift=2.5cm]
[scale=0.6, transform shape]
[->] (0.7, 0) – (3, 0) node[below] $n_s$;
[->] (0.7, 0) – (0.7, 2) node[above] $MSE$;
[scale=1, domain=0.8:2.5, smooth, variable=, blue, very thick] plot (, 1/(*));
[step=0.3,gray, very thin] (0.8,0) grid (2.7,1.7);

[arrow] (L_1_N_0) – (L_2_N_1);
[arrow] (L_1_N_0) – (L_2_N_2);

[arrow] (L_2_N_1) – (L_2_N_0);
[arrow] (L_2_N_2) – (L_2_N_3);
[arrow] (L_2_N_2) – (L_3_N_2);
[arrow] (L_2_N_3) – (L_3_N_2);
[arrow] (L_2_N_0) – (B_L_3_N_0);
[arrow] (L_2_N_1) – (B_L_3_N_0);

[arrow] (B_L_3_N_0) – (L_4_N_0);
[arrow] ($(L_3_N_2.west)$) -| ($(L_5_N_2.east)+(0.65,0)$) – ($(L_5_N_1.west)$);
[arrow] ($(L_5_N_2.east)+(0.65,0)$) – ($(L_5_N_2.east)$);

[arrow] (L_4_N_0) – (L_5_N_0);

[arrow] ($(B_L_3_N_0.east)+(0,-0.5)$) -| ($(L_6_N_0.east)+(0.5,0)$) – (L_6_N_0);

[arrow] (L_5_N_2) – (L_9_N_0);
[arrow] (L_6_N_0) – (L_7_N_0);
[arrow] (L_7_N_0) – (L_8_N_0);
[arrow] (L_8_N_0) – (L_10_N_0);
[arrow] (L_11_N_0) – (L_12_N_0);
[arrow] (L_12_N_0) – (L_13_N_0);
[arrow] (L_13_N_0) – (L_14_N_0);

[arrow] (L_4_N_1) – (L_5_N_1);
[arrow] (L_4_N_2) – (L_5_N_2);

[arrow] ($(L_9_N_0.east)$) node[anchor=south west]Yes $\rightarrow$ save $\left \{\left ( \mu_{t}^{y}, \Lambda_{t}^{y}\right )~\&~\Psi \right \}$ -| ($(B_L_3_N_0.west)+(-1.3,0)$) – (B_L_3_N_0);

[arrow] ($(L_9_N_0.west)$) node[anchor=south east]No-| ($(B_L_2_N_1.west)+(-4.8,0)$) –node[anchor=south]Update $\mathbf{m}_{t}^{y}\left ( \vartheta\right )$ (B_L_2_N_1);

[arrow] ($(L_10_N_0.east)$)  node[anchor=south west]No -| ($(B_L_3_N_0.east)+(2,0)$) – (B_L_3_N_0);

[arrow] ($(L_10_N_0.south)$) – node[anchor=west]Yes (L_11_N_0);

[arrow] ($(L_12_N_0.west)$) -| ($(L_13_N_0.west)+(-0.5,0)$) node[anchor=north east, align=center]$MSE_{n_{t}}\left ( \varepsilon_{\Psi}\right )=f\left ( n_{s}\right )$:

$N_{t}$ loop outer

$N_{s}$ loop inner– ($(L_13_N_0.west)$);

[arrow] ($(L_13_N_0.east)$) node[anchor=north west, align=center]Otherwise:

$N_{t}$ loop inner

$N_{s}$ loop outer-|($(L_12_N_0.east)+(0.5,0)$) – ($(L_12_N_0.east)$);

[arrow] (L_14_N_1) – (L_15_N_0);
[arrow] ($(L_14_N_0.west)$) –node[anchor=south west]Yes (L_14_N_1);
[arrow] ($(L_14_N_0.east)$) node[anchor=south west]No-| ($(L_1_N_0.east)+(7.2,0)$) – ($(L_1_N_0.east)$);

[arrow] ($(B_L_2_N_1.west)+(0,-0.5)$) -| (L_4_N_2);
[arrow] ($(B_L_2_N_1.south)+(1.5,0)$) |- ($(L_4_N_1.north)+(0,0.5)$) -| ($(L_4_N_1.north)$);

[arrow] ($(L_15_N_0.north)+(2.5,0)$) -| (L_16_N_0);
[arrow] ($(L_15_N_0.north)+(-2.5,0)$) -| (L_16_N_1);

Flow chart illustrating the simulation set-up based on synthetic datasets.

In the first set of experiments we start by drawing a joint sample $\left(\mathbf{\Lambda}_{t}^{y},\mathbf{\Lambda}_{s}^{y}\right)$ for each class $y\in \left\{0,1\right\}$ as previously described. Next, we iterate over the values of the hyperparameter $\vartheta$ to control $\mathbf{m}_{t}\left(\vartheta\right)$ through a dichotomic search to get a desired value $\tau$ of the Bayes error. This is achieved by drawing a sample $\mu_{t}^{y}\sim \mathcal{N}\left(\mathbf{m}_{t}\left(\vartheta\right), \left(\kappa_{t}^{y}\mathbf{\Lambda}_{t}^{y}\right)^{-1}\right)$ and then generating a test set based on the joint sample $\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)$. Using this test set, we determine the true error of the optimal QDA derived from $\left(\mu_{t}^{y}, \mathbf{\Lambda}_{t}^{y}\right)$. If the desired Bayes error (true error of the designed QDA) is attained then the iteration stops, otherwise we update $\vartheta$ and reiterate.
In our experiments, we set $\tau=0.2$ to mimic a moderate level of classification complexity. This step is indeed crucial as it maintains the same level of complexity across the experiments and guarantees a fair comparison across different levels of relatedness. We note that this procedure is valid for general covariances as it acts only on updating the value of the mean parameter without altering the structure of the covariances nor the random mean vectors. Obviously, this approach to specify the Bayes error maintains the Bayesian transfer learning framework intact. However, it is not guaranteed to find values of target parameters that correspond to the desired Bayes error especially for high dimensions and complex classification (large Bayes error) as we will discuss in the next section. Once the problem complexity is set and the classifier is fixed, we generate $N_d=10,000$ training datasets that we use to evaluate the MSE of the TL-based BEE as depicted in Fig. <ref>.
To estimate the TL-based BEE, we employ the importance sampling set-up previously described and we draw 1,000 MC samples from the importance density using HMC sampler.

In the second set of experiments, we follow a similar set-up using an LDA classifier designed based on the true model parameters. As before, we employ QDA to determine the Bayes error to maintain the same complexity level across different experiments. As in the first set of experiments, we use the TL-based BEE to estimate the true error of the designed LDA classifier.

In the last set of experiments on synthetic datasets, we conduct a comparative analysis study using an OBTL classifier designed using training datasets generated from the model parameters specified by the Bayes error. The error estimation task, in this scenario, aims at approximating the true error of the designed OBTL classifier determined using a large test set generated from the true feature-label distributions. As illustrated in Fig. <ref>, QDA and LDA classifiers are fixed and derived from the true model parameters while the OBTL classifier is designed based on training datasets collected from the underlying feature-label distributions that correspond to the specified Bayes error. In all simulations, the designed classifiers are fixed given the observed samples and the TL-based BEE estimator is safely applied.
Finally, regarding synthetic datasets, we note that the flow chart in Fig. <ref> is valid for all classifiers (QDA, LDA, and OBTL) and the notation $\Psi$ designates the classifier of interest in the corresponding set of experiments. For instance, in the second set of experiments, $\Psi$ refers to $\Psi_{LDA}$.

In addition to this in-depth analysis of the performance, behavior, and characteristics of our proposed TL-based BEE based on synthetic datasets, we also performed additional validation based on real-world biological datasets. By using RNA-seq datasets syn2759792 and syn4590909  taken from different brain regions for studying brain disorders, we train a QDA classifier using the  target data from the RNA-seq dataset syn2759792, and we leverage the source data from syn4590909 to evaluate the performance of the proposed TL-based BEE.

In the following section, we present the experimental results with in-depth discussions.

§ RESULTS

§.§ Performance on synthetic datasets

We start by evaluating the performance of the proposed TL-based BEE in estimating the Bayes error, that corresponds to the true error of the QDA in the target domain, for different levels of $\left|\alpha\right|$ and different size combinations of the utilized source and target datasets.

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,scaled y ticks=base 10:4,

legend columns=3,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named1,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,ymax=6.5e-4,
yticklabel style=
/pgf/number format/fixed,
/pgf/number format/fixed zerofill,
/pgf/number format/precision=1
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,scaled y ticks=base 10:4,

xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,
xlabel near ticks,
yticklabel style=
/pgf/number format/fixed,
/pgf/number format/fixed zerofill,
/pgf/number format/precision=1
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:3,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=20$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:3,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=50$

MSE deviation from true error for Gaussian distributions with respect to source sample size. The Bayes error is fixed at 0.2 in all figures.

In Fig. <ref> we investigate the behavior of the TL-based BEE when the target data is fixed while we vary the size of the source data. We show the results for $d=2$ in the top row, the results for $d=3$ in the second row, and the results for $d=5$ in the last row. The columns correspond to the results for target datasets with different sizes: $n_{t}=20$ on the left and $n_{t}=50$ on the right. The MSE curves show similar trends for all three values of $d$, where we can see that the deviation of the error estimate from the true error significantly decreases when highly related source data are employed. This behavior diminishes as the relatedness between the two domains decreases.
Notably, using large source datasets ($n_{s}\geq 200$) of moderate to small relatedness values ($\left|\alpha\right|\leq 0.7$) does not negatively impact the performance of the estimator for low dimensions ($d\in \left\{2,3\right\}$) as shown in rows 1 and 2 of Fig. <ref>. As the dimensionality further increases ($d=5$), large source datasets of moderate to small relatedness slightly increase the deviation of the estimated error from the true error (i.e. $\left|\alpha\right|=0.7$ in the third row). This tiny asymptotic deviation is explained by potential undesirable effects of relying on large source datasets of modest relatedness. However, it is important to note that the proposed TL-based BEE in the context of the given Bayesian TL framework suppresses this behaviour,  as it does not directly depend on the source data but the information transfer occurs through the joint prior. The joint prior acts like a bridge through which the useful knowledge passes from the source to the target domain. Effects of using source data in different TL settings (especially, a non-Bayesian setting) may require further investigation.

Moreover, the simulation results in different columns show that the MSE deviation decreases as we rely on larger target datasets. However, the gain in performance as we use additional source data is reduced when target data are more abundant. This is illustrated by the slope of the MSE graphs that flattens as $n_{t}$ increases.

Finally, Fig. <ref> shows that for higher dimensions, the MSE deviation tends to increase. This is expected as increasing the dimensionality generally leads to a more difficult error estimation problem.

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, ylabel near ticks,
xlabel near ticks,
xmin=5,xmax=50,xtick=5,10,20,30,40,50, scaled y ticks=base 10:4,

legend columns=3,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named2,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_s=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,  scaled y ticks=base 10:4,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_s=200$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, ylabel near ticks,
xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_s=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,
xlabel near ticks,

legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_s=200$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_s=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE, xmin=5,xmax=50,xtick=5,10,20,30,40,50,ymin=0,ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_s=200$

MSE deviation from true error for Gaussian distributions with respect to target sample size. The Bayes error is fixed at 0.2 in all figures.

Next, Fig. <ref> shows the MSE deviation with respect to the size of the target dataset for dimensions 2, 3, and 5.
The first column corresponds to the case of using source datasets of size $n_{s}=50$ and the second column shows the results for $n_{s}=200$. The  performance of the TL-based BEE estimator improves with the increasing availability of target data. We can also clearly see that the MSE deviation from the true error asymptotically converges to comparable values for all relatedness levels. When highly related source data are available, the TL-based estimator yields accurate estimation results even when the target dataset is small.
These results consolidate the findings in Fig. <ref> about the redundancy of source data in the presence of abundant target data. Across all graphs in Fig. <ref>, we can see that a relatedness coefficient $\left|\alpha\right|=0.95$ results in a nearly constant deviation from the true error as a function of target data size, which suggests that highly related source data $\left|\alpha\right|>0.95$ acts almost identically like the target data, regardless of the shift across the domains in terms of their means.

Similar to the trends shown in Fig. <ref>, results across different rows of Fig. <ref> demonstrate that the error estimation difficulty increases with the increase of dimensionality. This is clearly reflected in the MSE deviation from the true error in Fig. <ref>, which shows that as the dimension increases from $d=2$ (top row) to $d=5$ (last row), the MSE increases by one order of magnitude.

Now, we aim at investigating the effect of classification complexity on the performance of the proposed TL-based BEE. To this end, we conduct simulations, in which we vary the Bayes error through a wide range of possible values and evaluate the TL-based BEE at each given Bayes error for different sizes of target data while using source datasets of a fixed size $n_{s}=200$. In binary classification, the Bayes error has an upper bound specified by the true error of random classification that is $0.5$, as every data point can be randomly assigned one of the class labels. Ideally, we would vary the Bayes error across the interval $\left[0,~0.5\right]$ as in [17]. However, in our set-up, we do not impose any structure on the covariance matrices, nor do we assume they are scaled identities. This makes the control of the Bayes error much more difficult. Additionally, the joint sampling set-up within our Bayesian transfer learning framework inhibits any modification of the randomized parameters. Consequently, the only practical way to adjust the Bayes error is to tune the mean vector parameters $\mathbf{m}_{t}^{y}$ that specify the means for the class mean vectors $\mu_{t}^{y}$ with $y\in \left\{0,1\right\}$.
In our experiments, we were able to fully control the Bayes error for $d=2$ and we considered the following values $\left[0.05,~0.1,~0.15,~0.2,~0.25,~0.3,~0.35,~0.4,~0.45,~0.5\right]$. Achieving the same range of values for $d=3$ and $d=5$ was more challenging, and our implemented heuristic did not converge for high values of Bayes error as setting $\mathbf{m}_{t}^{0}=\mathbf{m}_{t}^{1}$ did not help in increasing the Bayes error. However, we were able to vary the Bayes error for $d=3$ within the range $\left[0.05,~0.1,~0.15,~0.2,~0.25,~0.3,~0.35,~0.4,~0.45\right]$, and for $d=5$, within $\left[0.05,~0.1,~0.15,~0.2,~0.25,~0.3,~0.35,~0.4\right]$, sufficient for observing the trends.

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
scaled y ticks=base 10:4,
legend columns=3,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named3,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=5.000000e-01,ymin=0,ymax=1.2e-03, scaled y ticks=base 10:4,
ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error,xtick=0.1,0.2,0.3,0.4, ylabel= MSE,ylabel near ticks, scaled y ticks=base 10:4,
xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.500000e-01,ylabel near ticks,ymin=0,xtick=0.1,0.2,0.3,0.4, scaled y ticks=base 10:4,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.000000e-01,ymin=0, scaled y ticks=base 10:3,
ylabel near ticks,xtick=0.1,0.2,0.3,0.4,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=20$

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.000000e-01,ymin=0, xtick=0.1,0.2,0.3,0.4,scaled y ticks=base 10:4,
ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=50$

MSE deviation from QDA true error with respect to Bayes error. Source sample size was set to $n_{s}=200$ in all figures.

Figure <ref> shows the MSE deviation with respect to the Bayes error for dimensions 2, 3, and 5. Results in the first column are obtained using target datasets of size 20 and those in the second column are obtained using target datasets of size 50. We can see that the Bayesian MMSE estimator performs best when using source data of high relatedness to the target domain as expected. For Bayes error in the range $\left [ 0.25, 0.35\right ]$, the MSE deviation from the true error is very high, which makes this range of Bayes error as the most challenging setting for error estimation. For a Bayes error of 0.2, the MSE deviation is average across all the experiments which confirms the validity of our previous assumption in selecting this value to investigate classification problems of moderate difficulty.
We note that the TL-based BEE shifts the performance in favor of low and high Bayes error levels. Indeed, the TL-based BEE performs well in this case because the estimated target parameters are sufficiently accurate, even with a small target sample.

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,scaled y ticks=base 10:3,
legend columns=3,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named4,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=10$, $n_s n_s$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:2,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=10$, $n_s n_t$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,scaled y ticks=base 10:3,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=30$, $n_s n_s$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:2,
xlabel near ticks,

legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=30$, $n_s n_t$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:3,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=50$, $n_s n_s$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:2,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$n_t=50$, $n_s n_t$

MSE deviation from true error

with respect to source sample size. The source class means are flipped with respect to target classes ($\mathbf{m}_{s}^{y} = \mathbf{m}_{t}^{1-y}$, for $y\in \left\{0,1\right\}$). In the first row, the source datasets are correctly considered as source samples. In the second row, the source datasets are intentionally considered as target samples. The Bayes error is fixed at 0.2 and $d=5$.

In addition to investigating the effect of different relatedness levels between source and target domains, we have examined in Fig. <ref> the performance of the TL-based BEE for the case  when the source class means are swapped between the two classes, such that they show opposite trends compared to the class means in the target domain. For this purpose, we reproduced the experiments in Fig. <ref> after flipping the class means of source datasets with respect to the target classes (i.e., $\mathbf{m}_{s}^{y} = \mathbf{m}_{t}^{1-y}$, for $y\in \left\{0,1\right\}$). In the first column of Fig. <ref>, we use the generated source datasets as observed samples from the source domain. Interestingly, the obtained results match those observed in Fig. <ref>. This postulates that the knowledge transfer across source and target domains in the context of the studied Bayesian TL framework does not depend on the arrangement of the class means in the source and target domains but only rests on the level of relatedness between the two domains. For verification, we have intentionally considered the same source datasets in the previous experiment as target datasets for estimating the TL-based BEE and we plotted the obtained results in the second column of Fig. <ref>.
Clearly, the TL-based BEE veers away from the true error as we consider additional source data points. This deviation is worse with poorly related source data ($\left|\alpha\right|=0.1$). These results confirm previous findings in [23] that the joint prior model in the utilized Bayesian TL framework acts like a bridge that distills the useful knowledge from the source domain and effectively transfers it to the target domain.

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
xmin=5.000000e-02,xmax=5.000000e-01,ymin=0,ymax=3e-03,scaled y ticks=base 10:4,

legend columns=3,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named6,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=5.000000e-01,ymin=0,scaled y ticks=base 10:4,

ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error,xtick=0.1,0.2,0.3,0.4, ylabel= MSE,ylabel near ticks, scaled y ticks=base 10:4,

xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.500000e-01,ylabel near ticks,ymin=0,xtick=0.1,0.2,0.3,0.4, scaled y ticks=base 10:4,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.000000e-01,scaled y ticks=base 10:4,ymax=2e-03,

ylabel near ticks,xtick=0.1,0.2,0.3,0.4,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=20$

[ width=1,font=,cycle list name=colors_qda,
xlabel= Bayes error, ylabel= MSE,
xmin=5.000000e-02,xmax=4.000000e-01, xtick=0.1,0.2,0.3,0.4,scaled y ticks=base 10:4,

ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=50$

MSE deviation from LDA true error for Gaussian distributions with respect to Bayes error. Source sample size is set to $n_{s}=200$ in all figures.

Results from the second set of experiments that use an LDA classifier were similar as the ones obtained using the QDA classifier except for some differences in the performance of the TL-based BEE with respect to the Bayes error that we report in Fig. <ref> (see Appendix E for additional results).
The TL-based BEE performance has similar trends with respect to small and moderate Bayes errors when compared to the presented results obtained using the QDA classifier.
A notable difference here is observed for large values of Bayes error where the TL-based BEE shows decreased performance in terms of MSE deviation from the true error, which is due to the fact that the employed LDA classifier is sub-optimal compared to the Bayes classifier. This is expected as linear decision boundaries tend to be more sensitive to deviations from true model parameters for highly overlapping class-conditional distributions.

spy using outlines=
connect spies,
[width=1,font=,cycle list name=colors_obtl,
xlabel= $n_t$, ylabel= MSE,ylabel near ticks,scaled y ticks=base 10:3,ymin=0,
xlabel near ticks,

legend columns=-1,
legend entries=resub, cv, loo, boot, bee,
legend to name=named8,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

[anchor=west] (source) at (axis cs:35,7e-03)$n_s=50$;
(destination) at (axis cs:22,0.0012);

[anchor=west] (source) at (axis cs:33.5,3.5e-03)$n_s=200$;
(destination) at (axis cs:30,0.0005);

(point) at (axis cs:12.2,4.9e-03);
(spy point) at (axis cs:22,8e-03);
on (point) in node [fill=white]  at (spy point);

spy using outlines=
connect spies,
[width=1,font=,cycle list name=colors_obtl,
xlabel= $n_t$, ylabel= MSE,
xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,scaled y ticks=base 10:3,ymin=0,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

[anchor=west] (source) at (axis cs:35,7e-03)$n_s=50$;
(destination) at (axis cs:22,0.00125);

[anchor=west] (source) at (axis cs:33.5,3.5e-03)$n_s=200$;
(destination) at (axis cs:30,0.0005);

(point) at (axis cs:13,4.9e-03);
(spy point) at (axis cs:22,8e-03);
on (point) in node [fill=white]  at (spy point);

spy using outlines=
connect spies,
[width=1,font=,cycle list name=colors_obtl,
xlabel= $n_t$, ylabel= MSE,
xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,scaled y ticks=base 10:3,ymin=0,ymax=14e-03,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

[anchor=west] (source) at (axis cs:35,7e-03)$n_s=50$;
(destination) at (axis cs:22,0.00125);

[anchor=west] (source) at (axis cs:33.5,3.5e-03)$n_s=200$;
(destination) at (axis cs:30,0.0005);

(point) at (axis cs:13.5,4.9e-03);
(spy point) at (axis cs:22.2,8e-03);
on (point) in node [fill=white]  at (spy point);

MSE deviation from true error with respect to target data size. The proposed TL-based BEE is compared with other widely used estimators. In all figures, the Bayes error is fixed at $0.2$, and $\left|\alpha\right|=0.9$.

In our final set of experiments using synthetic datasets, we compare the performance of the proposed TL-based BEE to standard error estimators for different dimensions and various source datasets of relatedness level $\left|\alpha\right|=0.9$ to the target domain for an OBTL classifier. We show in Fig. <ref> the MSE deviation with respect to different target dataset size. As clearly shown, our proposed TL-based BEE significantly outperforms all other standard error estimators by a substantial margin.
In agreement with previous findings in the literature, the standard error estimators perform comparably for low dimensions (i.e., $d=2$), where the bootstrap may show a slight advantage. As the dimensionality increases (i.e., $d=5$), the performance shift of the studied estimators becomes more apparent. For example, the resubstitution estimator performs poorly in the small sample regime while the bootstrap estimator outperforms LOO and CV. Furthermore, we could notice that increasing the size of the source dataset does not lead to any apparent performance improvement for the standard estimators. This is because these estimators do not directly depend on the source data for error estimation (as they are incapable of taking advantage of data from different yet relevant domains).

However, providing additional source data to the TL-based BEE considerably reduces the MSE deviation from the true error for all dimensions as shown in Fig. <ref>.

§.§ Performance on real-world RNA-seq datasets

To analyze the performance of the TL-based BEE on real-world data, we have trained a QDA classifier on a small target dataset that consists of five sample points per class extracted from syn2759792 in Table <ref>. Using different source datasets collected from syn4590909, we show in Fig. <ref> the MSE deviation of the TL-based BEE from the true error with respect to $\left|\alpha\right|$.

[width=1,font=,cycle list name=colors_scz,
xlabel= $\left|{\alpha}\right|$, ylabel= MSE,ylabel near ticks,xmin=0,xmax=1,xtick=0,0.2,0.4,0.6,0.8,1,
xlabel near ticks,
scaled y ticks=base 10:4,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=5$

[ width=1,font=,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
scaled y ticks=base 10:4,
grid style=dotted,

[blue!80!black,solid, thick]plot coordinates

$d=3$, $n_t=5$, $\left|{\alpha}\right|=0.99$

MSE deviation from QDA true error for normally distributed brain gene expression data with respect to $\left|\alpha\right|$ and $n_{s}$. a. Gene features from the FC brain region demonstrate high relatedness with those from DPLFC area ($\left|\alpha\right|=0.99$). b. Utilizing the data from source domain significantly reduces the MSE of the TL-based BEE in the target domain.

For all combinations and different sizes of source datasets, the FC brain region showed high relatedness to the DLPFC brain area where the optimal MSE deviation from the true error was obtained for $\left|\alpha\right|=0.99$.
Interestingly, findings in [39] also confirm that syn4590909 and syn2759792 are highly related, as independent gene expression assays for both brain regions have consistently replicated the gradient of transcriptomic severity observed for three different types of psychiatric disorders, including bipolar disorder and Schizophrenia [39]. We show in Fig. <ref> the increasing gain in accuracy of the TL-based BEE in estimating the classification error after using additional labeled observations from the source domain. These results again confirm the efficacy and advantages of our novel TL-based error estimation scheme, compared to other standard error estimation methods, when additional data are available from different source domains that are nevertheless relevant to the target domain.
From a practical perspective, our proposed TL-based BEE has the potential to facilitate the analysis of real-world datasets, in the context of small-sample classification. Challenges of designing and evaluating classifiers (e.g., for clinical diagnosis or prognosis) in a small-sample setting are prevalent in scientific studies in life sciences and physical sciences due to the formidable cost, time, and effort required for data acquisition. This is certainly the case for the example that we considered in this section, where invasive brain biopsies would be needed to get the data.

§ CONCLUSIONS

In this study, we have introduced a novel  Bayesian MMSE estimator that draws from concepts and theories in transfer learning to enable accurate estimation of classification error in the (target) domain of interest by utilizing samples from other closely related (source) domains. We have developed an efficient and robust importance sampling set-up that can be used for accurate error estimation in small-sample scenarios that often arise in many real-world scientific problems. Extensive performance analysis based on both synthetic and real biological data demonstrates the outstanding performance of the proposed TL-based BEE clearly outperforming conventional estimators.

In our proposed framework, Laplace approximations were used to alleviate the complexity associated with the exact evaluation of generalized hypergeometric functions that appear in the posterior distribution of the target parameters.
Beyond the Gaussian model assumed in the validation experiments, we also provided a general mathematical definition for the TL-based BEE that can directly be extended to applications with non-Gaussian distributions where the model parameters can be inferred through MCMC methods.

In this study, target and source domains were related through the joint prior of the model parameters that transfers useful knowledge across domains. A key property of the proposed TL-based BEE is its elegant ability to handle the uncertainty about the model parameters by integrating this prior with data, deducing robust estimates by accounting for all possible parameter values.

Paramount practical challenges for the TL-based BEE include the identification of suitable source domains that share similar families of distributions as the target domain of interest. This is crucial as the relatedness across domains is mathematically modeled assuming the similarity of the feature-label distributions across domains. Furthermore, learning the joint prior for the distributions and modeling the relatedness between different domains may also present an engineering challenge. While techniques for knowledge-driven prior construction have been developed [41, 42], such techniques yet have to be developed for joint prior construction for relevant domains, which is an important future research direction.

§ ACKNOWLEDGMENT

This work was supported in part by the Department of Energy (DOE) under Award DE-SC0019303.

Portions of this research were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.

§ BINARY CLASSIFICATION

A standard binary classification problem involves a sample space $\mathcal{X}$ associated with the real-valued $d$-dimensional feature space $\mathbb{R}^{d}$ where the class labels are modeled by a random variable $y$ in the output space $\mathcal{Y} = \left \{ 0,1 \right \}$.
Classification aims at designing an optimal function $\psi$ that assigns each feature vector of the form $\mathbf{x}=\left ( x_{1}, x_{2},\cdots,x_{d}  \right ) \in \mathcal{X}$ to its corresponding true class-label $y \in \mathcal{Y}$. The joint relationship between $\mathbf{x}$ and $y$ is fully characterised by the feature-label distribution $f\left ( \mathbf{x},y \right )$.

Using the standard Bayesian terminology, the class-conditional density is denoted by $f\left ( \mathbf{x}~|~y \right )$ and the prior probability by $P\left ( y\right )$. In binary classification, the class prior probabilities are fully characterized by one parameter $c=P\left ( y=0\right )=1-P\left ( y=1\right )$.
The classification error $\varepsilon$ of a classifier $\psi(\cdot)$ is determined either using the feature-label distribution $f\left ( \mathbf{x},y \right )$ or the class-conditional distribution $f\left ( \mathbf{x}~|~y \right )$.

Using $f\left ( \mathbf{x},y \right )$, $\varepsilon$ represents the probability of misclassification given by

\begin{equation}
P\left ( \psi\left ( \mathbf{x} \right )\neq y \right )=E\left [ \left | y- \psi\left ( \mathbf{x} \right )\right | \right ].
\end{equation}

Using the class-conditional densities $f\left ( \mathbf{x}~|~y \right )$, $\varepsilon$ is decomposed into two error components contributed by each class and denoted by $\varepsilon^{y}$.

Thus, $\varepsilon$ is expressed as

\begin{equation}
\varepsilon = c~\varepsilon^{0}+\left ( 1-c\right )~\varepsilon^{1},
\end{equation}

\begin{align}
\begin{split}
\varepsilon^{y}&=P\left( \psi\left(\mathbf{x}\right)=1-y~|~y\right)\\
&=\int_{\psi\left(\mathbf{x}\right)=1-y}f\left( \mathbf{x}~|~y\right)~d\mathbf{x}.
\end{split}
\end{align}

When the feature-label distribution is unknown, a sample $\mathcal{S}_{n}=\left \{ \left ( \mathbf{x}_{1},y_{1}\right ), \left ( \mathbf{x}_{2},y_{2}\right ), \cdots, \left ( \mathbf{x}_{n},y_{n}\right )\right\}$ of size $n$ is drawn from $f\left ( \mathbf{x},y \right )$ through measurements collected from the data generating process.

Instead of finding $\left (\psi, \varepsilon \right )$, the classification problem turns out to be a search for another pair $\left (\psi_{n}, \varepsilon_{n} \right )$ that is a function of the observed sample $\mathcal{S}_{n}$. As $\varepsilon_{n}$ is also unknown for unknown feature-label distributions, an estimate $\hat{\varepsilon}$ that accurately approximates $\varepsilon_{n}$ is considered. Designing $\psi_{n}$ is generally possible through a classification rule ($\Psi\left( \mathcal{S}_{n}\right )=\psi_{n}$) defined by

\begin{equation}
\Psi :\left [ \mathbb{R}^{d} \times \left \{ 0,1 \right \}\right ]^{n}\rightarrow \mathcal{C},
\end{equation}

\begin{equation}
\mathcal{C}=\left \{ \psi\ |\ \psi : \mathbb{R}^{d}\rightarrow \left \{ 0,1 \right \} \right \}
\end{equation}

is the class containing all possible classifiers.

Similarly, deriving $\hat{\varepsilon}$ is possible through an estimation rule ($\Upsilon\left( \mathcal{S}_{n}\right )=\hat{\varepsilon}$) defined by

\begin{equation}
\Upsilon :\left [ \mathbb{R}^{d} \times \left \{ 0,1 \right \}\right ]^{n}\rightarrow \left [ 0,1\right ].
\end{equation}

From an error estimation perspective, we lock our focus on evaluating the performance of $\hat{\varepsilon}$ regardless of the design of $\psi_{n}$.

§ BAYESIAN MMSE ESTIMATOR

The BEE has been introduced for the first time in [3, 17] as a parametric estimator that assumes uncertainty in the parameters of a family of feature-label distributions $f\left ( \mathbf{x},y \right )$.
For a specific sample realization $S_{n} \in \mathcal{S}_{n}$, and a parameter vector $\mathbf{\theta}$, for which the notation of the feature-label distribution reduces to $f_{\mathbf{\theta}}\left ( \mathbf{x},y \right )$, the BEE is given by

\begin{equation}\label{BEE_long}
\hat{\varepsilon}\left (S_{n} \right )=E_{\mathbf{\theta}}\left[\varepsilon_{n}\left( \mathbf{\theta},S_{n}~|~S_{n}\right) \right].
\end{equation}

After omitting the sample $S_{n}$ from the notation of the posterior density $\pi^{*}\left( \theta~|~S_{n} \right)$ for simplicity, the BEE is shortly expressed as

\begin{equation}
\hat{\varepsilon}=E_{\pi^{*}}\left[\varepsilon_{n} \right].
\end{equation}

The BEE has proven to be optimal when averaged over a given family of feature-label distributions and unbiased when averaged over a given family and all samples [3].
In binary classification the parameter vector $\theta$ can be expressed using three components ($\theta=\left [\theta_{0}, \theta_{1}, c \right ]$), where $\theta_{0}$ and $\theta_{1}$ model the class conditional distributions and $c$ represents the prior probability of class 0.

Assuming that $c$, $\theta_{0}$, and $\theta_{1}$ are independent prior to observing the data, the BEE can be further decomposed as

\begin{equation}
\hat{\varepsilon}=E_{\pi^{*}}\left[c \right]~E_{\pi^{*}}\left[\varepsilon_{n}^{0} \right]+\left(1-E_{\pi^{*}}\left[c \right]\right)~E_{\pi^{*}}\left[\varepsilon_{n}^{1} \right],
\end{equation}

\begin{equation}
\end{equation}

with $\Theta_{y}$ being the parameter space that contains all possible values for $\theta_{y}$.

§ LAPLACE APPROXIMATIONS OF CONFLUENT AND GAUSS HYPERGEOMETRIC
FUNCTIONS OF MATRIX ARGUMENT

Definition 1 [31]: The generalized hypergeometric function of one matrix argument is defined by

\begin{equation}\label{HG_series}
\pFq{p}{q}{a_{1}, \cdots, a_{p}}{b_{1}, \cdots, b_{q}}{\mathbf{X}}=\sum_{k=0}^{\infty}\sum_{\kappa \vdash k}\frac{\left(a_{1}\right)_{\kappa}\cdots \left(a_{p}\right)_{\kappa}}{\left(b_{1}\right)_{\kappa}\cdots \left(b_{q}\right)_{\kappa}}\frac{C_{\kappa}\left(\mathbf{X}\right)}{k!},
\end{equation}

where $\left(a_{i}\right)_{i=1\cdots p}$ and $\left(b_{j}\right)_{j=1\cdots q}$ are arbitrary complex numbers, $C_{\kappa}\left(\mathbf{X}\right)$ is the zonal polynomial of the complex symmetric matrix $\mathbf{X}$ of order $d$ that corresponds to the ordered partition $\kappa$ of $k$ defined as $\kappa=\left(k_{1}, \cdots, k_{d}\right),~k_{1}\geq \cdots \geq k_{d} \geq 0,~\sum_{i=1}^{d}k_{i}=k$, and $\sum_{\kappa \vdash k}$ denotes the summation over all ordered partitions $\kappa$ of $k$. The generalized hypergeometric coefficient $\left(a\right)_{\kappa}$ is defined by

\begin{equation}
\left(a\right)_{\kappa}=\prod_{i=1}^{d}\left(a-\frac{i-1}{2}\right)_{k_{i}},
\end{equation}

where $\left(a\right)_{r}=a\left(a+1\right)\cdots\left(a+r-1\right),~r=1,~2\cdots,$ with $\left(a\right)_{0}=1$.

Conditions for convergence of the series in (<ref>) are covered in [32]. Some special cases of matrix-variate generalized hypergeometric functions include

\begin{align}
\begin{split}
\pFq{0}{1}{.}{b}{\mathbf{X}}&=\sum_{k=0}^{\infty}\sum_{\kappa \vdash k}\frac{C_{\kappa}\left(\mathbf{X}\right)}{\left(b\right)_{\kappa}~k!}; \\
\pFq{1}{1}{a}{b}{\mathbf{X}}&=\sum_{k=0}^{\infty}\sum_{\kappa \vdash k}\frac{\left(a\right)_{\kappa}}{\left(b\right)_{\kappa}}\frac{C_{\kappa}\left(\mathbf{X}\right)}{k!}; \\
\pFq{2}{1}{a,b}{c}{\mathbf{X}}&=\sum_{k=0}^{\infty}\sum_{\kappa \vdash k}\frac{\left(a\right)_{\kappa}~\left(b\right)_{\kappa}}{\left(c\right)_{\kappa}}\frac{C_{\kappa}\left(\mathbf{X}\right)}{k!},~\left \| \mathbf{X} \right \|<1,
\end{split}
\end{align}

where $\left \| \mathbf{X} \right \|$ denotes the maximum of the absolute values of the eigenvalues of $\mathbf{X}$. $\pFq{1}{1}{a}{b}{\mathbf{X}}$ and $\pFq{2}{1}{a,b}{c}{\mathbf{X}}$ are respectively called confluent and Gauss matrix-variate hypergeometric functions.

The confluent hypergeometric function has the following integral representation:

\begin{equation}
\pFq{1}{1}{a}{b}{\mathbf{Z}}=B_{d}^{-1}\left(a,b-a\right)\int_{0_{d}<\mathbf{Y}<\mathbf{I}_{d}}\mbox{etr}\left(\mathbf{Z}\mathbf{Y}\right)\left|\mathbf{Y}\right|^{a-\frac{d+1}{2}}\left|\mathbf{I}_{d}-\mathbf{Y}\right|^{b-a-\frac{d+1}{2}}d\mathbf{Y},
\end{equation}

which is valid under the following conditions:

* $\mathbf{Z}\in \mathbb{C}^{d\times d}$ is symmetric;

* $\Re\left(a\right)>\frac{d-1}{2}$;

* $\Re\left(b-a\right)>\frac{d-1}{2}$,

where $B_{d}\left(\alpha,\beta\right)$ is the multivariate beta function defined by:

\begin{equation}
\end{equation}

with $\Gamma\left(a\right)$ denoting the multivariate gamma function defined by:

\begin{equation}
\Gamma\left(a\right)=\int_{0_{d}<\mathbf{Y}<\mathbf{I}_{d}}\mbox{etr}\left(\mathbf{-Y}\right)\left|\mathbf{Y}\right|^{a-\frac{d+1}{2}}d\mathbf{Y}.
\end{equation}

Laplace approximations provide accurate evaluations of integrals of the form:

\begin{equation}
I=\int_{y\in\mathbf{D}}h\left(y\right)\mbox{e}^{-\lambda g\left(y\right)}dy,
\end{equation}

where $\mathbf{D}\subseteq\mathbb{R}^{d}$ is an open set and $\lambda$ is a real-valued parameter, when $g\left(y\right)$ has a unique minimum over the closure of $\mathbf{D}$, and this minimum occurs at stationary point $\hat{y}\in \mathbf{D}$ of $g\left(y\right)$.

This approximation is given by:

\begin{equation}\label{Laplace}
\tilde{I}=\left(2\pi\right)^{\frac{d}{2}}\lambda^{\frac{-d}{2}}\left|{g}''\left(\hat{y}\right)\right|^{\frac{-1}{2}}h\left(\hat{y}\right)\mbox{e}^{-\lambda g\left(\hat{y}\right)},
\end{equation}

where ${g}''\left(y\right)=\frac{\partial^{2}g\left(y\right)}{\partial y\partial y^{T}}$ is the Hessian of $g\left(y\right)$.
In the case of matrix-variate hypergeometric functions, an important invariance property that facilitates the selection of $\left(g;h\right)$ is that $\pFq{1}{1}{a}{b}{\mathbf{Z}}$ and $\pFq{2}{1}{a,b}{c}{\mathbf{Z}}$ depend only on the eigenvalues of $\mathbf{Z}$ when $\mathbf{Z}$ is symmetric. As such, $\mathbf{Z}$ can be assumed diagonal. Furthermore, we focus on the case where $\mathbf{Z}$ is real and we write instead $\mathbf{X}$.

A favorable representation of $\left(g;h\right)$ as discussed in [43] is given by:

\begin{equation}
\end{equation}

\begin{equation}
\end{equation}

Using (<ref>), the Laplace approximation to $\pFq{1}{1}{a}{b}{\mathbf{X}}$ is given by:

\begin{equation}
\pFqtilde{1}{1}{a}{b}{\mathbf{X}}=2^{\frac{d}{2}}\pi^{\frac{d\left(d+1\right)}{4}}B_{d}^{-1}\left(a,b-a\right)J_{1,1}^{-\frac{1}{2}}\prod_{i=1}^{d}\left\{\hat{y}_{i}^{a}\left(1-\hat{y}_{i}\right)^{b-a}\mbox{e}^{x_{i}\hat{y}_{i}}\right\},
\end{equation}

where $\hat{y}_{i}$ is given by:

\begin{equation}
\hat{y}_{i}=\frac{2a}{b-x_{i}+\sqrt{\left(x_{i}-b\right)^{2}+4ax_{i}}},
\end{equation}

\begin{equation}
\end{equation}

A more accurate approximation emerges from the calibrated expression given by:

\begin{equation}
\pFqhat{1}{1}{a}{b}{\mathbf{X}}=\frac{\pFqtilde{1}{1}{a}{b}{\mathbf{X}}}{\pFqtilde{1}{1}{a}{b}{\mathbf{0}}}
\end{equation}

since $\pFq{1}{1}{a}{b}{\mathbf{0}}=1$. Consequently, the calibrated Laplace approximation to $\pFq{1}{1}{a}{b}{\mathbf{X}}$ is given by:

\begin{equation}
\pFqhat{1}{1}{a}{b}{\mathbf{X}}=b^{bd-\frac{d\left(d+1\right)}{4}}R_{1,1}^{-\frac{1}{2}}\prod_{i=1}^{d}\left\{\left(\frac{\hat{y}_{i}}{a}\right)^{a}\left(\frac{1-\hat{y}_{i}}{b-a}\right)^{b-a}\mbox{e}^{x_{i}\hat{y}_{i}}\right\},
\end{equation}

\begin{equation}
R_{1,1}=\prod_{i=1}^{d}\prod_{j=1}^{d}\left\{\frac{\hat{y}_{i}\hat{y}_{j}}{a}+\frac{\left ( 1-\hat{y}_{i} \right )\left ( 1-\hat{y}_{j} \right )}{b-a}\right\}.
\end{equation}

Similarly, the Gauss hypergeometric function has the following representation:

\begin{equation}
\pFq{2}{1}{a,b}{c}{\mathbf{Z}}=B_{d}^{-1}\left(a,c-a\right)\int_{0_{d}<\mathbf{Y}<\mathbf{I}_{d}}\left|\mathbf{Y}\right|^{a-\frac{d+1}{2}}\left|\mathbf{I}_{d}-\mathbf{Y}\right|^{c-a-\frac{d+1}{2}}\left|\mathbf{I}_{d}-\mathbf{Z}\mathbf{Y}\right|^{-b}d\mathbf{Y},
\end{equation}

which is valid under the following conditions:

* $\mathbf{Z}\in \mathbb{C}^{d\times d}$ is symmetric;

* $\Re\left(\mathbf{Z}\right)<\mathbf{I}_{d}$;

* $\Re\left(a\right)>\frac{d-1}{2}$;

* $\Re\left(c-a\right)>\frac{d-1}{2}$.

By following similar steps as in [23], the calibrated Laplace approximation to $\pFq{2}{1}{a,b}{c}{\mathbf{X}}$ is given by:

\begin{align}
\pFqhat{2}{1}{a,b}{c}{\mathbf{X}}=c^{cd-\frac{d\left(d+1\right)}{4}}R_{2,1}^{-\frac{1}{2}}\prod_{i=1}^{d}\left\{\left(\frac{\hat{y}_{i}}{a}\right)^{a}\left(\frac{1-\hat{y}_{i}}{c-a}\right)^{c-a}\left(1-x_{i}\hat{y}_{i}\right)^{-b}\right\},
\end{align}

\begin{equation}
\end{equation}

Detailed discussions in [43] show that the relative errors of the approximations to $\pFq{1}{1}{a}{b}{\mathbf{X}}$ and $\pFq{2}{1}{a,b}{c}{\mathbf{X}}$ are uniformly bounded. These bounds are as follows:

\begin{equation}
\operatorname*{sup}_{b\geq b_{0},~a\in \mathbb{R},~\mathbf{X}\in \mathbb{R}^{d\times d}}\left|\log\left(\pFqhat{1}{1}{a}{b}{\mathbf{X}}\right)-\log\left(\pFq{1}{1}{a}{b}{\mathbf{X}}\right)\right|<\infty
\end{equation}

\begin{equation}
\operatorname*{sup}_{c\geq c_{0},~a,b\in \mathbb{R},~\mathbf{0}_{d}\leq\mathbf{X}< \left(1-\varepsilon\right)\mathbf{I}_{d}}\left|\log\left(\pFqhat{2}{1}{a,b}{c}{\mathbf{X}}\right)-\log\left(\pFq{2}{1}{a,b}{c}{\mathbf{X}}\right)\right|<\infty
\end{equation}

given that $d$ is fixed, $b_{0}>\frac{d-1}{2}$, $c_{0}>\frac{d-1}{2}$, and $\forall \varepsilon \in \left(0,1\right)$.

Similar to the examples provided in [43] and the illustrations presented in [23], we provide simulations in Fig. <ref> to show the good numerical accuracy of the discussed Laplace approximations as compared to the exact evaluations.

xlabel= $\tau$, ylabel= Function Value,ylabel near ticks,
xlabel near ticks,xmin=0,xmax=0.8,ymin=1,ytick=1,10,xtick=0,0.2,0.4,0.6,0.8,
legend style=nodes=right,legend style=nodes=scale=1, transform shape, legend pos= north west,
grid style=dotted,

[smooth,blue,thick,mark=o,mark options=solid,every mark/.append style=solid, fill=white,densely dashed] plot coordinates

[smooth, thick,mark=star,red,only marks]plot coordinates
Exact confluent HG function,Laplace approximation;

$d=5$, $a=3$, $b=4$.

[ font=,ymode=log,
xlabel= $\tau$, xmin=0,xmax=0.8,ymin=1,ymax=1e07,ylabel= Function Value,ylabel near ticks,xtick=0,0.2,0.4,0.6,0.8,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,legend pos= north west,
grid style=dotted,

[smooth,blue,thick,mark=o,mark options=solid,every mark/.append style=solid, fill=white,densely dashed] plot coordinates

[smooth, thick,mark=star,red,only marks]plot coordinates
Exact Gauss HG function,Laplace approximation;

$d=5$, $a=3$, $b=4$, $c=6$.

[ font=,ymode=log,
scaled ticks=false,ytick=1.05,1.1,1.15,1.2,1.25,1.3,yticklabels=1.05,1.1,1.15,1.2,1.25,1.3,
xlabel= $b$, xmin=10,xmax=50,ymin=1.05,ylabel= Function Value,ylabel near ticks,
scaled y ticks = false,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

[smooth,blue,thick,mark=o,mark options=solid,every mark/.append style=solid, fill=white,densely dashed] plot coordinates

[smooth, thick,mark=star,red,only marks]plot coordinates
Exact confluent HG function,Laplace approximation;

$d=10$, $a=30$, $\tau=0.01$.

[ font=,ymode=log,
xlabel= $c$,xmin=10,xmax=50,ymin=10, ylabel= Function Value,ylabel near ticks,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

[smooth,blue,thick,mark=o,mark options=solid,every mark/.append style=solid, fill=white,densely dashed] plot coordinates

[smooth, thick,mark=star,red,only marks]plot coordinates
Exact Gauss HG function,Laplace approximation;

$d=10$, $a=30$, $b=50$, $\tau=0.01$.

Laplace approximations and exact evaluations of confluent and Gauss hypergeometric functions of matrix argument

As the hypergeometric functions of matrix argument depend only on eigenvalues of their argument, we set $\mathbf{X}=\tau\mathbf{I}_{d}$ and we evaluate in Figs <ref> and <ref> the accuracy of the approximations to $\pFq{1}{1}{a}{b}{\tau\mathbf{I}_{d}}$ and $\pFq{2}{1}{a,b}{c}{\tau\mathbf{I}_{d}}$, respectively, for $d=5,~a=3,~b=4,~c=6$ and w.r.t $\tau$ where $0<\tau<1$ guarantees the convergence as mentioned in [43].

Figs. <ref> and <ref> show the exact and approximate values of confluent and Gauss hypergeometric functions for variable $b$ and $c$ values, respectively, when $d=10, a=30, \left(b=50,~\mbox{for}~\pFq{2}{1}{a,b}{c}{\tau\mathbf{I}_{d}}\right)$, and $\tau=0.01$. As the integral representations of $\pFq{1}{1}{a}{b}{\tau\mathbf{I}_{d}}$ and $\pFq{2}{1}{a,b}{c}{\tau\mathbf{I}_{d}}$ are only valid when $\Re\left(b-a\right)>\frac{d-1}{2}$ and $\Re\left(c-a\right)>\frac{d-1}{2}$, respectively, the simulation results show that the Laplace approximations still give good accuracy even for values of $b$ and $c$ where the integral is not valid. This behaviour has been discussed thoroughly in [43] and was explained as an advantage of the calibrated approximations $\pFqhat{1}{1}{a}{b}{\tau\mathbf{I}_{d}}$ and $\pFqhat{2}{1}{a,b}{c}{\tau\mathbf{I}_{d}}$ where the singularities of the invalid integrals are removed.

§ MONTE CARLO SAMPLING

We consider the problem of approximating:

\begin{equation}
\end{equation}

for an integrand $f\left(x\right)$ and a probability density $p\left(x\right)$ that are defined on $\mathcal{D}\subseteq \mathbb{R}^{d}$ where $E_{p}\left[\cdot\right]$ denotes the expectation for $\mathbf{x}\sim p$.

§.§ Importance sampling

If q is a positive probability density function on $\mathbb{R}^{d}$, we can write:

\begin{equation}
I=\int_{\mathcal{D}}\frac{f\left(x\right)p\left(x\right)}{q\left(x\right)}q\left(x\right) dx=E_{q}\left[\frac{f\left(\mathbf{x}\right)p\left(\mathbf{x}\right)}{q\left(\mathbf{x}\right)}\right].
\end{equation}

By sampling $\left(\mathbf{x}_{i}\right)_{i=1\cdots n}\sim q$, the importance sampling estimate of $I$ is given by:

\begin{equation}\label{eq_Iq}
\hat{I}_{q}= \frac{1}{n}\sum_{i=1}^{n}\frac{f\left(\mathbf{x}_{i}\right)p\left(\mathbf{x}_{i}\right)}{q\left(\mathbf{x}_{i}\right)}.
\end{equation}

Theorem 1 [44]: Let $\hat{I}_{q}$ be given by (<ref>) where $q\left(x\right)>0$ whenever $f\left(x\right)p\left(x\right)\neq0$. Then $E_{q}\left[\hat{I}_{q}\right] = I$, and $\mbox{Var}_{q}\left(\hat{I}_{q}\right)=\frac{\sigma_{q}^{2}}{n}$ where:

\begin{equation}
\sigma_{q}^{2} = \int_{\mathcal{Q}}\frac{\left(f\left(x\right)p\left(x\right)\right)^{2}}{q\left(x\right)}dx-I^{2}=\int_{\mathcal{Q}}\frac{\left(f\left(x\right)p\left(x\right)-Iq\left(x\right)\right)^{2}}{q\left(x\right)}dx,
\end{equation}

where $\mathcal{Q}=\left\{x~|~q\left(x\right)>0\right\}$.

§.§ Self-normalized importance sampling

When $p$ or $q$ has an unknown normalization constant, we resort to estimate the ratio $\frac{p\left(x\right)}{q\left(x\right)}$. Suppose we can compute an unnormalized ratio $\frac{p_{u}\left(x\right)}{q_{u}\left(x\right)}$ where $p_{u}\left(x\right) = a~p\left(x\right)$ and $q_{u}\left(x\right) = b~q\left(x\right)$ for $a,b>0$.

Then, we compute the ratio $\mathcal{L}_{u}\left(x\right)=\frac{p_{u}\left(x\right)}{q_{u}\left(x\right)}=\frac{a}{b}\frac{p\left(x\right)}{q\left(x\right)}$ and we consider the self-normalized importance sampling estimate given by:

\begin{equation}\label{eq_I_self_norm}
\tilde{I}_q=\frac{\sum_{i=1}^{n}f\left(\mathbf{x}_{i}\right)\mathcal{L}_{u}\left(\mathbf{x}_{i}\right)}{\sum_{i=1}^{n}\mathcal{L}_{u}\left(\mathbf{x}_{i}\right)},
\end{equation}

where $\mathbf{x}_{i}\sim q$ are independent. The factor $\frac{a}{b}$ cancels from the numerator and denominator leading to:

\begin{equation}
\tilde{I}_q=\frac{\sum_{i=1}^{n}f\left(\mathbf{x}_{i}\right)\mathcal{L}\left(\mathbf{x}_{i}\right)}{\sum_{i=1}^{n}\mathcal{L}\left(\mathbf{x}_{i}\right)},
\end{equation}

where $\mathcal{L}\left(x\right)=\frac{p\left(x\right)}{q\left(x\right)}$.

Theorem 2 [44]: Let p be a probability density function on $\mathbb{R}^{d}$ and let $f\left(x\right)$ be a function such that $I=\int f\left(x\right)p\left(x\right)dx$ exists. Suppose that $q\left(x\right)$ is a probability density function on $\mathbb{R}^{d}$ with $q\left(x\right)>0$ whenever $p\left(x\right)>0$. Let $\left(\mathbf{x}_{i}\right)_{i=1\cdots n}\sim q$ be independent and let $\tilde{I}_{q}$ be the self-normalized importance estimate in (<ref>). Then

\begin{equation}
P\left(\operatorname{lim}_{n\to \infty}\tilde{I}_{q}=I\right)=1.
\end{equation}

§.§ Control variates in importance sampling

A control variate function $h\left(x\right)$ can be usefully combined with importance sampling if it satisfies $\int h\left(x\right)q\left(x\right)dx=\delta$, where $\delta$ is a known constant.

A more stable estimate to $I$ is given by:

\begin{equation}
\tilde{I}_{q,\beta}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{f\left(\mathbf{x}_{i}\right)p\left(\mathbf{x}_{i}\right)}{q\left(\mathbf{x}_{i}\right)}- \beta h\left(\mathbf{x}_{i}\right)\right]+\beta \delta,
\end{equation}

where $\left(\mathbf{x}_{i}\right)_{i=1\cdots n}\sim q$ and $\beta$ is a hyperparameter that controls the variance reduction.

Theorem 3 [44]: Let $q$ be a probability density function with $q\left(x\right)>0$ whenever $h\left(x\right)\neq0$ or $f\left(x\right)p\left(x\right)\neq0$. Then $E_{q}\left[\tilde{I}_{q,\beta}\right]=I$ for any $\beta \in \mathbb{R}$.

Let $\varsigma\left(x\right)=\frac{f\left(x\right)p\left(x\right)}{q\left(x\right)}$. The optimal value of $\beta$ is chosen in a way to minimize the variance of the random variable $Z = \varsigma\left(\mathbf{x}\right)-\beta\left(h\left(\mathbf{x}\right)-\delta\right)$.

We have:

\begin{equation}
\mbox{Var}\left(Z\right) = \mbox{Var}\left(\varsigma\right) - 2\beta \mbox{Cov}\left[\varsigma, h\right]+\beta^{2}\mbox{Var}\left(h\right).
\end{equation}

Solving for $\beta$ we get:

\begin{equation}
\beta^{\star}=\frac{\mbox{Cov}\left[\varsigma, h\right]}{\mbox{Var}\left(h\right)}.
\end{equation}

The reduced variance is then given by:

\begin{equation}
\mbox{Var}\left(Z\right)=\mbox{Var}\left(\varsigma\right)-\frac{\mbox{Cov}^{2}\left[\varsigma, h\right]}{\mbox{Var}\left(h\right)} = \mbox{Var}\left(\varsigma\right)\left(1-\rho^{2}\left[\varsigma, h\right]\right),
\end{equation}

where $\rho\left[\varsigma, h\right]$ is a correlation coefficient given by:

\begin{equation}
\rho\left[\varsigma, h\right]=\mbox{Corr}\left[\varsigma, h\right]=\frac{\mbox{Cov}\left[\varsigma, h\right]}{\mbox{Var}\left(\varsigma\right)\mbox{Var}\left(h\right)}.
\end{equation}

In practice $\beta^{\star}$ is unknown. Instead, it can be estimated from the Monte Carlo samples. For $n$ independent samples $\left(\mathbf{x}_{i}\right)_{i=1\cdots n} \sim q$ we can compute:

\begin{equation}
\left\{\begin{matrix}
\widehat{\mbox{Var}}\left(h \right )=\frac{1}{n}\sum_{i=1}^{n}\left(h\left(\mathbf{x}_{i}\right )-\delta \right )^{2}\\
\widehat{I}=\frac{1}{n}\sum_{i=1}^{n}\varsigma\left(\mathbf{x}_{i}\right)\\
\widehat{\mbox{Cov}}\left[\varsigma,h \right ]=\frac{1}{n}\sum_{i=1}^{n}\left(\varsigma\left(\mathbf{x}_{i}\right)- \widehat{I}\right )\left(h\left(\mathbf{x}_{i}\right)- \delta\right )
\end{matrix}\right.
\end{equation}

and thus we determine:

\begin{equation}
\widehat{\beta}^{\star}=\frac{\widehat{\mbox{Cov}}\left[\varsigma, h\right]}{\widehat{\mbox{Var}}\left(h\right)}.
\end{equation}

§.§ Integral evaluation of LDA error expectation as a control variate function

From [17], we know that the posterior density of model parameters $\left(\mu^{y}, \Sigma^{y}\right)$ after observing the sample $S_{n}$ is a Gaussian inverse-Wishart density of covariance $\Sigma^{y}\sim \mathcal{W}^{-1}\left(\mathbf{S}^{y}, \nu^{y}\right)$ and mean $\mu^{y}\sim \mathcal{N}\left(\mathbf{m}^{y}, \frac{\Sigma^{y}}{\kappa^{y}}\right)$.

For an LDA classifier defined by:

\begin{equation}
\psi_{n}\left(x \right )=\left\{\begin{matrix}
0, &\mbox{if }g\left(x \right )\leq0 \\
1, &otherwise
\end{matrix}\right.
\end{equation}

where $g\left(x\right)=\mathbf{a}^{T}x+b$ with $\mathbf{a}$ and $b$ are a constant vector and a constant scalar, respectively, the standard BEE has been derived in closed-form as:

\begin{equation}
E_{\pi^{*}}\left[\varepsilon _{n}^{y}\right]=\frac{1}{2}+\frac{\mbox{sgn}\left(A\right)}{2}\mathcal{I}\left(\frac{A^{2}}{A^{2}+{\mathbf{a}}^{T}\mathbf{S}^{y}\mathbf{a}};\frac{1}{2},\frac{\nu^{y}-d+1}{2}\right),
\end{equation}

where $\mbox{sgn}\left(\cdot\right)$ is the sign function,

\begin{equation}
\end{equation}

and $\mathcal{I}\left(\cdot;\cdot,\cdot\right)$ denotes the regularized incomplete beta function given by

\begin{equation}
\mathcal{I}\left(x;a,b\right)=\frac{\Gamma\left(a+b\right)}{\Gamma\left(a\right)\Gamma\left(b\right)}\int_{0}^{x}t^{a-1}\left(1-t\right)^{b-1}dt,
\end{equation}

with $\Gamma\left(\cdot\right)$ being the regular univariate gamma function.

In our TL Bayesian framework we have: $\Phi^{*}= \pi^{*}\left(\mu_{t}^{y}, \Lambda_{t}^{y}\right)=\pi^{*}\left(\mu_{t}^{y}| \Lambda_{t}^{y}\right)\pi^{*}\left(\Lambda_{t}^{y}\right)$.

Clearly, $\Phi^{*}$ is a Gaussian Wishart posterior density of model parameters $\left(\mu_{t}^{y}, \Lambda_{t}^{y}\right)$ with a precision matrix $\Lambda_{t}^{y} \sim \mathcal{W}\left(\mathbf{M}_{t,n}^{y}, \nu^{y}+n_{t}^{y}\right)$ and a mean $\mu_{t}^{y}\sim \mathcal{N}\left(\mathbf{m}_{t,n}^{y}, \left(\kappa_{t,n}^{y}\Lambda_{t}^{y}\right)^{-1}\right)$.

After making the change of variable $\Sigma_{t}^{y}=\left(\Lambda_{t}^{y}\right)^{-1}$ we get:

$\Sigma_{t}^{y} \sim \mathcal{W}\left(\left[\mathbf{M}_{t,n}^{y}\right]^{-1}, \nu^{y}+n_{t}^{y}\right)$ and a mean $\mu_{t}^{y}\sim \mathcal{N}\left(\mathbf{m}_{t,n}^{y}, \frac{\Sigma_{t}^{y}}{\kappa_{t,n}^{y}}\right)$.
By analogy, the rest of the derivation follows from [17].

§ ADDITIONAL RESULTS FOR LINEAR CLASSIFIERS

In Fig. <ref> we consider an LDA classifier and we investigate the behavior of the TL-based BEE when the target data are fixed while we vary the size of the source data. We show in the left column results for $d=2$, in the middle column results for $d=3$, and in the right column results for $d=5$. The rows correspond to results for different sizes of target datasets. Obviously, The MSE has similar trends across the three dimensions as compared to the QDA results in the main text. The deviation from the true error significantly decreases when highly related source data are employed.

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,

legend columns=-1,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named5,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=20$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,ylabel near ticks,scaled y ticks=base 10:4,

xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=20$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:4,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=20$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_t=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,
xlabel near ticks,

legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_t=50$

[width=1, font=,cycle list name=colors_qda,
xlabel= $n_s$, ylabel= MSE,
xmin=10,xmax=500,xtick=10,100,200,300,400,500,ylabel near ticks,scaled y ticks=base 10:4,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_t=50$

MSE deviation from true error for Gaussian distributions with respect to source sample size. The Bayes error is fixed at 0.2 in all figures.

We show in Fig. <ref> the MSE deviation with respect to the size of target samples for dimensions 2, 3, and 5.
First column corresponds to results when using source datasets of size $n_{s}=50$ and second row shows results for $n_{s}=200$. In addition to the enhanced performance of the TL-based BEE estimator with the increasing availability of target data, we can clearly see that the MSE deviation from the true error asymptotically converges to comparable values. This convergence is observed for all relatedness levels. These results match the observed behaviour with QDA.

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,

legend columns=-1,
legend entries=$\left|\alpha\right|=0.1$,$\left|\alpha\right|=0.3$,$\left|\alpha\right|=0.5$,$\left|\alpha\right|=0.7$,$\left|\alpha\right|=0.9$,$\left|\alpha\right|=0.95$,
legend to name=named7,
legend style=nodes=right,legend style=nodes=scale=0.8, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_s=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE,ylabel near ticks,
xlabel near ticks,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=3$, $n_s=50$

[ width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE,
xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,scaled y ticks=base 10:3,
xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=5$, $n_s=50$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE,
xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,

xlabel near ticks,
legend style=nodes=right,legend style=nodes=scale=1, transform shape,
grid style=dotted,

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

plot coordinates

$d=2$, $n_s=200$

[width=1,font=,cycle list name=colors_qda,
xlabel= $n_t$, ylabel= MSE,
xmin=5,xmax=50,xtick=5,10,20,30,40,50,ylabel near ticks,
xlabel near ticks,

legend style=nodes=right,legend style=nodes=scale=1, transform shape,