# The Pearson Bayes factor: An analytic formula for computing evidential value
from minimal summary statistics

Thomas J. Faulkenberry
Department of Psychological Sciences, Tarleton State University, Stephenville,
TX, USA Corresponding author<EMAIL_ADDRESS>
###### Abstract

In Bayesian hypothesis testing, evidence for a statistical model is quantified
by the Bayes factor, which represents the relative likelihood of observed data
under that model compared to another competing model. In general, computing
Bayes factors is difficult, as computing the marginal likelihood of data under
a given model requires integrating over a prior distribution of model
parameters. This paper builds on prior work that uses the BIC to compute
approximate Bayes factors directly from summary statistics from common
experimental designs (i.e., $t$-test and analysis of variance). Here, I
capitalize on a particular choice of prior distribution that allows the Bayes
factor to be expressed without integral representation (i.e., an analytic
expression), leading to a relatively simple formula (the Pearson Bayes factor)
that requires only minimal summary statistics commonly reported in scientific
papers, such as the $t$ or $F$ score and the degrees of freedom. This provides
applied researchers with the ability to compute exact Bayes factors from
minimal summary data, and thus easily assess the evidential value of any data
for which these summary statistics are provided, even when the original data
is not available.

Keywords: Bayesian statistics; Bayes factor; Pearson Type VI distribution;
Summary statistics; $t$-test; analysis of variance.

Word count: 5,229

Across many fields of study, the hypothesis test remains one of the primary
tools used by researchers to distinguish signal from noise in empirical
observations. Two of the most commonly used tools are the $t$-test and the
analysis of variance. Though classically considered within a frequentist
paradigm, Bayesian versions of these tests are becoming increasingly popular.
Bayesian hypothesis testing makes use of the Bayes factor (Kass and Raftery,,
1995), which represents the relative likelihood of observed data under one
model – say, a model hypothesizing a nonzero group difference – compared to a
competing model (e.g., a model hypothesizing a null difference). Bayes factors
have many desirable properties for inference, but their difficulty of
computation has likely hindered their widespread adoption across the sciences.

One method that has been proposed for computing Bayes factors is the BIC
method (Kass and Raftery,, 1995; Wagenmakers,, 2007), which can be even be
used when only summary statistics are available (Faulkenberry,, 2018;
Faulkenberry, 2019a, ). This popular and well-cited method only requires the
use of a pocket calculator, but it suffers from one primary drawback – it is
an approximate method whose performance suffers in small samples. In this
paper, I will present a new formula for computing Bayes factors from minimal
summary statistics – the Pearson Bayes Factor (PBF). This formula is based on
the work of Min Wang and colleagues (Wang and Sun,, 2014; Wang and Liu,, 2016)
and retains the easy-to-calculate nature of the BIC method while providing the
user with an exact Bayes factor computation, regardless of sample size. The
formula is summarized in the following theorem:

######  Theorem 1.

Given an analysis of variance summary reported in the form $F(x,y)$, where $x$
equals the between-treatments degrees of freedom and $y$ equals the residual
(within-treatments) degrees of freedom, the Bayes factor can be expressed
exactly as

${\operatorname{PBF}}_{10}=\frac{\Gamma\Bigl{(}\frac{x}{2}+\alpha+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{x+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{y}{y+xF}\Biggr{)}^{\alpha-\frac{y}{2}+1}.$

The paper is roughly organized as follows. After first providing some general
statistical background, I will review some of the recent methods used to
compute Bayes factors for $t$-tests and analysis of variance designs and
illustrate the various problems associated with their use. Then, I will
outline the development of the Pearson Bayes factor and provide examples of
its computation. Finally, I will present a simulation study which shows that
the Pearson Bayes factor compares favorably against other well-known Bayes
factors, and in some cases, it has properties which deem it a better choice to
use.

##  1 Background

###  1.1 $t$-test

Originally devised by William Sealy Gosset, who published his work under the
pseudonym Student, (1908), the $t$-test uses the mean and variance of two
observed samples to detect differences in the underlying population means. The
$t$-test works by first assuming a null hypothesis $\mathcal{H}_{0}$, and then
calculating a $t$-score, which indexes the likelihood of obtaining some sample
of observed data under the null hypothesis. If this probability is small, we
reject the null hypothesis $\mathcal{H}_{0}$ in favor of some alternative
hypothesis $\mathcal{H}_{1}$.

Specifically, let us consider the following general setup. Let $x_{ij}$ denote
measurements for the $i^{th}$ participant ($i=1,\dots,N$) in the $j^{th}$
experimental condition ($j=1,2$). Further, we assume that the $x_{ij}$ are
drawn from independent and normally distributed populations with mean
$\mu_{j}$ and variance $\sigma^{2}$. Then we can test the hypotheses

$\mathcal{H}_{0}:\mu_{1}=\mu_{2}\text{ versus
}\mathcal{H}_{1}:\mu_{1}\neq\mu_{2}$

by computing a test statistic

$t=\frac{\overline{x}_{1}-\overline{x}_{2}}{\hat{\sigma}_{p}/\sqrt{N_{\delta}}}$

whose components are computed as follows. $\overline{x}_{j}$ represents the
sample mean of the measurements in group $j$. $\hat{\sigma}_{p}$ is the
“pooled” estimate of $\sigma$, defined by the relationship

$\hat{\sigma}^{2}_{p}=\frac{\hat{\sigma}_{1}^{2}(N_{1}-1)+\hat{\sigma}_{2}^{2}(N_{2}-1)}{N_{1}+N_{2}-2},$

where each $\hat{\sigma}_{j}$ is the sample standard deviation of the
measurements in group $j$. Finally, $N_{\delta}=(1/N_{1}+1/N_{2})^{-1}$; this
is often called the effective sample size for the experiment.

Under the null hypothesis $\mathcal{H}_{0}$, the distribution of these $t$
scores is well known as Student’s $t$ distribution, with density function

$f_{\nu}(x)=\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{\nu\pi}\Gamma\left(\frac{\nu}{2}\right)}\left(1+\frac{x^{2}}{\nu}\right)^{-\frac{\nu+1}{2}},$

where $\nu$ represents degrees of freedom of the test and
$x\in(-\infty,\infty)$. The cumulative distribution function
$F_{\nu}(x)=\int_{-\infty}^{x}f_{\nu}(u)du$ can then be used to index the
probability of observing data at least as extreme as that which we observed
under the null hypothesis $\mathcal{H}_{0}$. Specifically, for an observed
$t$-statistic $t_{\text{obs}}$, we compute
$p(|x|>t_{\text{obs}})=1-F_{\nu}(t_{\text{obs}})$, a quantity commonly known
as a $p$-value. If this probability is small (say, less than 5%), then one may
decide to reject $\mathcal{H}_{0}$ in favor of $\mathcal{H}_{1}$ and conclude
that $\mu_{1}\neq\mu_{2}$, thus implying that the two populations from which
we sampled are indeed different.

###  1.2 Analysis of variance

In 1925, Sir Ronald Fisher introduced the analysis of variance (Fisher,,
1925). As a generalization of the $t$-test procedure that works for comparing
multiple population means, the analysis of variance (or ANOVA) has become one
of the most enduring techniques of hypothesis testing in the experimental
sciences. In its simplest form, the analysis of variance is designed to test
for differences among multiple group means, making it useful for a wide
variety of scientific investigations. Since its introduction, the analysis of
variance has become widespread in its popularity and a core topic in most
introductory statistics textbooks. In fact, its use is so ubiquitous in the
psychological and behavioral sciences that Rouder et al. (Rouder et al.,,
2016) referred to the analysis of variance as the “workhorse” of experimental
psychology.

Roughly, the analysis of variance works by partitioning the total variance in
a set of observed data $\mathcal{D}$ into two sources: the variance between
experimental treatment groups, and the residual, or left over, variance. Then,
one calculates an $F$ statistic, which is defined as the ratio of the between-
groups variance to the residual variance. Like the $t$-test above, inference
about differences between the treatment groups is then performed by
quantifying the likelihood of the observed data $\mathcal{D}$ under a null
hypothesis $\mathcal{H}_{0}$. Specifically, this is done by computing the
probability of obtaining the observed $F$ statistic (or greater) under
$\mathcal{H}_{0}$. If this $p$-value is small, this indicates that the data
$\mathcal{D}$ are rare under $\mathcal{H}_{0}$, so the researcher may reject
$\mathcal{H}_{0}$ in favor of an alternative hypothesis $\mathcal{H}_{1}$
which posits that there are significant differences among the treatment
groups.

###  1.3 Alternatives to $p$-values?

Despite the popularity of the previously-described approaches to testing group
differences, there have been many recent criticisms against their use, and
more generally, against null hypothesis significance testing (see
Wagenmakers,, 2007)). In fact, the American Statistical Association has
recently recommended against the use of $p$-values and significance testing
for scientific inference (Wasserstein and Lazar,, 2016). One alternative that
has been recommended in its place is the Bayes factor (Jeffreys,, 1961;
Raftery,, 1995), which indexes the extent to which the observed data
$\mathcal{D}$ are more likely under one hypothesis (e.g., $\mathcal{H}_{1}$)
than another (e.g., $\mathcal{H}_{0}$). The Bayes factor has several
advantages over the $p$-value as a tool for inference, including the ability
to index support for either the null hypothesis $\mathcal{H}_{0}$ or the
alternative hypothesis $\mathcal{H}_{1}$. However, the widespread adoption of
Bayes factors has likely been prevented due to the general difficulty of their
computation.

In the next section, I will give a conceptual definition of Bayes factors and
outline several approaches to computing Bayes factors for $t$-tests and
ANOVAs.

##  2 Bayes factors

###  2.1 Conceptual definition

To introduce the Bayes factor, let us first recall Bayes’ theorem. Here, we
will relax our notation and let $p$ represent both prior probabilities and
likelihoods. In later sections, we will denote priors and likelihoods by the
more conventional notations $\pi$ and $f$, respectively. Bayes’ theorem states

$\underbrace{p(\mathcal{H}\mid\mathcal{D})}_{\begin{subarray}{c}\text{Posterior
beliefs}\\\ \text{about
model}\end{subarray}}=\underbrace{p(\mathcal{H})}_{\begin{subarray}{c}\text{Prior
beliefs}\\\ \text{about
model}\end{subarray}}\times\underbrace{\frac{p(\mathcal{D}\mid\mathcal{H})}{p(\mathcal{D})}}_{\text{predictive
updating factor}}.$ (1)

which tells us that our posterior belief about model $\mathcal{H}$ can be
found by taking our _prior_ belief about $\mathcal{H}$ and multiplying it by
the fraction $p(\mathcal{D}\mid\mathcal{H})/p(\mathcal{D})$. This fraction is
a _predictive updating factor_ in the sense that it indexes how well our data
is predicted by $\mathcal{H}$ compared to how well it is predicted by all
other possible models. In hypothesis testing, our goal is to directly compare
the predictive adequacy of two models $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$.
To this end, we can use Equation 1 to derive the following:

$\displaystyle\frac{p(\mathcal{H}_{0}\mid\mathcal{D})}{p(\mathcal{H}_{1}\mid\mathcal{D})}$
$\displaystyle=\frac{p(\mathcal{H}_{0})\cdot\frac{p(\mathcal{D}\mid\mathcal{H}_{0})}{p(\mathcal{D})}}{p(\mathcal{H}_{1})\cdot\frac{p(\mathcal{D}\mid\mathcal{H}_{1})}{p(\mathcal{D})}}$

which we can simplify to

$\underbrace{\frac{p(\mathcal{H}_{0}\mid\mathcal{D})}{p(\mathcal{H}_{1}\mid\mathcal{D})}}_{\begin{subarray}{c}\text{posterior
beliefs}\\\ \text{about
models}\end{subarray}}=\underbrace{\frac{p(\mathcal{H}_{0})}{p(\mathcal{H}_{1})}}_{\begin{subarray}{c}\text{prior
beliefs}\\\ \text{about
models}\end{subarray}}\times\underbrace{\frac{p(\mathcal{D}\mid\mathcal{H}_{0})}{p(\mathcal{D}\mid\mathcal{H}_{1})}}_{\text{predictive
updating factor}}.$ (2)

The predictive updating factor

$\operatorname{BF}_{01}=\frac{p(\mathcal{D}\mid\mathcal{H}_{0})}{p(\mathcal{D}\mid\mathcal{H}_{1})}$

is called the _Bayes factor_ , and as we can see in Equation (2), it
represents the factor by which our relative belief in $\mathcal{H}_{0}$ over
$\mathcal{H}_{1}$ can be updated after observing data. Note that we can easily
express our relative belief in $\mathcal{H}_{1}$ over $\mathcal{H}_{0}$ by
taking reciprocals; that is,
$\operatorname{BF}_{10}=1/\operatorname{BF}_{01}$.

One interesting consequence of Equation 2 is that we can use the Bayes factor
to compute the posterior probability of $\mathcal{H}_{0}$ as a function of the
prior model probabilities. To see this, consider the following. If we solve
Equation 2 for the posterior probability $p(\mathcal{H}_{0}\mid\mathcal{D})$
and then use Bayes’ theorem, we see

$\displaystyle p(\mathcal{H}_{0}\mid\mathcal{D})$
$\displaystyle=\operatorname{BF}_{01}\cdot\frac{p(\mathcal{H}_{0})}{p(\mathcal{H}_{1})}\cdot
p(\mathcal{H}_{1}\mid\mathcal{D})$
$\displaystyle=\frac{\operatorname{BF}_{01}\cdot p(\mathcal{H}_{0})\cdot
p(\mathcal{D}\mid\mathcal{H}_{1})\cdot
p(\mathcal{H}_{1})}{p(\mathcal{H}_{1})\cdot p(\mathcal{D})}$
$\displaystyle=\frac{\operatorname{BF}_{01}\cdot p(\mathcal{H}_{0})\cdot
p(\mathcal{D}\mid\mathcal{H}_{1})}{p(\mathcal{D}\mid\mathcal{H}_{0})\cdot
p(\mathcal{H}_{0})+p(\mathcal{D}\mid\mathcal{H}_{1})\cdot
p(\mathcal{H}_{1})}.$

Dividing both numerator and denominator by the marginal likelihood
$p(\mathcal{D}\mid\mathcal{H}_{1})$ gives us

$p(\mathcal{H}_{0}\mid\mathcal{D})=\frac{\operatorname{BF}_{01}\cdot
p(\mathcal{H}_{0})}{\operatorname{BF}_{01}\cdot
p(\mathcal{H}_{0})+p(\mathcal{H}_{1})}.$

Since $\operatorname{BF}_{10}=1/\operatorname{BF}_{01}$, it can then be shown
similarly that

$p(\mathcal{H}_{1}\mid\mathcal{D})=\frac{\operatorname{BF}_{10}\cdot
p(\mathcal{H}_{1})}{\operatorname{BF}_{10}\cdot
p(\mathcal{H}_{1})+p(\mathcal{H}_{0})}.$

In practice, researchers often assume both models are a priori equally likely,
and thus set both $p(\mathcal{H}_{0})=p(\mathcal{H}_{1})=0.5$. In this case,
we obtain the simplified forms

$p(\mathcal{H}_{0}\mid\mathcal{D})=\frac{\operatorname{BF}_{01}}{\operatorname{BF}_{01}+1},\hskip
28.45274ptp(\mathcal{H}_{1}\mid\mathcal{D})=\frac{\operatorname{BF}_{10}}{\operatorname{BF}_{10}+1}.$
(3)

To motivate the work in this paper, let us consider the general problem of
computing Bayes factors. From the conceptual definition given in Equation 2,
we see that the Bayes factor $\operatorname{BF}_{01}$ is the ratio of
likelihoods for $\mathcal{D}$ under $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$,
respectively. Technically, they are marginal likelihoods, where the likelihood
function
$f_{i}(\mathcal{D}\mid\theta_{i})=f(\mathcal{D}\mid\theta_{i},\mathcal{H}_{i})$
($i=0,1$) is integrated over the prior distribution
$\pi_{i}(\theta_{i})=\pi(\theta_{i},\mathcal{H}_{i})$ for the vector of model
parameters $\theta_{i}$ associated with model $\mathcal{H}_{i}$. That is,

$\operatorname{BF}_{01}=\frac{m_{0}(\mathcal{D})}{m_{1}(\mathcal{D})}$

where

$m_{0}(\mathcal{D})=\int
f_{0}(\mathcal{D}\mid\theta_{0})\pi_{0}(\theta_{0})d\theta_{0}$

and

$m_{1}(\mathcal{D})=\int
f_{1}(\mathcal{D}\mid\theta_{1})\pi_{1}(\theta_{1})d\theta_{1}.$

In general, computing these marginal likelihoods can be difficult. One reason
is that as sample size increases, the likelihood becomes highly peaked at its
maximum. Without prior knowledge of the location of this maximum, numerical
integration methods can fail to converge around this highly peaked area.
Related to this is the problem of dimensionality: the parameter vector
$\theta_{i}$ is often of high dimension. Combined with the peaked nature of
the integrand, numerical methods often have difficulty – a problem that
Thiele, Rouder, and Haaf refer to as “finding a needle in a haystack” (Thiele
et al.,, 2017). Thus, when possible, it is desirable to make simplifying
assumptions which lead to analytic solutions (i.e., without integral
representation) for these marginal likelihoods.

###  2.2 BIC approximation

The most popular approach to this problem is the BIC approximation, described
in detail by various papers of Kass and Raftery (e.g., Raftery,, 1995; Kass
and Raftery,, 1995) and introduced to the applied psychological science
community by Wagenmakers, (2007) and Masson, (2011). The basic idea is to
construct a second-order Taylor approximation of $\ln m_{i}(\mathcal{D})=\ln
p(\mathcal{D}\mid\theta_{i},\mathcal{H}_{i})$ about the posterior mode for
$\theta_{i}$. Of course, by removing all terms beyond the second-derivative,
the result is an approximation of the log marginal likelihood. However,
Raftery, (1995) pointed out that with a certain choice of noninformative prior
on $\theta_{i}$ (the unit information prior), the order of the approximation
is $\mathcal{O}(1/\sqrt{n})$. This results in the approximation

$\ln m_{i}(\mathcal{D})=\ln
p(\mathcal{D}\mid\hat{\theta}_{i},\mathcal{H}_{i})-\frac{k_{i}}{2}\ln
n+\mathcal{O}\Bigl{(}\frac{1}{\sqrt{n}}\Bigr{)}$ (4)

where $\hat{\theta}_{i}$ is the maximum likelihood estimate for $\theta_{i}$
under $\mathcal{H}_{i}$, $k_{i}$ is the number of parameters in
$\mathcal{H}_{i}$, and $n$ is the total number of observations in
$\mathcal{D}$. Equation 4 has a fortuitous relation to the Bayesian
information criterion (BIC) of Schwarz, (1978):

$\operatorname{BIC}(\mathcal{H}_{i})=-2\ln L_{i}+k_{i}\ln n,$ (5)

where $L$ is the maximum likelihood estimate for model $\mathcal{H}_{i}$.
Combining Equations 4 and 5, we have

$\ln
m_{i}(\mathcal{D})=-\frac{1}{2}\operatorname{BIC}(\mathcal{H}_{i})+\mathcal{O}\Bigl{(}\frac{1}{\sqrt{n}}\Bigr{)},$

or equivalently

$m_{i}(\mathcal{D})\approx\exp\Bigl{(}-\frac{1}{2}\operatorname{BIC}(\mathcal{H}_{i})\Bigr{)}.$

It is then easy to derive the following expression for the BIC approximation
to the Bayes factor:

$\displaystyle\operatorname{BF}_{01}=\frac{m_{0}(\mathcal{D})}{m_{1}(\mathcal{D})}$
$\displaystyle\approx\frac{\exp\Bigl{(}-\frac{1}{2}\operatorname{BIC}(\mathcal{H}_{0})\Bigr{)}}{\exp\Bigl{(}-\frac{1}{2}\operatorname{BIC}(\mathcal{H}_{1})\Bigr{)}}$
(6)
$\displaystyle=\exp\Biggl{(}\frac{\operatorname{BIC}(\mathcal{H}_{1})-\operatorname{BIC}(\mathcal{H}_{0})}{2}\Biggr{)}.$

Compared to computing marginal likelihoods by integrating the prior-weighted
likelihoods, Equation 6 requires only that we know the BIC values for models
$\mathcal{H}_{0}$ and $\mathcal{H}_{1}$. For ANOVA models, the calculation is
straightforward. In the ANOVA context, the BIC can be calculated (Raftery,,
1995) as

$\displaystyle BIC$ $\displaystyle=n\ln(1-R^{2})+k\ln n$
$\displaystyle=n\ln\Biggl{(}\frac{SSR}{SST}\Biggr{)}+k\ln n$

where SST represents the total sum of squares for the model (i.e.,
$\sum_{i}\sum_{j}(Y_{ij}-\overline{Y}_{\cdot\cdot})^{2})$, SSR represents the
residual sum of squares left over after accounting for the model’s postulated
treatment effects (i.e., $SST-SSA$, where
$SSA=r\sum_{i}(\overline{Y}_{i\cdot}-\overline{Y}_{\cdot\cdot})^{2})$, and $k$
equals the number of distinct treatment conditions in the model.

The BIC approach from Equation 6 requires that we have enough “raw” data
available in order to compute $SSR$ and $SST$. One improvement to Equation 6
is to remove this dependency and make the Bayes factor computation accessible
with only summary statistics. Faulkenberry, (2018) did exactly this, further
refining the BIC approximation to produce a formula that can be used by the
applied researcher who has access to the commonly reported summary statistics
from the ANOVA, usually reported in the standard form $F(x,y)$, where $F$ is
the observed $F$-statistic, $x$ is the degrees of freedom between treatments,
and $y$ is the residual (within treatments) degrees of freedom. This results
in the simple closed form expression

$\operatorname{BF}_{01}\approx\sqrt{n^{x}\Biggl{(}1+\frac{Fx}{y}\Biggr{)}^{-n}}.$
(7)

The advantage of Faulkenberry’s formula in Equation 7 is that researchers may
compute approximations to Bayes factors from summary statistics only, whereas
the approach of Wagenmakers, (2007) and Masson, (2011) that uses the BIC
approximation in Equation 6 requires knowing the underlying sum of squares
terms in the ANOVA calculation. Thus, Equation 7 can be useful to researchers
who are trying to assess evidential value of observed data from published
summary statistics, especially where raw data (or even ANOVA summary tables)
are not available (Faulkenberry, 2019b, ).

The Bayes factor approximation in Equation 7 can be easily adapted to use with
the $t$-test by simply noting that $t^{2}=F$ and $x=1$. Thus, for a $t$ test
with degrees of freedom equal to $\nu$, we have the approximation

$\operatorname{BF}_{01}\approx\sqrt{n\Biggl{(}1+\frac{t^{2}}{\nu}\Biggr{)}^{-n}}.$

.

Before proceeding further, it will be instructive to consider a simple example
of computing the BIC Bayes factor for an ANOVA model.

######  Example 1.

Consider the following data representing 6 measurements in each of 3 treatment
groups:

Treatment 1 | Treatment 2 | Treatment 3
---|---|---
2 | 5 | 8
3 | 9 | 6
8 | 10 | 12
6 | 13 | 11
5 | 8 | 11
6 | 9 | 12

From here, it is a routine exercise to complete the following ANOVA summary
table:

Source | $SS$ | $df$ | $MS$ | $F$
---|---|---|---|---
Treatments | 84 | 2 | 42 | 7.16
Residual | 88 | 15 | 5.87 |
Total | 172 | 17 |  |

To use Equation 6, we need to compute the BIC of $\mathcal{H}_{0}$ and
$\mathcal{H}_{1}$. First, we compute

$\displaystyle\operatorname{BIC}(\mathcal{H}_{1})$
$\displaystyle=n\ln\Biggl{(}\frac{SSR}{SST}\Biggr{)}+k\ln n$
$\displaystyle=18\ln\Biggl{(}\frac{88}{172}\Biggr{)}+3\ln(18)$
$\displaystyle=-3.392.$

Under $\mathcal{H}_{0}$, we note that all data are predicted by the grand mean
$\mu$ only. Thus, all variation to be considered as residual (i.e.,
$SSR=172$), and there is $k=1$ treatment condition. This gives us

$\displaystyle\operatorname{BIC}(\mathcal{H}_{0})$
$\displaystyle=18\ln\Biggl{(}\frac{172}{172}\Biggr{)}+1\ln(18)$
$\displaystyle=2.890.$

Thus,

$\displaystyle\operatorname{BF}_{01}$
$\displaystyle\approx\exp\Biggl{(}\frac{\operatorname{BIC}(\mathcal{H}_{1})-\operatorname{BIC}(\mathcal{H}_{0})}{2}\Biggr{)}$
$\displaystyle=\exp\Biggl{(}\frac{-3.392-2.890}{2}\Biggr{)}$
$\displaystyle=0.0432$

Taking the reciprocal, we see

$\operatorname{BF}_{10}=\frac{1}{BF_{01}}\approx\frac{1}{0.0432}=23.15$

indicating that the observed data $\mathcal{D}$ are 23.15 times more likely
under $\mathcal{H}_{1}$ than $\mathcal{H}_{0}$. Note that this also agrees
with the formula using only the summary statistics from the ANOVA table:
$F(2,15)=7.16$. Taking from this expression $x=2$ and $y=15$, we evaluate
Equation 7 as

$\displaystyle\operatorname{BF}_{01}$
$\displaystyle\approx\sqrt{n^{x}\Biggl{(}1+\frac{Fx}{y}\Biggr{)}^{-n}}$
$\displaystyle=\sqrt{18^{2}\Biggl{(}1+\frac{7.16\cdot 2}{15}\Biggr{)}^{-18}}$
$\displaystyle=0.0432,$

thus demonstrating the equivalence of the approaches from Equations 6 and 7.

Though Equations 6 and 7 provide simple ways to compute Bayes factors from
ANOVA summaries, it is important to remember that the BIC formulas represent
an approximation to the Bayes factor. This fact is underscored by the
following computation. Sellke et al., (2001) showed that under a reasonable
class of prior distributions for $p$-values, an upper bound for the Bayes
factor can be computed directly from the $p$-value as

$\operatorname{BF}_{10}\leq-\frac{1}{e\cdot p\ln(p)}.$

Applying this bound to our previous example, we first note that for
$F(2,15)=7.16$, the associated $p$-value is $p=0.0066$. This gives us the
upper bound

$\displaystyle\operatorname{BF}_{10}$ $\displaystyle\leq-\frac{1}{e\cdot
0.0066\cdot\ln(0.0066)}$ $\displaystyle=11.10.$

Thus, it is clear that our BIC Bayes factor of 23.15 is quite an overestimate
of the actual Bayes factor, so our desire for exact Bayes factors is quite
warranted.

##  3 Analytic methods for computing Bayes factors

Since its introduction by Kass and Raftery, (1995), the BIC Bayes factor has
been a popular method for estimating evidential value from empirical data. One
of the attractive features of the BIC Bayes factor is that it can be computed
with a simple calculator. This makes it easy to implement not only for the
researcher, but also in the context of a beginning statistics course. However,
it is a large-sample approximation, and the discussion above calls into
question its computational stability for smaller examples, especially those
which might be used in a traditional statistics course.

Thus, it is desirable to find analytic methods – that is, solutions without
the need for approximating integrals – to replace the BIC approximation. In
other words, we want a method which removes the integral representation, but
unlike the BIC, gives an exact value for the Bayes factor. In the following
sections, I will discuss some historical approaches to developing these
methods for the $t$-test and the analysis of variance, followed by a
discussion of their limitations.

###  3.1 Analytic Bayes factors for the $t$-test

One approach to developing a Bayesian $t$-test with an analytic solution was
provided by Gönen et al., (2005). Gönen et al. reparameterized
$\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ in terms of “effect size”
$\delta=\mu_{1}-\mu_{2}$. Their goal was to test

$\mathcal{H}_{0}:\delta=0\text{ versus }\mathcal{H}_{1}:\delta\neq 0,$

using a parameterization $(\mu=(\mu_{1}+\mu_{2})/2,\delta,\sigma^{2})$.
Following Jeffreys, (1961), they placed a noninformative prior on the mean and
variance:

$\pi(\mu,\sigma^{2})\propto 1/\sigma^{2}$

and assumed, under $\mathcal{H}_{1}$, that $\delta$ was drawn from a normal
distribution centered at 0 with a variance parameter $\sigma^{2}_{a}$. This
variance term serves as as hyperparameter that needed to be specified in
advance by the analyst. Given this setup, Gönen et al., (2005) showed

$\operatorname{BF}_{10}=\Biggl{(}\frac{1+\frac{t^{2}}{\nu}}{1+\frac{t^{2}}{\nu(1+n_{\delta}\sigma^{2}_{a})}}\Biggr{)}^{\frac{\nu+1}{2}}(1+n_{\delta}\sigma^{2}_{a})^{-\frac{1}{2}}.$
(8)

Certainly, Equation 8 meets our requirements of (1) being an exact Bayes
factor, and (2) being solvable by means of a calculator, but it does have some
drawbacks. While the ability for the researcher to specify the prior variance
of the effect size under $\mathcal{H}_{1}$ may be initially appealing, it is
fraught with some unfortunate consequences.

Issue 1: It is a common practice to specify a large value for
$\sigma^{2}_{a}$, as such a choice would minimize the information provided by
the prior and thus potentially better serve as an “objective” choice of prior.
However, this causes a problem. Consider a fixed data summary; suppose we
observed $t=3.00$ with $N_{1}=N_{2}=20$. Applying Equation 8 with increasing
values of $\sigma^{2}_{a}$, we see a troubling pattern, shown in the following
table:

$\sigma^{2}_{a}$ | 1 | 10 | 100 | 1000 | 10000
---|---|---|---|---|---
$\operatorname{BF}_{10}$ | 8.99 | 6.01 | 2.08 | 0.67 | 0.21

Indeed, this table indicates that the Gönen Bayes factor exhibits Bartlett’s
paradox; that is, for fixed information from the data, the Bayes factor tends
to 0 (and thus always favors $\mathcal{H}_{0}$) with increasing prior variance
(see also Wang and Liu,, 2016).

Issue 2: Fix a set of observed data (samples) and suppose these observed data
are always generated under $\mathcal{H}_{1}$. Then as $t$ increases with fixed
$\sigma^{2}_{a}$, the Bayes factor in Equation 8 tends to
$(1+n_{\delta}\sigma^{2}_{a})^{\eta/2}$ (see Figure 1). However, since the
samples are drawn from $\mathcal{H}_{1}$, we would expect that as $t$
increases without bound, the Bayes factor should also increase without bound.
The unexpected asymptote in Figure 1 is referred to as the information
paradox.

Figure 1: Information paradox: with fixed data under $\mathcal{H}_{1}$, the
Gönen Bayes factor is bounded with increasing $t$.

###  3.2 Analytic Bayes factors for ANOVA

Similar to the Gönen et al., (2005) $t$-test, it is also possible to calculate
analytic solutions for Bayes factors with ANOVA summaries. One approach is to
use Zellner’s $g$-prior approach, which was originally developed for the
context of linear regression models (Zellner,, 1986). In the $g$-prior
approach, the standardized regression coefficients $\beta$ are given a
multivariate normal prior centered at 0 with a scaled covariance matrix that
is proportional to the inverse Fisher information matrix for $\beta$. This
scaling factor, denoted $g$, must be chosen by the analyst. Given the
construction of the prior, $g$ serves as a scale factor for the analyst to
weigh the relative “importance” of the prior and the sample in the resulting
inference. If $g=1$, then the prior and the sample have the same importance in
the inference. If $g=2$, then the prior has half the importance as the sample
in the inference. This pattern continues; if $g=n$, then the prior has
$1/n^{\text{th}}$ the importance as the sample in the inference.

The $g$-prior Bayes factor can be computed as

$\operatorname{BF}_{10}=(1+g)^{\frac{n-k-1}{2}}\bigl{(}1+g(1-R^{2})\bigr{)}^{-\frac{n-1}{2}}$
(9)

where $n$ equals the number of observations, $k$ equals the number of
predictors in the regression model, and $R^{2}$ equals the coefficient of
determination for the regression model. The analyst must choose the value of
$g$, as discussed above. Unfortunately, this leads to the same issues
discussed with the Gönen et al., (2005) above. First, assuming fixed $n$ and
$k$, increasing values of $g$ in Equation 9 will cause the Bayes factor to
approach 0 (i.e., Bartlett’s paradox). Second, increasing values of $R^{2}$
from data drawn under $\mathcal{H}_{1}$ will only result in Bayes factors that
approach an asymptote, $(1+g)^{(n-k-1)/2}$ (i.e., the information paradox).

In light of these paradoxes (both the Gönen and Zellner approaches), a common
approach for resolution is to place a prior distribution on the relevant
scaling parameter (e.g., the $g$ in the Zellner Bayes factor). Two solutions
have been recently proposed: (1) the hyper-$g$ approach (Liang et al.,, 2008);
and (2) the JZS approach (Rouder et al.,, 2012). Both of these approaches
place prior priors on the scale factor $g$ and thus avoid the paradoxes
described above. However, the drawback to both of these approaches is that the
resulting Bayes factors have integral representation. That is, their
computation requires evaluation of at least a one-dimensional integral. With
software, this is no problem, but these approaches are inaccessible to
beginners (especially students) with no calculus background or programming
knowledge.

###  3.3 Analytic Bayes factors without integral representation

Since its introduction by Kass and Raftery, (1995), the BIC Bayes factor has
been a popular method for estimating evidential value from empirical data.
However, since it is based on a large sample approximation, its utility is
limited, particularly in a small sample context. It is especially limited for
use in a teaching context, where sample sizes for hand-worked examples tend to
be small. Thus, it is desirable to find a calculation method to replace the
BIC approximation. Specifically, we want a method which, like the BIC, removes
the integral representation, but unlike the BIC, gives an exact value for the
Bayes factor.

Recent work by Wang and Sun, (2014) has made possible such calculation. In
what follows, we consider our data to be a balanced set of observations
$\mathcal{D}=Y_{ij}$ where $i=1,\dots,p$ represents the experimental unit
(i.e., the treatment group) and $j=1,\dots,r$ represents the specific
replicate within each unit. This gives a total of $n=pr$ observations. On
these data, we place a random effects linear model:

$Y_{ij}=\mu+a_{i}+\varepsilon_{ij},$

assuming that $a_{i}\sim\mathcal{N}(0,\sigma^{2}_{a})$ and
$\varepsilon\sim\mathcal{N}(0,\sigma^{2})$. In this context, the main goal of
analysis of variance is to test whether the random effects term $a_{i}$ is
identically 0. We do this by restricting the variability of the random effects
term. That is, we want to test

$\mathcal{H}_{0}:\sigma^{2}_{a}=0\text{ versus
}\mathcal{H}_{1}:\sigma_{a}^{2}\neq 0.$

Wang and Sun extended the approach of García-Donato and Sun, (2007), who
assumed common, noninformative priors on $\mu$ and $\sigma$ in both
$\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ and considered a proper prior on the
ratio of variance components $\tau=\sigma^{2}_{a}/\sigma^{2}$ under
$\mathcal{H}_{1}$. Under such prior specification, Garcia-Donato and Sun
showed

$BF_{10}=\int_{0}^{\infty}(1+\tau r)^{\frac{1-p}{2}}\Biggl{(}1-\frac{\tau
r}{1+\tau r}\cdot\frac{SSA}{SST}\Biggr{)}^{\frac{1-n}{2}}\cdot\pi(\tau)d\tau$
(10)

To use Equation 10, the user must place a prior distribution on the variance
components $\tau$. Motivated by the work of Maruyama, (2009), Wang and Sun
used a Pearson Type VI distribution to serve as the prior for $\tau$. The
Pearson Type VI distribution has three parameters: two shape parameters
$\alpha>-1$ and $\beta>-1$ and a scale parameter $\kappa>0$. The density
function for the Pearson Type VI prior is given by

$\pi^{PT}(\tau)=\frac{\kappa(\kappa\tau)^{\beta}(1+\kappa\tau)^{-\alpha-\beta-2}}{\mathcal{B}(\alpha+1,\beta+1)}I_{(0,\infty)}(\tau)$

where $\mathcal{B}(x,y)=\int_{0}^{1}t^{x-1}(1-t)^{y-1}dt$ is the standard Beta
function. Wang and Sun further restricted $\pi^{PT}$ to one parameter
$\alpha\in[-\frac{1}{2},0]$ by taking $\kappa=r$ and
$\beta=\frac{n-p}{2}-\alpha-2$. A plot of this prior can be seen in Figure 2;
in this figure, we take the values $r=6$, $p=3$, and $n=18$ from Example 1,
thus setting $\kappa=6$ and $\beta=\frac{18-3}{2}-\alpha-2$, where $\alpha$
takes specific values $-\frac{1}{2},-\frac{1}{4},-\frac{1}{10},0$. Further,
Figure 2 shows the effect of varying $\alpha$ on the prior distribution for
$\tau$. As we can see, as $\alpha$ increases to 0, $\tau$ becomes more
dispersed and less peaked around the mode. This places more prior mass on
larger treatment effects than we would see for values of $\alpha$ closer to
$-\frac{1}{2}$.

Figure 2: A Pearson Type VI prior for $\tau$, plotted as a function of shape
parameter $\alpha$.

With the above prior specification, Wang and Sun, (2014) proved that the Bayes
factor of Garcia-Donato and Sun (Equation 10) simplifies to an analytic
expression without integral representation:

$\operatorname{BF}_{10}=\frac{\Gamma\Bigl{(}\frac{p}{2}+\alpha+\frac{1}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\frac{n-p}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{n-1}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\cdot\Biggl{(}\frac{SSR}{SST}\Biggr{)}^{\alpha-\frac{n-p-2}{2}},$
(11)

With this equation, the Bayes factor can be computed exactly without any need
for numerical integration routines. Note: though the gamma function is itself
defined by an integral, its exact values can often be found, particularly for
positive integers (i.e., $\Gamma(n)=(n-1)!$ for $n\in\mathbb{N}$). In other
cases, extremely precise approximations can be had with most standard
programming libraries).

######  Example 2.

Let us now compute a Bayes factor using the Wang and Sun, (2014) approach of
Equation 11. We will consider the same data set given in Example 1, so we have
$n=18$ and $p=3$. Equation 11 requires the user to specify
$\alpha\in[-\frac{1}{2},0]$, which governs the shape of the Pearson Type VI
prior on $\tau$, the ratio of variance components $\sigma^{2}_{a}/\sigma^{2}$.
To see the effect of this choice on the resulting Bayes factor, we will
consider two cases on the boundary: $\alpha=0$ and $\alpha=-\frac{1}{2}$.

Case 1: Set $\alpha=0$. Then by Equation 11

$\displaystyle\operatorname{BF}_{10}$
$\displaystyle=\frac{\Gamma(\frac{3}{2}+0+\frac{1}{2})\cdot\Gamma(\frac{18-3}{2})}{\Gamma(\frac{18-1}{2})\cdot\Gamma(0+1)}\Biggl{(}\frac{88}{172}\Biggr{)}^{0-\frac{18-3-2}{2}}$
$\displaystyle=\frac{\Gamma(2)\cdot\Gamma(\frac{15}{2})}{\Gamma(\frac{17}{2})\cdot\Gamma(1)}\Biggl{(}\frac{88}{172}\Biggr{)}^{-\frac{13}{2}}$
$\displaystyle=\frac{1\cdot 1871.254}{14034.41\cdot
1}\Biggr{(}\frac{88}{172}\Biggr{)}^{-\frac{13}{2}}$ $\displaystyle=10.393.$

Case 2: Set $\alpha=-\frac{1}{2}$. Then by Equation 11

$\displaystyle\operatorname{BF}_{10}$
$\displaystyle=\frac{\Gamma(\frac{3}{2}-\frac{1}{2}+\frac{1}{2})\cdot\Gamma(\frac{18-3}{2})}{\Gamma(\frac{18-1}{2})\cdot\Gamma(-\frac{1}{2}+1)}\Biggl{(}\frac{88}{172}\Biggr{)}^{-\frac{1}{2}-\frac{18-3-2}{2}}$
$\displaystyle=\frac{\Gamma(\frac{3}{2})\cdot\Gamma(\frac{15}{2})}{\Gamma(\frac{17}{2})\cdot\Gamma(\frac{1}{2})}\Biggl{(}\frac{88}{172}\Biggr{)}^{-7}$
$\displaystyle=\frac{0.8862269\cdot 1871.254}{14034.41\cdot
1.772454}\Biggr{(}\frac{88}{172}\Biggr{)}^{-7}$ $\displaystyle=7.265.$

Notice that for both values of $\alpha$, the obtained Bayes factor is quite
different from the value obtained from the BIC approximation in Example 1.
However, both are well under the Sellke et al. bound of 11.10 (which we
calculated earlier).

##  4 Analytic Bayes factor from summary statistics

We are now ready to prove the primary new result in this paper; namely, the
development of the Pearson Bayes factor, which is summarized in Theorem 1. The
Pearson Bayes factor provides researchers with an analytic Bayes factor (i.e.,
without integral representation) that can be obtained directly from the
summary statistics typically reported from a $t$-test or analysis of variance.
To do this, we simply need to start with Equation 11 of Wang and Sun, (2014).
As we previously saw, their Bayes factor requires knowing the residual sum of
squares term $SSR$, the total sum of squares $SST$, the total number of
observations $n$, and the number of treatment groups $p$. Additionally, the
user must specify a hyperparameter $\alpha$, which controls the shape of the
prior distribution for $\tau=\sigma^{2}_{a}/\sigma^{2}$, the ratio of variance
components.

To begin the proof, let us first consider two facts. First, $x$ is the
between-treatments degrees of freedom, so $x=p-1$. This implies $p=x+1$.
Second, $y$ is the residual (within-treatments) degrees of freedom, so
$y=n-p$. Putting these two facts together, we have $x+y=(p-1)+(n-p)=n-1$.
Substituting these into Equation 11 gives us

$\displaystyle\operatorname{BF}_{10}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{p}{2}+\alpha+\frac{1}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\frac{n-p}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{n-1}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\cdot\Biggl{(}\frac{SSR}{SST}\Biggr{)}^{\alpha-\frac{n-p-2}{2}}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{x+1}{2}+\alpha+\frac{1}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{x+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\cdot\Biggl{(}\frac{SSR}{SST}\Biggr{)}^{\alpha-\frac{y-2}{2}}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{x}{2}+\alpha+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{x+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\cdot\Biggl{(}\frac{SSR}{SST}\Biggr{)}^{\alpha-\frac{y}{2}+1}.$

The last step is to express the fraction $SSR/SST$ in terms of $F$, $x$, and
$y$. To this end, we note that the reciprocal can be written as

$\frac{SST}{SSR}=\frac{SSR+SSA}{SSR}=1+\frac{SSA}{SSR}.$

By definition, the $F$ ratio is written as
$F=\frac{SSA/x}{SSR/y}=\frac{SSA}{SSR}\cdot\frac{y}{x}$. Combining these gives
us

$\frac{SST}{SSR}=1+\frac{xF}{y}=\frac{y+xF}{y}.$

Reciprocating and substituting back into the equation above puts the Bayes
factor into the form of Theorem 1, thus completing the proof.

######  Example 3.

Again, we consider the same data set as Example 1, where we found
$F(2,15)=7.16$. From these summary statistics, we apply Theorem 1 using
$F=7.16$, $x=2$, and $y=15$. As before, we consider the two boundary cases for
the value of $\alpha$:

Case 1: $\alpha=0$. By Theorem 1, we have

$\displaystyle\operatorname{PBF}_{10}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{2}{2}+0+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{15}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{2+15}{2}\Bigr{)}\cdot\Gamma(0+1)}\Biggl{(}\frac{15}{15+2\cdot
7.16}\Biggr{)}^{0-\frac{15}{2}+1}$ $\displaystyle=\frac{1\cdot
1871.254}{14034.41\cdot 1}(0.5116)^{-6.5}$ $\displaystyle=10.397$

Case 2: $\alpha=-1/2$. Here, we have

$\displaystyle\operatorname{PBF}_{10}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{2}{2}-\frac{1}{2}+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{15}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{2+15}{2}\Bigr{)}\cdot\Gamma(-\frac{1}{2}+1)}\Biggl{(}\frac{15}{15+2\cdot
7.16}\Biggr{)}^{-\frac{1}{2}-\frac{15}{2}+1}$
$\displaystyle=\frac{0.8862269\cdot 1871.254}{14034.41\cdot
1.772454}(0.5116)^{-7}$ $\displaystyle=7.268.$

In both cases, Theorem 1 produces Bayes factors reasonably close to the
answers we produced using the Wang and Sun formula (Equation 11) in Example 2.
Indeed, the inaccuracy results from the inevitable two-digit rounding that we
often see with summary statistics when reported in papers. If instead we had
used $F=7.159$ instead of $F=7.16$, our answers would have matched exactly.

Finally, let us consider the case where $p=2$; that is, we have two treatment
conditions. In this case, the between-treatments degrees of freedom is equal
to 1. Thus, if we apply Theorem 1 with $x=1$, we have the following:

$\displaystyle\operatorname{PBF}_{10}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{x}{2}+\alpha+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{x+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{y}{y+xF}\Biggr{)}^{\alpha-\frac{y}{2}+1}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{1}{2}+\alpha+1\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{1+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{y}{y+F}\Biggr{)}^{\alpha-\frac{y}{2}+1}$
$\displaystyle=\frac{\Gamma\Bigl{(}\alpha+\frac{3}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\frac{y}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{1+y}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{y}{y+F}\Biggr{)}^{(2\alpha-y+2)/2}.$

If we substitute the identity $F=t^{2}$ and replace $y$ with $\nu$ (the
typical symbol used to represent degrees of freedom in a $t$-test), we get the
following:

$\displaystyle\operatorname{PBF}_{10}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{\nu}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\alpha+\frac{3}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{\nu+1}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{\nu}{\nu+t^{2}}\Biggr{)}^{(2\alpha-\nu+2)/2}$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{\nu}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\alpha+\frac{3}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{\nu+1}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}\frac{\nu+t^{2}}{\nu}\Biggr{)}^{(-2\alpha+\nu-2)/2}.$
$\displaystyle=\frac{\Gamma\Bigl{(}\frac{\nu}{2}\Bigr{)}\cdot\Gamma\Bigl{(}\alpha+\frac{3}{2}\Bigr{)}}{\Gamma\Bigl{(}\frac{\nu+1}{2}\Bigr{)}\cdot\Gamma(\alpha+1)}\Biggl{(}1+\frac{t^{2}}{\nu}\Biggr{)}^{(\nu-2\alpha-2)/2}$

which exactly matches the Pearson Bayes factor derived by Wang and Liu, (2016)
for the $t$-test. Thus, our proposed Pearson Bayes factor formula from Theorem
1 serves as a generalization that works for summary statistics from both the
$t$-test and analysis of variance.

##  5 Simulation study

In this final section, I report the results of a simulation study designed to
benchmark the performance of the Pearson Bayes factor against other popular
Bayes factors. In the simulation, I generated random datasets that each
reflected the balanced one-factor designs that we have discussed throughout
this paper. Specifically, data were generated as

$Y_{ij}=\mu+a_{j}+\varepsilon_{ij};\hskip 14.22636pti=1,\dots,r;\hskip
8.53581ptj=1,\dots,p,$

where $\mu$ represents a grand mean, $a_{j}\sim\mathcal{N}(0,\sigma_{a})$
represent each of the $p$ randomly drawn treatment effects, and
$\varepsilon_{ij}\sim\mathcal{N}(0,\sigma^{2})$ represent the normally-
distributed error terms. For convenience I set $p=3$, though similar results
were obtained with other values of $p$ (not reported here). Also, without loss
of generality I set $\mu=0$ and $\sigma=1$. I then systematically varied the
following three components of the model:

  1. 1.

The number of replicates per treatment condition $r$ was set to either $n=10$,
$n=30$, or $n=80$;

  2. 2.

The size of the overall treatment effect was manipulated by setting
$\tau=\sigma_{a}^{2}/\sigma^{2}$ to be either $\tau=0$, $\tau=0.5$, or
$\tau=1$. Note that data generated under the condition $\tau=0$ reflect the
null model $\mathcal{H}_{0}$, whereas data generated under $\tau>0$ reflect
the alternative model $\mathcal{H}_{1}$.

  3. 3.

The prior setting $\alpha$ used in the computation of the Pearson Bayes factor
was set at both ends of its consistency range; that is, $\alpha=-1/2$ and
$\alpha=0$.

For each combination of number of replicates ($r=10,30,80$), treatment effect
($\tau=0,0.5,1.0$), and prior setting $\alpha$ ($\alpha=-0.5,0$), I generated
1000 simulated datasets. For each of the datasets, I computed three different
Bayes factors: the Pearson Bayes factor proposed in this paper, the JZS Bayes
Factor (Rouder et al.,, 2012), and the BIC Bayes factor (Faulkenberry,, 2018).
These Bayes factors were then converted to posterior probabilities via
Equation 3. To compare the performance of the three Bayes factors in the
simulation, I considered three analyses for each simulated dataset: (1) a
visualization of the distribution of posterior probabilities
$p(\mathcal{H}_{0}\mid\mathcal{D})$; (2) a calculation of the proportion of
simulated trials for which the correct model was chosen (i.e., model choice
accuracy); and (3) a calculation of the proportion of simulated trials for
which the Pearson Bayes factor chose the same model as the JZS Bayes factor
(i.e., model choice consistency).

To begin, let us consider the distribution of posterior probabilities
$p(\mathcal{H}_{0}\mid\mathcal{D})$. Boxplots of these posterior probabilities
can be seen in Figure 3. From the figure, we can see that the Pearson Bayes
factor performs well compared to the JZS and BIC Bayes factors, but there are
some distinct differences. First, when $\tau=0$ (and thus, $\mathcal{H}_{0}$
is the correct model) both the Pearson and BIC Bayes factors produce larger
posterior probabilities for $\mathcal{H}_{0}$ than the JZS Bayes factor. When
the prior parameter $\alpha$ is set to $-1/2$, the Pearson Bayes factor
produce a smaller range of posterior probabilities than the BIC Bayes factor.
When $\tau>0$ ($\mathcal{H}_{1}$ is the correct model), the overall advantage
goes to the JZS Bayes factor, though the distribution of posterior
probabilities becomes indistinguishable across methods as $\tau$ increases and
$r$ increases. Note that in the $\tau=0.5$ case (where treatment effect
variance is 1/2 of residual variance), model selection is a bit more
scattered. In this case, the Pearson Bayes factor outperforms (i.e., smaller
range and smaller posterior probability for $\mathcal{H}_{0}$) the BIC Bayes
factor when the prior parameter $\alpha$ is set to 0. This is not surprising,
as the Pearson Type VI prior for $\tau$ is more dispersed when $\alpha=0$,
allowing for a greater prior mass being placed on larger effects.

Figure 3: Results from our simulation. Each boxplot depicts the distribution
of the posterior probability $p(\mathcal{H}_{0}\mid\mathcal{D})$ for 1000
Monte Carlo simulations. White boxes represent posterior probabilities derived
from the Pearson Bayes factor. Medium gray boxes represent posterior
probabilities derived from the JZS Bayes factor. Dark gray boxes represent
posterior probabilities derived from the BIC Bayes factor.

Though the distributions of posterior probabilities tend to follow similar
patterns, it is not clear to what extent the three methods provide the user
with an accurate inference. To compare the accuracy of the three methods, I
calculated model choice accuracy, defined as the proportion of simulated
datasets for which the correct model was chosen. Model choice was defined by
considering $\mathcal{H}_{0}$ to be chosen whenever $\operatorname{BF}_{01}>1$
and $\mathcal{H}_{1}$ to be chosen whenever $\operatorname{BF}_{01}<1$. The
results are displayed in Table 1. From the table, we can see confirmation of
the intuition gained from the distribution of posterior probabilities in
Figure 3. Again, for null effects, both the PBF and the BIC consistently
outperform the JZS Bayes factor (particularly when the prior shape parameter
is set to $\alpha=-1/2$). When $\tau>0$, the the JZS Bayes factor chooses the
correct model ($\mathcal{H}_{1}$) more often than either the PBF or the BIC.
Again, the PBF outperforms BIC when setting $\alpha=0$. Accuracy of both PBF
and BIC increase as the number of replicates per treatment condition
increases.

|  | $\alpha=-1/2$ |  | $\alpha=0$
---|---|---|---|---
|  | PBF | JZS | BIC |  | PBF | JZS | BIC
$\tau=0$: |  |  |  |  |  |  |  |
$r=10$ |  | .980 | .913 | .965 |  | .931 | .886 | .978
$r=30$ |  | .983 | .933 | .982 |  | .978 | .930 | .982
$r=80$ |  | .995 | .973 | .995 |  | .987 | .966 | .994
$\tau=0.5$: |  |  |  |  |  |  |  |
$r=10$ |  | .497 | .655 | .551 |  | .584 | .652 | .555
$r=30$ |  | .732 | .824 | .742 |  | .771 | .820 | .733
$r=80$ |  | .872 | .914 | .872 |  | .892 | .913 | .874
$\tau=1.0$: |  |  |  |  |  |  |  |
$r=10$ |  | .732 | .830 | .765 |  | .751 | .790 | .727
$r=30$ |  | .855 | .908 | .859 |  | .868 | .900 | .856
$r=80$ |  | .932 | .956 | .933 |  | .945 | .959 | .938
Table 1: Model choice accuracies for the Pearson Bayes factor (PBF), JZS Bayes
factor (JZS), and BIC Bayes factor (BIC), calculated as the proportion of
simulated datasets for which the correct model was chosen.

Finally, I assessed the degree to which the model chosen by the Pearson Bayes
factor matches the model chosen by the JZS Bayes factor. This model choice
consistency is displayed in Table 2. Overall, model choice consistency is very
good (all above 84%). Certainly, consistency improves as the number of
replicates per treatment condition increases. However, it is also notable that
consistency improves across the board when $\alpha$ is set to 0.

|  | $\tau=0$ | $\tau=0.5$ | $\tau=1.0$
---|---|---|---|---
$\alpha=-1/2$ |  |  |  |
$r=10$ |  | .933 | .842 | .902
$r=30$ |  | .950 | .908 | .947
$r=80$ |  | .978 | .958 | .976
$\alpha=0$ |  |  |  |
$r=10$ |  | .955 | .932 | .961
$r=30$ |  | .952 | .951 | .968
$r=80$ |  | .979 | .979 | .986
Table 2: Model choice consistency for the PBF Bayes factor and the JZS Bayes
factor, calculated as the proportion of simulated datasets for which both
methods chose the same model.

##  6 Conclusion

In this paper, I have developed an analytic Bayes factor (the Pearson Bayes
factor) which allows a researcher to obtain Bayes factors using the minimal
summary statistics (e.g., the test statistic and the degrees of freedom) from
both $t$-tests and analysis of variance designs. This formula improves upon
the BIC Bayes factor formula of Faulkenberry, (2018) by not only providing the
user with an exact Bayes factor instead of an approximation, but by also
giving the user the ability to specify prior knowledge about the size of
effects expected. Our simulation study shows that the Pearson Bayes factor
performs comparably to the JZS Bayes factor of Rouder et al., (2012)
(especially when the prior parameter $\alpha$ is set to 0). Also, the Pearson
Bayes factor outperforms JZS on data generated under a null model. The Pearson
Bayes factor did not perform as well as the JZS Bayes on data generated under
the assumption of nonzero treatment effects, but this limitation is outweighed
by the Pearson Bayes factor’s unique ability to be computed from summary
statistics alone and without the need for computing a multi-dimensional
integral. For this reason, the Pearson Bayes factor is an excellent tool for
researchers who wish to assess the evidential value of data from published
studies in which raw data is not readily available. In all, the Pearson Bayes
factor should be a valuable tool in the scientist’s statistical toolbox.

## References

  * Faulkenberry, (2018) Faulkenberry, T. J. (2018).  Computing Bayes factors to measure evidence from experiments: An extension of the BIC approximation.  Biometrical Letters, 55(1):31–43.
  * (2) Faulkenberry, T. J. (2019a).  Estimating Bayes factors from minimal ANOVA summaries for repeated-measures designs.  arXiv:1905.05569.
  * (3) Faulkenberry, T. J. (2019b).  Estimating evidential value from analysis of variance summaries: A comment on Ly (2018).  Advances in Methods and Practices in Psychological Science, 2(4):406–409.
  * Fisher, (1925) Fisher, R. A. (1925).  Statistical Methods for Research Workers.  Oliver & Boyd, Edinburgh.
  * García-Donato and Sun, (2007) García-Donato, G. and Sun, D. (2007).  Objective priors for hypothesis testing in one-way random effects models.  Canadian Journal of Statistics, 35(2):303–320.
  * Gönen et al., (2005) Gönen, M., Johnson, W. O., Lu, Y., and Westfall, P. H. (2005).  The Bayesian two-sample $t$ test.  The American Statistician, 59(3):252–257.
  * Jeffreys, (1961) Jeffreys, H. (1961).  The Theory of Probability (3rd ed.).  Oxford University Press, Oxford, UK.
  * Kass and Raftery, (1995) Kass, R. E. and Raftery, A. E. (1995).  Bayes factors.  Journal of the American Statistical Association, 90(430):773.
  * Liang et al., (2008) Liang, F., Paulo, R., Molina, G., Clyde, M. A., and Berger, J. O. (2008).  Mixtures of $g$ priors for Bayesian variable selection.  Journal of the American Statistical Association, 103(481):410–423.
  * Maruyama, (2009) Maruyama, Y. (2009).  A Bayes factor with reasonable model selection consistency for ANOVA model.  arXiv:0906.4329v2.
  * Masson, (2011) Masson, M. E. J. (2011).  A tutorial on a practical Bayesian alternative to null-hypothesis significance testing.  Behavior Research Methods, 43(3):679–690.
  * Raftery, (1995) Raftery, A. E. (1995).  Bayesian model selection in social research.  Sociological Methodology, 25:111–163.
  * Rouder et al., (2016) Rouder, J. N., Engelhardt, C. R., McCabe, S., and Morey, R. D. (2016).  Model comparison in ANOVA.  Psychonomic Bulletin & Review, 23(6):1779–1786.
  * Rouder et al., (2012) Rouder, J. N., Morey, R. D., Speckman, P. L., and Province, J. M. (2012).  Default Bayes factors for ANOVA designs.  Journal of Mathematical Psychology, 56(5):356–374.
  * Schwarz, (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * Sellke et al., (2001) Sellke, T., Bayarri, M. J., and Berger, J. O. (2001).  Calibration of $p$-values for testing precise null hypotheses.  The American Statistician, 55(1):62–71.
  * Student, (1908) Student (1908).  The probable error of a mean.  Biometrika, 6(1):1.
  * Thiele et al., (2017) Thiele, J. E., Haaf, J. M., and Rouder, J. N. (2017).  Is there variation across individuals in processing? Bayesian analysis for systems factorial technology.  Journal of Mathematical Psychology, 81:40–54.
  * Wagenmakers, (2007) Wagenmakers, E.-J. (2007).  A practical solution to the pervasive problems of $p$ values.  Psychonomic Bulletin & Review, 14(5):779–804.
  * Wang and Liu, (2016) Wang, M. and Liu, G. (2016).  A simple two-sample Bayesian $t$-test for hypothesis testing.  The American Statistician, 70(2):195–201.
  * Wang and Sun, (2014) Wang, M. and Sun, X. (2014).  Bayes factor consistency for one-way random effects model.  Communications in Statistics - Theory and Methods, 43(23):5072–5090.
  * Wasserstein and Lazar, (2016) Wasserstein, R. L. and Lazar, N. A. (2016).  The ASA statement on $p$-values: Context, process, and purpose.  The American Statistician, 70(2):129–133.
  * Zellner, (1986) Zellner, A. (1986).  On assessing prior distributions and Bayesian regression analysis with $g$-prior distributions.  In Goel, P. K. and Zellner, A., editors, Bayesian inference and decision techniques: Essays in Honor of Bruno de Finetti, pages 233–243. Elsevier.


# Meta-learning for RIS-assisted NOMA Networks

Yixuan Zou1,Yuanwei Liu1,Kaifeng Han2,Xiao Liu1,Kok Keong Chai1
1 Queen Mary University of London, London, UK
2 China Academy of Information and Communications Technology, Beijing, China
{yixuan.zou, yuanwei.liu, x.liu<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

A novel reconfigurable intelligent surfaces (RISs)-based transmission
framework is proposed for downlink non-orthogonal multiple access (NOMA)
networks. We propose a quality-of-service (QoS)-based clustering scheme to
improve the resource efficiency and formulate a sum rate maximization problem
by jointly optimizing the phase shift of the RIS and the power allocation at
the base station (BS). A model-agnostic meta-learning (MAML)-based learning
algorithm is proposed to solve the joint optimization problem with a fast
convergence rate and low model complexity. Extensive simulation results
demonstrate that the proposed QoS-based NOMA network achieves significantly
higher transmission throughput compared to the conventional orthogonal
multiple access (OMA) network. It can also be observed that substantial
throughput gain can be achieved by integrating RISs in NOMA and OMA networks.
Moreover, simulation results of the proposed QoS-based clustering method
demonstrate observable throughput gain against the conventional channel
condition-based schemes.

##  I Introduction

Multiple-input-multiple-output (MIMO) has been recognized as a promising
technology to enhance the capacity and the spectral efficiency of fifth-
generation (5G) wireless networks [1]. However, deploying and configuring a
large number of antennas can lead to severe hardware impairment, heavy
computational cost, and substantial power consumption. To overcome these
limitations, reconfigurable intelligent surfaces (RISs) have emerged as a
cost-effective solution [2]. An RIS composes of a large number of low-cost
reflecting elements that can proactively reconfigure the propagation of
incident signals. By intelligently adjusting the phase shift of each RIS
element, the communication channels can be effectively manipulated to enhance
the spectral efficiency and the network coverage [3, 4]. Benefiting from the
low-cost meta-materials, RISs can be seamlessly integrated with emerging
technologies, such as MIMO, to further improve transmission throughput in a
cost-effective manner.

As a promising technology for supporting massive connectivity, non-orthogonal
multiple access (NOMA) [5] can be integrated into RIS-aided networks to
further enhance the spectral efficiency. Moreover, the integration of NOMA and
RIS can effectively improve the design flexibility of NOMA schemes [6].
However, the implementation of NOMA in RIS-aided networks increases the
resource allocation difficulty. In particular, RISs can alter the channel
quality of individual devices, which directly influences the decoding orders
and the clustering results of NOMA systems.

Extensive research contributions have been devoted to investigating the
integration of RIS and NOMA techniques [7, 8, 9]. In particular, the authors
of [7] derived extensive analytical results, including ergodic rates, energy
efficiency and spectral efficiency, for RIS-assisted NOMA networks. The sum
rate maximization problem of RIS-aided NOMA systems was investigated in [8],
where the passive beamforming at the RIS was jointly optimized with the active
beamforming at the base station (BS) under both the ideal and non-ideal RIS
elements. To maximize sum rate while ensuring user fairness, the authors of
[9] formulated a max-min problem for RIS-enhanced NOMA networks by jointly
optimizing the power allocation and the RIS phase shift. The aforementioned
research contributions revealed the potentials of RIS when integrated into
NOMA networks and established a foundation for solving various challenges in
RIS-assisted networks. However, these contributions mainly investigated the
implementations of conventional convex optimization techniques, which often
suffer from high computational complexity and poor scalability. Moreover, the
objective functions are often non-convex, which can not be directly tackled by
convex optimization methods. Hence, alternative cost-effective non-convex
optimization schemes are necessary to fulfil the requirements of massive
connectivity in next-generation wireless networks.

In recent years, artificial intelligence (AI) has emerged as a tremendous
technology to address the problems of exploding data volume, non-convex
optimization, and computational complexity [10]. In particular, deep learning
(DL) techniques utilize the extensive offline training phase to reduce the
algorithm complexity in the application and have received overwhelming
research interests in the optimization of RIS-assisted wireless communication
systems [11, 12, 13]. In [11], deep transfer learning was employed to solve
the beamforming optimization problem in multiple-input-single-output (MISO)
networks, based on a small amount of training data. The problem was further
extended into the discrete phase shift cases to accommodate hardware
limitations. The authors of [12] and [13] utilized neural networks to learn
the interactions between the receiver locations and the optimal RIS phase
shift to achieve maximal communication throughput. These aforementioned
contributions demonstrated the outstanding performance of DL-based techniques
when solving high-dimensional and non-convex optimization problems in RIS-
enhanced wireless networks.

However, existing DL-based RIS optimization methods [11, 12, 13] are all based
on orthogonal multiple access (OMA) systems, where the RIS phase shift is the
only optimization variable. To the best of our knowledge, there does not exist
a DL-based solution for RIS-aided NOMA networks, which motivates this study.
In this paper, we investigate the sum rate optimization problem in RIS-aided
downlink MISO-NOMA networks, where both the RIS phase shift and the BS power
allocation are optimized to maximize the total transmission sum rate. We adopt
the zero-forcing (ZF) precoding method and the successive interference
cancellation (SIC) decoding method to eliminate the effect of multi-user
interference on the strong users. However, this approach causes the weak users
to suffer from both inter-cluster and intra-cluster interference, leading to
poor achievable rates. To improve the resource efficiency, we propose a
quality-of-service (QoS)-based NOMA clustering method, which aims to maximize
the QoS deviation within each cluster. In terms of the DL model, we design a
neural network to output optimized power allocation given the RIS phase shift,
resulting in a low-complexity model. Meanwhile, the phase shift is optimized
through a gradient descent algorithm, given the trained network. We further
employ meta-learning in the training process to improve the convergence rate
of the phase shift optimization. The main contributions are as follows:

  1. 1.

We propose a RIS-enhanced NOMA downlink framework and formulate the sum rate
maximization problem by jointly optimizing the phase shift of the RIS and the
power allocation of the BS. To improve the resource efficiency, we propose a
QoS-based NOMA clustering scheme, which maximizes the QoS deviation within
clusters.

  2. 2.

We propose a MAML-based DL algorithm to solve the joint optimization problem.
The algorithm can output optimized solutions in as few as five iterations and
the model has lower network complexity compared to the conventional design.

  3. 3.

Simulation results indicate that the implementation of RIS can induce
approximately 5% to 25% throughput gain as the number of RIS elements
increases from 8 to 64, in both NOMA and OMA networks. Results also show that
the proposed QoS-based clustering method achieves higher throughput than the
conventional channel condition-based approach.

##  II System Model and Problem Formulation

###  II-A System Model

Figure 1: Illustration of the RIS-assisted downlink MISO-NOMA network.

As illustrated in Fig. 1, we consider a downlink MISO system with one BS and
$K$ mobile users (MUs). The BS is equipped with $M$ antenna elements and each
MU is equipped with one single antenna. The communication between the BS and
the MUs is assisted by a RIS of $N$ reflecting elements, whose phase shift and
amplitude absorption can be adjusted by a controller.

The channels between the BS and the RIS are modelled as Rician fading
channels. The BS-MU channels and the RIS-MU channels are modelled as Rayleigh
fading channels. The path loss of a particular MU $k$ is modelled as
$\text{PL}_{k}=d_{k}^{-\alpha}$ where $d_{k}$ is the distance, calculated in
meters, between the MU and the BS, and $\alpha$ denotes the path loss
exponent.

###  II-B NOMA Signal Model

In this subsection, we formulate the NOMA-based signal model and introduce the
proposed QoS-based clustering method.

####  II-B1 Signal model

The signal received at each MU is a composition of the signals derived from
the direct link between BS and MU, and the signals derived from the reflecting
link. In particular, for MU $i$ in the $l$-th cluster, we denote the RIS-MU
link and BS-MU link by $\mathbf{h}_{R,l,i}^{H}\in{\mathbb{C}^{1\times N}}$,
and $\mathbf{h}_{B,l,i}^{H}\in{\mathbb{C}^{1\times M}}$, respectively. We
further denote the BS-RIS link by ${\mathbf{H}}_{BR}\in{\mathbb{C}^{N\times
M}}$. The phase shift of the RIS is denoted by
$\bm{\theta}=[\theta_{1},\cdots,\theta_{n},\cdots,\theta_{N}]$ where
$\theta_{n}\in[0,2\pi)$. The diagonal phase-shifting matrix is expressed as
$\mathbf{\Theta}=\text{diag}(\beta_{1}e^{j\theta_{1}},\cdots,\beta_{n}e^{j\theta_{n}},\cdots,\beta_{N}e^{j\theta_{N}})$,
where $\beta_{n}\in[0,1]$ is the amplitude reflection coefficient. For
simplicity, we assume that all amplitude coefficients are ones, i.e.,
$\beta_{n}=1,\forall n$.

We assume that each cluster is formed by two users and we denote the strong MU
as MU $s$ and the weak MU as MU $w$. The signals transmitted to the strong and
the weak MU in the $l$-th cluster are denoted by ${s_{l,s}}$ and ${s_{l,w}}$,
respectively and we denote the transmit power allocated to the strong and the
weak MUs in the $l$-th cluster by $p_{l,s}$ and $p_{l,w}$, respectively.
Hence, the transmit signal of the $l$-th cluster is formulated as
${x_{l}}=\sqrt{{p_{l,s}}}{s_{l,s}}+\sqrt{{p_{l,w}}}{s_{l,w}}$. The
corresponding signal received by MU $i$ in the $l$-th cluster can be expressed
as

$\displaystyle{{y_{l,i}}=\left({{\mathbf{h}}_{B,l,i}^{H}+{\mathbf{h}}_{R,l,i}^{H}{{\mathbf{\Theta}}}{{\mathbf{H}}_{BR}}}\right)\sum\limits_{l=1}^{K/2}{{{\mathbf{w}}_{l}}{x_{l}}}+{n_{l,i}}},$
(1)

where ${{{\mathbf{w}}_{l}}}$ denotes the beamforming vector of the $l$-th
cluster and $n_{l,i}$ denotes the additive white Gaussian noise (AWGN),
modelled as $n_{l,i}\sim\mathcal{C}\mathcal{N}\left({0,\sigma^{2}}\right)$.

To decode the received symbols from the multiplexed signal $x_{l}$, each
strong MU employs the SIC technique which eliminates the intra-cluster
interference. The weak MU, however, decodes the signal directly without SIC.
Additionally, we eliminate the inter-cluster interference for the strong MUs
through the ZF beamforming technique. The normalized ZF precoding vector is
given by

$\displaystyle{\mathbf{w}_{l}}=\frac{{\mathbf{h}_{l,s}}{\left({{{\mathbf{h}^{H}_{l,s}}}{\mathbf{h}_{l,s}}}\right)^{-1}}}{\rho_{l,s}},$
(2)

where
$\mathbf{h}_{l,s}^{H}=\mathbf{h}^{H}_{B,l,s}+\mathbf{h}^{H}_{R,l,s}\mathbf{\Theta}\mathbf{H}_{BR}$
denotes the combined channel and $\rho_{l,s}$ denotes the normalizing constant
formulated as
$\rho_{l,s}=|{{\mathbf{h}_{l,s}}{({{{\mathbf{h}_{l,s}}^{H}}{\mathbf{h}_{l,s}}})^{-1}}}|^{2}$.
The corresponding ZF precoding constraints are expressed as follows

$\displaystyle\left\\{\begin{array}[]{*{20}{c}}\mathbf{h}_{j,s}^{H}\mathbf{w}_{l}&=0,&{\kern
1.0pt}\forall j\neq l,{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}j=1,\cdots,{K/2},\\\
\mathbf{h}_{j,s}^{H}\mathbf{w}_{l}&=\frac{1}{\rho_{l,s}},&j=l.\end{array}\right.$
(5)

Hence, without the interference, the signal received at the strong MU in the
$l$-th cluster can be simplified into

$\displaystyle{{y_{l,s}}{\kern
1.0pt}=\mathbf{h}^{H}_{l,s}{{\mathbf{w}}_{l}}\sqrt{{p_{l,s}}}{s_{l,s}}+{n_{l,s}}}.$
(6)

Based on (5), we can calculate the received SINR of the strong MU in the
$l$-th cluster as

$\displaystyle{{\gamma_{l,s}}=\frac{{{{\left|{\mathbf{h}^{H}_{l,s}{{\mathbf{w}}_{l}}\sqrt{{p_{l,s}}}{s_{l,s}}}\right|}^{2}}}}{{\sigma^{2}}}=\frac{p_{l,s}}{{\rho_{l,s}\sigma^{2}}}}.$
(7)

Since both inter-cluster interference and intra-cluster interference exists in
the weak MU’s received signals, the received SINR of the weak MU in the $l$-th
cluster is derived as

$\displaystyle{{\gamma_{l,w}}=\frac{{{{\left|{{{\mathbf{h}}_{l,w}}{{\mathbf{w}}_{l}}}\right|}^{2}}{p_{l,w}}}}{{{{\left|{{{\mathbf{h}}_{l,w}}{{\mathbf{w}}_{l}}}\right|}^{2}}{p_{l,s}}+{{\left|{{{\mathbf{h}}_{l,w}}\sum\limits_{j=1,j\neq
l}^{K/2}{{{\mathbf{w}}_{j}}{x_{j}}}}\right|}^{2}}+\sigma^{2}}}}.$ (8)

####  II-B2 QoS-based clustering scheme

When both the ZF precoding and the SIC decoding techniques are employed in
NOMA, the weak MUs suffer from both the inter-cluster and the intra-cluster
interference, resulting in low SINR and low achievable rate compared to the
strong MUs, who are served in an interference-free manner. Conventional
clustering methods allocate MUs by exploiting the difference between their
channel conditions. However, when MUs have different QoS requirements, we may
notice a weak MU acquiring a high QoS, which is challenging to achieve given
the multiuser interference. Moreover, due to the low SINR, a great amount of
transmit power has to be allocated to the weak MU to fulfil the QoS. Hence, it
is more sensible to assign MUs with lower QoS requirements as the weak MUs to
improve the resource efficiency and to enhance the network throughput.
Therefore, we propose a QoS-based clustering scheme, which assigns the MUs
with higher or lower QoS requirements as the strong or weak MUs, respectively.

To be specific, the objective of the QoS-based clustering method is to
maximize the minimum QoS deviation among all clusters. The clustering problem
can be formulated as
$\mathop{\max}\mathop{\min}\limits_{l={1,\cdots,K/2}}(R^{l,s}_{QoS}-R^{l,w}_{QoS})$,
where $R^{l,s}_{QoS}$ and $R^{l,w}_{QoS}$ denote the QoS requirements of the
strong and the weak MUs in the $l$-th cluster, respectively.

To achieve the maximal QoS deviation, we propose a simple clustering method
when each cluster consists of two MUs. We assume that $K$ is an even number
and all MUs are ordered in terms of their QoS requirements, namely, the $k$-th
MU has the $k$-th highest QoS requirement. The maximal QoS deviation can be
achieved by assigning the $k$-th MU and the $(k+K/2)$-th MU into the same
cluster, for all $k\leq K/2$.

###  II-C Problem Formulation

The optimization goal is to maximize the total throughput of the network by
jointly optimizing the RIS phase shift
${{\bm{\theta}}}=[{\theta_{1}},\cdots,{\theta_{n}},\cdots,{\theta_{N}}]$ and
the BS power allocation vector $\mathbf{P}=[\mathbf{P}_{s},\mathbf{P}_{w}]$,
where $\mathbf{P}_{s}=[p_{1,s},\cdots,p_{K/2,s}]$ and
$\mathbf{P}_{w}=[p_{1,w},\cdots,p_{K/2,w}]$. The optimization problem is
formulated as

$\displaystyle\mathop{\max}\limits_{\bm{\theta},\mathbf{P}}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}R=\sum\nolimits_{l=1}^{K/2}(R_{l,s}+R_{l,w})$ (9a)
$\displaystyle{\text{s}}{\text{.t}}{\text{.}}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{R_{l,i}}\geq R_{\text{QoS}}^{l,i},\forall l,\forall
i\in\\{s,w\\}$ (9b) $\displaystyle\left|{{e^{j\theta_{n}}}}\right|=1,\forall
n$ (9c) $\displaystyle\sum\nolimits_{l=1}^{K/2}(p_{l,s}+p_{l,w})\leq P_{max}$
(9d)

where $R_{l,i}=B_{l}\log_{2}(1+\gamma_{l,i})$ denotes the throughput achieved
by MU $i$ in cluster $l$ and $R_{\text{QoS}}^{l,i}$ denotes the minimal QoS
requirement of the given MU. Hence, (9b) represents the minimum transmit rate
constraint. Moreover, (9c) denotes the phase shift constraint of the RIS and
(8d) qualifies the total transmit power constraint of the BS. Due to the non-
convex constraint (9c), the optimization problem can not be directly solved by
conventional approaches. Hence, we proposed to tackle the joint optimization
problem utilizing machine learning techniques.

##  III DL-Based Power Allocation and Phase Shift Optimization

In this section, we introduce the proposed meta-learning enabled DL algorithm
that jointly optimizes the power allocation and the RIS phase shift, given
MUs’ QoS requirements.

###  III-A Proposed MAML-based Phase Shift and Power Allocation Optimization
Algorithm

The main idea of DL algorithms is to extensively train a neural network such
that, given any inputs, the outputs of the network achieve minimal loss. A
conventional design is to construct a neural network that outputs all
optimization variables, namely, the RIS phase shift and the power allocation.
However, this design will result in an extremely large input space that
consists of all channel information and QoS information. In particular, all
channel matrices contribute $(2K\times N+2K\times M+2N\times M)$ to the input
dimension, leading to exceedingly expensive computational costs. Moreover, the
phase shift and the power allocation have vastly different value ranges and
distributions, which greatly increase the training difficulty.

######  Remark 1.

Optimizing the phase shift requires the knowledge of all channels among the
BS, the RIS and the MUs. However, the optimization of the power allocation
only requires the information of the combined channel.

Inspired by the fact that the combined channel of all MUs, denoted by
$\mathbf{H}=[\mathbf{h}_{1,s},\mathbf{h}_{1,w},\cdots,\mathbf{h}_{K/2,s},\mathbf{h}_{K/2,w}]$,
provides sufficient channel information for optimizing $\mathbf{P}$ but not
$\bm{\theta}$, we propose to design the neural network to input the combined
channel and output the optimized power allocation. The real and the imaginary
parts of the combined channel $\mathbf{H}$ contributes $(2K\times M)$ to the
input dimension, which is significantly smaller than the input dimension of
the intuitive design. The neural network $G_{\eta}$ is formulated as

$\mathbf{P}=G_{\eta}(\mathbf{H}(\bm{\theta}),\mathbf{R}_{QoS},\mathbf{L}_{\text{path}}),$
(10)

where $\mathbf{H}(\bm{\theta})$ denotes the combined channel calculated using
the phase shift $\bm{\theta}$, $\mathbf{R}_{QoS}\in\mathbb{R}^{K}$ denotes the
QoS requirement vector, and $\mathbf{L}_{\text{path}}\in\mathbb{R}^{K}$ is the
path loss vector. The phase shift $\bm{\theta}$ is optimized separately using
a gradient descent algorithm. Two optimization algorithms are connected in an
alternating structure, by using the output of the other algorithm as the
input.

In contrast to the conventional alternating optimization approach, we train
the neural network $G_{\eta}$ to output the optimized power allocation given
any phase shift. Hence, given a trained $G_{\eta}$, we can find the optimized
pair of $\bm{\theta}$ and $\mathbf{P}$ by solely performing the optimization
on $\bm{\theta}$. To further improve the convergence rate of the gradient
descent algorithm, the network $G_{\eta}$ is trained using MAML, such that,
the optimized pair of $\bm{\theta}$ and $\mathbf{P}$ can be obtained in as few
as five iterations.

MAML is a meta-learning technique, which is designed to optimize the model
parameters such that a few gradient steps will produce a maximally effective
performance on a new task [14]. As demonstrated in [15], MAML can be employed
to reduce the number of gradient descent steps required to optimize the
network input space. In our model, the network inputs are optimized by
adjusting the phase shift $\bm{\theta}$. Moreover, the gradient descent steps
on $\bm{\theta}$ are performed by back-propagating through the weights of
$G_{\eta}$. Hence, we propose to train $G_{\eta}$ with MAML to obtain a set of
network weights that can greatly reduce the number of gradient steps required
to update $\bm{\theta}$.

###  III-B Loss Functions

As described in (9), both $\bm{\theta}$ and $\mathbf{P}$ need to be optimized
to maximize the system throughput given the constraints. Hence, they share the
same loss function, denoted by $\mathcal{L}(\bm{\theta},\bm{\eta})$, where
$\bm{\eta}$ is the weights of the neural network $G_{\eta}$ that outputs
$\mathbf{P}$. The loss function consists of two parts, the total throughput
and the constraint term for enforcing the QoS requirements, given by

$\mathcal{L}(\bm{\theta},\bm{\eta})=w_{1}\sum_{l=1}^{K/2}\sum_{i=s,w}R_{l,i}(\bm{\theta},\bm{\eta})+\\\
w_{2}\sum_{l=1}^{K/2}\sum_{i=s,w}\text{max}\big{(}R_{l,i}(\bm{\theta},\bm{\eta})-R^{l,i}_{QoS},0\big{)},$
(11)

where $R_{l,i}(\bm{\theta},\bm{\eta})$ is the sum rate calculated using
$\bm{\theta}$ and $\bm{\eta}$, and
$\text{max}\big{(}R_{l,i}(\bm{\theta},\bm{\eta})-R^{l,i}_{QoS},0\big{)}$
indicates the QoS deficiency of MU $i$ in cluster $l$. The weights $w_{1}$ and
$w_{2}$ are tuned during training. Since the sum rate is positive-valued and
the QoS deficiency is negative-valued, $w_{1}$ and $w_{2}$ should be negative
and positive, respectively.

Suppose we aim to optimize $\bm{\theta}$ in $J$ gradient steps, we can derive
the gradient descent formula in the $j$-th gradient step of the $p$-th
training episode as

$\bm{\theta}^{(j)}\leftarrow\bm{\theta}^{(j-1)}-\gamma_{\theta}\frac{\partial}{\partial\bm{\theta}^{(j-1)}}\mathcal{L}(\bm{\theta}^{(j-1)},\bm{\eta}^{(p-1)}),$
(12)

where $\gamma_{\theta}$ denotes the step size and $\bm{\eta}^{(p-1)}$ denotes
the neural network weights obtained in the previous training episode. In order
to satisfy the phase shift constraint in (9c), we clip the values of
$\bm{\theta}$ to $[0,2\pi]$ after each update. Based on (12), the loss
function after completing the $J$-th gradient step is therefore
$\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta}^{(p-1)})$. Then, we optimize the
neural network $G_{\eta}$ to further minimize
$\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta}^{(p-1)})$ through the following
update formula

$\bm{\eta}^{(p)}\leftarrow\bm{\eta}^{(p-1)}-\gamma_{\eta}\frac{\partial}{\partial\bm{\eta}^{(p-1)}}\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta}^{(p-1)}),$
(13)

where $\gamma_{\eta}$ denotes the learning rate.

Algorithm 1 Meta-learning Based Training Algorithm

0: Channel matrix $\mathbf{H}$, QoS vector $\mathbf{R}_{QoS}$, MU locations,
neural network $G_{\eta}$, number of phase shift update steps $J$, phase shift
learning rate $\gamma$

0: Trained neural network $G_{\hat{\eta}}$ Initialize $\bm{\eta}$

1: repeat

2: for each episode do

3: Initialize phase shift
$\theta_{1},...,\theta_{N}\overset{\text{iid}}{\sim}\mathcal{U}(0,2\pi)$

4: Calculate path loss vector $\mathbf{L}_{\text{path}}$

5: for $j=0$ to $J-1$ do

6: Obtain power allocation
$\mathbf{P}^{(j)}=G_{\eta}(\mathbf{H}(\bm{\theta}^{(j)}),\mathbf{R}_{QoS},\mathbf{L}_{\text{path}})$

7: Calculate loss function $\mathcal{L}(\bm{\theta}^{(j)};\bm{\eta})$

8: Update phase shift using (12)

9: end for

10: Given the optimized phase shift $\bm{\theta}^{(J)}$, calculate the
optimized power allocation
$\mathbf{P}^{(J)}=G_{\eta}(\mathbf{H}(\bm{\theta}^{(J)}),\mathbf{R}_{QoS},\mathbf{L}_{\text{path}})$

11: Calculate loss function $\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta})$ using
the optimized phase shift

12: Update network weights using (13)

13: end for

14: until reaches the maximum training steps

15: Return $G_{\hat{\eta}}$

To employ MAML, we calculate (13) by implicitly performing the second order
differentiation with respect to the loss function
$\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta}^{(p-1)})$ and back-propagating
through the $J$ phase shift optimization steps in (12). What’s more, we also
perform MAML on the hyper-parameter $\gamma_{\theta}$ against the loss
function $\mathcal{L}(\bm{\theta}^{(J)},\bm{\eta}^{(p-1)})$ to reduce the need
for further hyper-parameter tuning. The update equation of $\gamma_{\theta}$
can be derived in the same way as in (13).

Figure 2: An illustration of the MAML-based training framework.

###  III-C Training Algorithm

As shown in Fig. 2, each training epoch can be divided into two stages,
corresponding to the inner and outer MAML steps:

  1. 1.

Phase shift optimization (inner step): The initial phase shift is sampled
according to a random uniform distribution, i.e.
$\bm{\theta}^{(0)}\sim\mathcal{U}(0,2\pi)$. In the $j$-th gradient loop, the
corresponding power allocation $\mathbf{P}^{(j)}$ is obtained based on (10),
using $\bm{\theta}^{(j)}$. Then, $\bm{\theta}^{(j)}$ is optimized with respect
to the loss function $\mathcal{L}(\bm{\theta}^{(j)},\bm{\eta})$, as in (12).
We repeat (12) for $J$ iterations. The final optimized phase shift is thus
$\bm{\theta}^{(J)}$.

  2. 2.

Power allocation optimization (outer step): After completing $J$ gradient
descent loops, the current optimal power allocation $\mathbf{P}^{(J)}$ can be
computed using $\bm{\theta}^{(J)}$ and (10). Then, the network weights
$\bm{\eta}$ are updated according to (13), by backpropagating through all $J$
gradient descent iterations.

The pseudocode of the training algorithm is presented in Algorithm 1, where
lines 2-8 correspond to the phase shift optimization steps and lines 9-11
correspond to the power allocation optimization steps. To apply the trained
network on test datasets, we only need to perform the phase shift optimization
procedure for $J$ times, after which the optimized phase shift
$\bm{\theta}^{(J)}$ and the corresponding power allocation $\mathbf{P}^{(J)}$,
are the solutions to our joint optimization problem in (9).

###  III-D Complexity Analysis

The computational complexity of the proposed joint optimization algorithm
mainly depends on three factors, namely, the number of phase shift update
steps $J$, the complexity of the loss function, and the complexity of the
neural network. Trivially, the complexity of (11) is dominated by the
calculation of individual MU’s combined channel vector $\mathbf{h}^{H}_{l,i}$,
of which the complexity is $\mathcal{O}(NM)$. Since we compute the combined
channel for each MU, the total complexity induced by calculating the loss
function is $\mathcal{O}(NKM)$. Then, we note that a fully-connected neural
network of $D$ layers, including input and output layers, has a computational
complexity of $\mathcal{O}(\sum_{i=1}^{D-1}n_{i}n_{i+1})$, where $n_{i}$ is
the number of neurons in layer $i$. Therefore, the proposed algorithm has a
complexity of $\mathcal{O}(JNKM\sum_{i=1}^{D-1}n_{i}n_{i+1})$.

##  IV Simulations

In this section, we present the simulation results for the proposed MAML-based
optimization algorithm for RIS-assisted NOMA downlink networks. We consider
MUs moving in a square area of width 10 meters. The BS is located at a corner
of the area and the RIS is randomly placed in the area. In all simulations,
the path loss exponent is set to $\alpha=3$, the noise power spectral
frequency is -169 dBm/Hz, and the total bandwidth is 4 MHz. The number of
phase shift optimization steps is $J=5$ unless otherwise stated.

Figure 3: Sum rate versus the number of MUs, $K$, in NOMA and OMA systems
using the QoS-based or the channel condition-based clustering schemes. Figure
4: Sum rate versus the number of reflecting elements $N$ for NOMA and OMA
cases, given $M=16$ antennas and 20 dBm transmit power at BS. Figure 5: Sum
rate versus total transmit power at BS for NOMA and OMA cases, given $M=16$ BS
antennas and $N=16$ reflecting elements.

###  IV-A QoS-based versus Channel Condition-based Clustering Methods

In Fig. 3, we compare the network throughput given different clustering
methods, namely, the proposed QoS-based clustering scheme and the conventional
channel condition-based clustering scheme. Simulations are performed with
$M=24$ BS antennas and $P_{\text{max}}=20$ dBm maximum transmit power. It can
be observed that, for smaller numbers of MUs, i.e., $K\leq 8$, the performance
difference between two clustering schemes is small because their clustering
results are likely to be similar. However, when there are more MUs in the
system, i.e., $K\geq 12$, the proposed QoS-based method starts to achieve
higher throughput than the conventional channel condition-based approach. The
performance gains further increases as the number of MUs increases.

###  IV-B Sum Rate versus the Number of RIS Elements

In Fig. 4, we observe that the proposed NOMA network with the QoS-based
clustering scheme outperforms the conventional OMA network of around 4 dBm/Hz
of sum rate, without the enhancement of RIS. In both NOMA and OMA networks,
the deployment of RIS induces approximately 5% to 25% throughput gain as the
number of reflecting elements ranges from $N=8$ to $N=64$. Higher performance
gain can be attained by increasing the number of RIS elements, however, the
cost of optimization complexity and the deployment cost increases as well.

###  IV-C Sum Rate versus BS Total Transmission Power

Fig. 5 illustrates the throughput performance between OMA and NOMA systems as
the BS power varies between 0 dBm and 30 dBm. It can be observed that the NOMA
system outperforms the OMA system for all values of BS transmit power, given
the same number of MUs. We also notice that the NOMA network with 4 MUs
achieves higher throughput than the OMA network with 6 MUs when the BS power
is less than 15 dBm. Moreover, as the number of MUs increases, the throughput
of the NOMA networks increases by a larger amount compared to the throughput
of the OMA networks.

##  V Conclusions

In this article, we proposed a QoS-based clustering method to improve the
resource efficiency in RIS-assisted NOMA networks. We formulated the sum rate
maximization problem by jointly optimizing the RIS phase shift and the BS
power allocation. The proposed DL solution utilized a low-complexity network
architecture and employed MAML to improve the convergence rate. Simulation
results demonstrated that higher transmission throughput was achieved by the
proposed QoS-based clustering method than the baseline method. Results also
illustrated that the proposed QoS-based NOMA model achieved higher throughput
compared to the conventional OMA models for a wide range of values of the BS
transmit power and the RIS elements. Moreover, the implementation of RIS
improved network throughput by a substantial amount, which further grows as
the number of reflecting elements increases.

## References

  * [1] L. Sanguinetti, E. Björnson, and J. Hoydis, “Toward massive MIMO 2.0: Understanding spatial correlation, interference suppression, and pilot contamination,” _IEEE Trans. Commun._ , vol. 68, no. 1, pp. 232–257, Jan. 2020.
  * [2] Y. Liu, X. Liu, X. Mu, T. Hou, J. Xu, M. Di Renzo, and N. Al-Dhahir, “Reconfigurable intelligent surfaces: Principles and opportunities,” _IEEE Commun. Surv. & Tut._, May 2021, accepted to appear, doi:10.1109/COMST.2021.3077737.
  * [3] X. Liu, Y. Liu, and Y. Chen, “Machine learning empowered trajectory and passive beamforming design in UAV-RIS wireless networks,” _IEEE J. Sel. Areas Commun._ , Dec. 2020, doi:10.1109/JSAC.2020.3041401.
  * [4] X. Liu, Y. Liu, Y. Chen, and H. V. Poor, “RIS enhanced massive non-orthogonal multiple access networks: Deployment and passive beamforming design,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 4, pp. 1057–1071, Apr. 2021.
  * [5] Y. Liu, Z. Qin, M. Elkashlan, Z. Ding, A. Nallanathan, and L. Hanzo, “Nonorthogonal multiple access for 5G and beyond,” _Proc. IEEE_ , vol. 105, no. 12, pp. 2347–2381, Dec. 2017.
  * [6] Y. Liu, X. Mu, X. Liu, M. Di Renzo, Z. Ding, and R. Schober, “Reconfigurable intelligent surface (RIS) aided multi-user networks: Interplay between NOMA and RIS,” _arXiv:2011.13336_ , Nov. 2020.
  * [7] T. Hou, Y. Liu, Z. Song, X. Sun, Y. Chen, and L. Hanzo, “Reconfigurable intelligent surface aided NOMA networks,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2575–2588, Jul. 2020.
  * [8] X. Mu, Y. Liu, L. Guo, J. Lin, and N. Al-Dhahir, “Exploiting intelligent reflecting surfaces in NOMA networks: Joint beamforming optimization,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6884–6898, Jul. 2020\.
  * [9] G. Yang, X. Xu, and Y.-C. Liang, “Intelligent reflecting surface assisted non-orthogonal multiple access,” in _IEEE Wireless Commun. Netw. Conf. (WCNC 2020)_ , May 2020, pp. 1–6.
  * [10] S. Sun, Z. Cao, H. Zhu, and J. Zhao, “A survey of optimization methods from a machine learning perspective,” _IEEE Trans. Cybern._ , vol. 50, no. 8, pp. 3668–3681, Nov. 2020.
  * [11] Y. Ge and J. Fan, “Beamforming optimization for intelligent reflecting surface assisted MISO: A deep transfer learning approach,” _IEEE Trans. Veh. Technol._ , pp. 1–1, Mar. 2021.
  * [12] C. Huang, G. C. Alexandropoulos, C. Yuen, and M. Debbah, “Indoor signal focusing with deep learning designed reconfigurable intelligent surfaces,” in _IEEE Int. Workshop Signal Process. Adv. Wireless Commun. (SPAWC)_ , Jul. 2019, pp. 1–5.
  * [13] B. Sheen, J. Yang, X. Feng, and M. M. U. Chowdhury, “A deep learning based modeling of reconfigurable intelligent surface assisted wireless communications for phase shift configuration,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 262–272, Jan. 2021.
  * [14] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” _Proc. Int. Conf. Mach. Learn. (ICML 2017)_ , pp. 1126–1135, Jul. 2017.
  * [15] Y. Wu, M. Rosca, and T. Lillicrap, “Deep compressed sensing,” _Proc. Int. Conf. Mach. Learn._ , vol. 97, pp. 6850–6860, May 2019.