# Variational Bayesian hierarchical regression for data analysis

Dennis Becker<EMAIL_ADDRESS>Department of Information Systems, Leuphana
University

###### Abstract

Collected data, which is used for analysis or prediction tasks, often have a
hierarchical structure, for example, data from various people performing the
same task. Modeling the data’s structure can improve the reliability of the
derived results and prediction performance of newly unobserved data. Bayesian
modeling provides a tool-kit for designing hierarchical models. However,
Markov Chain Monte Carlo methods which are commonly used for parameter
estimation are computationally expensive. This often renders its use for many
applications not applicable. However, variational Bayesian methods allow to
derive an approximation with much less computational effort. This document
describes the derivation of a variational approximation for a hierarchical
linear Bayesian regression and demonstrates its application to data analysis.

Paper Draft - arxiv.org

##  1 Introduction

Bayesian methods allow to develop models that describe the data generation
process, derive confidence bounds on its estimated parameters and predictions
of new observations. The Bayesian inference to derive the model parameters,
however, especially for hierarchical models, quickly becomes intractable. But
the analyzed data often have a hierarchical structure. For example, in social
science data can be observed from different people that perform the same task.
Therefore, it appears natural to model such a structure to derive an estimate
for each subject, provide individual predictions of future observations, and
conclude about the general population.

For the estimation of such models usually, Markov Chain Monte Carlo (MCMC)
methods are used, which approximate the model by sampling from a Markov chain
with the stationary distribution of the posterior distribution [6, 4].
Although these methods provide guarantees about the samples that are taken
from the targeted density [11], they are computationally expensive even for
small data sets.

Thus, variational inference can be suited for larger data sets and scenarios
where the model or a variety of models have to be estimated more quickly. It,
however, only provides an approximation of the posterior distribution, does
not provide the same guarantees such as MCMC methods, and underestimates the
variance of the posterior distribution [7].

In the following, we will describe the derivation of a variational Bayesian
hierarchical regression model. The model is very closely related to the
Bayesian linear regression model that has been described by Drugowitsch (2013)
and Bishop (2006), to which we add a hierarchical prior distribution. After
deriving the model, we will demonstrate the use on a freely available data set
[5]. An implementation of the algorithm in R is
available:https://github.com/dennisthemenace2/hBReg/

##  2 Variational Inference

Probabilistic models enable to describe how the observed variables
$\mathbf{X}=x_{i},...,x_{n}$ have been generated by the influence of a number
of latent variables $\mathbf{Z}$. Under these latent variables, we summarize
all parameters that are used to model the relations in the observed data. The
derived probabilistic model is specified by the joint density
$p(\mathbf{X},\mathbf{Z})$. To derive an estimate of the latent variables
given the observed data we aim to compute $p(\mathbf{Z}|\mathbf{X})$ which is
called the posterior distribution. Using Bayes’ rule the posterior
distribution is stated as:

$p(\mathbf{Z}|\mathbf{X})=\frac{p(\mathbf{X},\mathbf{Z})}{p(\mathbf{X})}=\frac{p(\mathbf{X}|\mathbf{Z})p(\mathbf{Z})}{p(\mathbf{X})}=\frac{p(\mathbf{X}|\mathbf{Z})p(\mathbf{Z})}{\int
p(\mathbf{X},\mathbf{Z})d\mathbf{Z}}.$

The denominator is called the marginal likelihood or model evidence.
Typically, the calculation of the evidence is unavailable in closed-form,
therefore inference in such models is often based on Markov Chain Monte-Carlo
methods [6, 4]. Alternatives are approximative methods such as variational
inference [1, 7].

Variational inference aims to approximate the posterior distribution $p$ with
a distribution $q$ from a set of tractable distribution $Q$. It, therefore,
states the inference problems as an optimization problem, which allows to
estimate an approximate solution $q(\mathbf{Z})\approx
p(\mathbf{Z}|\mathbf{X})$.

To solve this optimization problem, we need a measure of the similarity
between $q$ and $p$. The most common type of variational Bayes to describe the
difference between these two distributions is the Kullback-Leibler ($KL$)[8]
divergence which measures the differences in the information contained within
two distributions:

$KL(q||p)=\sum_{x}q(x)\log\frac{q(x)}{p(x)}.$ (1)

The difference between both distributions is always greater or equal to zero
$KL(q||p)\geq 0$ for all $q$,$p$, and only equal to zero if $q=p$. Next, we
plug the posterior distribution $p(\mathbf{Z}|\mathbf{X})$ into the $KL$
divergence:

$\displaystyle KL(q(\mathbf{Z})||p(\mathbf{Z}|\mathbf{X}))$
$\displaystyle=\sum_{\mathbf{Z}}q(\mathbf{Z})\log\frac{q(\mathbf{Z})}{p(\mathbf{Z}|\mathbf{X})}$
$\displaystyle=\sum_{\mathbf{Z}}q(\mathbf{Z})\log\frac{q(\mathbf{Z})}{p(\mathbf{X},\mathbf{Z})}+\log
p(\mathbf{X})$ $\displaystyle=KL(q(\mathbf{Z})||p(\mathbf{X},\mathbf{Z}))+\log
p(\mathbf{X})$

First, we pulled the normalizing constant, or marginal probability
$p(\mathbf{X})$ out, and recognize the remaining term as the $KL$ divergence
of the unnormalized posterior distribution. By maximizing the unnormalized
$KL$ divergence, we minimize the above defined $KL$ divergence. With respect
to the variational distribution $q(\mathbf{Z})$, the log model evidence $\log
p(\mathbf{X})$ is constant. Since
$KL(q(\mathbf{Z})||p(\mathbf{Z}|\mathbf{X}))\geq 0$, we get by rearranging
terms that

$\log
p(\mathbf{X})=KL(q(\mathbf{Z})||p(\mathbf{Z}|\mathbf{X}))-KL(q(\mathbf{Z})||p(\mathbf{X},\mathbf{Z})).$

We notice that the log model evidence $\log p(\mathbf{X})$ is equal to the
difference between the normalized and unnormalized $KL$ divergence.
Furthermore, the $KL$ divergence of the unnormalized posterior distribution is
a lower bound on the model evidence. Due to this property, the negative
$KL(q(\mathbf{Z})||p(\mathbf{X},\mathbf{Z}))$ is called the variational lower
bound or the evidence lower bound (ELBO). Minimizing the lower bound amounts
to maximizing a lower bound on the model evidence.

To complete the specification of the optimization problem, we need to describe
the variational family $Q$. A widely used class of distributions is the mean-
field approximation, which assumes an independent factorization:

$q(\mathbf{Z})=\prod_{m=1}^{M}q_{m}(Z_{m}).$

Each latent variables $Z_{m}$ can be governed by its own variational factor
$q_{m}(Z_{m})$, which renders the individual factors mutually independent. For
the mean-field choice of $Q$, we can optimize the problem using coordinate
descent. We iterate over the variational factors $q_{m}(Z_{m})$ and for each
$m$ we optimize the evidence lower bound over $q_{m}$ while keeping the other
variational factors $q_{j}(Z_{j})$ constant. This results into the log of the
optimal solution for each factor:

$\ln q_{m}(Z_{m})=\mathbb{E}_{j\neq m}\Big{[}\ln
p(\mathbf{X},\mathbf{Z})\Big{]}+\text{const}.$ (2)

This procedure iteratively fits the fully-factored
$q(\mathbf{Z})=q_{1}(Z_{1})q_{2}(Z_{2})\cdots q_{m}(Z_{m})$ approximation of
$p(\mathbf{Z}|\mathbf{X})$.

###  2.1 Hierarchical Variational Bayesian Regression

An often encountered scenario is that the analyzed data has a hierarchical
structure, for example, people are undergoing the same treatment type. Then we
can consider the data from each client as an independent sample from the same
population, which naturally suggest a hierarchical structure of the data.
Typically, the observed data consist of some target variable
$\mathbf{y}_{i}\in\mathbb{R}$ and independent variables
$\mathbf{x}_{i}\in\mathbb{R}^{D}$ for each client. If we assume a linear
relationship between those, the suggested model is a hierarchical linear
regression.

We further assume that we want the benefits of a Bayesian model and that it
has to be estimated quickly, therefore we want to derive a hierarchical
Bayesian linear regression model. The derived model is very similar to the
Bayesian linear regression model derived by Drugowitsch (2013) and Bishop
(2006). In addition to these models, we have the hierarchical prior for the
individual people or clients’ weights. The complete model is shown in plate
notation in Figure 1.

Figure 1: Plate notation of the hierarchical linear regression model

The model is shown with all components, which makes it appear quite confusing,
but this allows to talk about the design of the probabilistic graphical model
and about the reasoning behind it. The filled nodes are observed or known and
circles typically represent distributions, whereas boxes represent constants.
The rectangles with names in it are plates, which indicate that their elements
are observed mutable times. We start from the bottom of the illustration with
the observed data $\mathbf{Y}$, regarding the nature of the target variable we
assume a continues variable for the regression model. This suggests a normal
distribution $\mathcal{N}(y\mid\mu,\sigma)$, which has two parameters, the
mean value $\mu$ and a variance parameter $\sigma$. These parameters have both
to be specified by either a distribution or a constant value. That is why
these models naturally appear to grow in one direction. Next, we discuss the
plates surrounding the data. These illustrate that we have multiple clients or
subjects from which we observed repeated measures or observations. We specify
the inverse of the variance parameter using the conjugate prior a gamma
distribution $\operatorname{Gamma}(\sigma|a_{0},b_{0})$. The Gamma
distribution has two parameters that we need to specify. For these values, we
assume both small values, which represents an uninformative prior. The
specification of this prior distribution on the precision leads to an
estimation of the noise from the data. We continue with specifying the $\mu$
parameter of the normal distribution which is the target variable. Based on
the linear relationship between the observed data and the target variable we
use matrix multiplication of the observed data $\mathbf{x}_{ci}$ and a client-
specific weight vector $\beta_{c}$ which is also represented by the conjugate
prior of the mean value a normal distribution. This completes the description
of the observed data and we sate the factors:

$\displaystyle p(\mathbf{Y}\mid\mathbf{x}\beta_{i},\sigma)$
$\displaystyle=\prod_{i=1}^{C}\prod_{c=1}^{M}\mathcal{N}(y_{c,m}\mid\mathbf{x}_{c,m}\beta_{i},\sigma^{-1}),$
$\displaystyle p(\sigma)$ $\displaystyle=\operatorname{Gamma}(\sigma\mid
a_{0},b_{0}).$

In an iterative fashion, we continue describing the client individual
parameters $\beta_{i}$. The prior on the mean value of the client individual
parameters which is titled $\Delta$ represents the population of which all
individuals are samples. The precision $s$ describes the influence of
population or hierarchical prior on the individual. We describe the clients
specific factors as:

$\displaystyle p(\beta_{i}\mid\Delta,s)$
$\displaystyle=\prod_{i=1}^{C}\mathcal{N}\left(\beta_{i}\mid\Delta,s^{-1}\right)$
$\displaystyle p(s)$ $\displaystyle=\operatorname{Gamma}(s\mid c_{0},d_{0})$

For the prior on the mean of the hierarchical prior $\Delta$, we use the
constant 0, which encourages the weights to become small. It shrinks the
weights towards 0 similar to a quadratic regularisation term in ridge
regression [2]. For the prior on the precision, we chose a Gamma distribution
for each observed variable (or Dimension) individually. The reasoning is that
during optimization of the model irrelevant coefficients will shrink
automatically. This process is known as automatic relevance determination
(ARD)[9, 13, 14]. This concludes the model construction with the final
factors:

$\displaystyle p(\Delta\mid 0,\mathbf{w})$
$\displaystyle=\mathcal{N}(\Delta\mid 0,\ \mathbf{w}^{-1})$ $\displaystyle
p(\mathbf{w})$
$\displaystyle=\prod_{d=1}^{D}\operatorname{Gamma}(\mathbf{w}_{d}\mid
e_{0},f_{0})$

In order to continue with the model development and variational approximation,
we state the models joint probability:

$\displaystyle p(\mathbf{Y},\mathbf{x},\beta_{i},\sigma,\Delta,s,\mathbf{w})$
$\displaystyle=p(\mathbf{Y}\mid\mathbf{x}\beta_{i},\sigma)p(\beta_{i}\mid\Delta,s)p(\Delta\mid
0,w)p(\sigma)p(s)p(\mathbf{w}),$

and assume that the posterior distribution is approximated by the factored
variational posterior distributions,

$\displaystyle p(\beta_{i},\Delta,\sigma,s,w\mid\mathbf{D})\approx
q(\beta_{i})q(\Delta)q(\sigma)q(s)q(\mathbf{w}).$

Now, we have to derive the update equations for the variational factors using
the coordinate descent algorithm by applying Equation 2.

###  2.2 Update factor q($\beta_{i}$)

To derive the update equations, we select the factors that depend on beta,
since the others are constant and absorbed into a constant term. We replace
the factors with their actual distributions and multiply everything out to
separate the terms that not depend on $\beta,\Delta$ or $\sigma$ out into the
constant term. The terms are rearranged and grouped until we receive a
distribution that we recognize. Since we used the conjugate prior and the
product of two normal distributions, the result will also be a normal
distribution.

$\displaystyle\ln q_{\beta_{i}}^{*}(\beta_{i})$
$\displaystyle=\operatorname{E}_{\Delta,\sigma,s}\left[\ln
p(\mathbf{y}\mid\mathbf{x}\beta_{i},\sigma)+\ln
p(\beta_{i}\mid\Delta,s)\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\sigma}\left[\sum_{c=1}^{M}\frac{1}{2}\ln\sigma-\frac{\sigma}{2}(y_{c}-\mathbf{x}_{c}\beta_{i})^{2}\right]+\operatorname{E}_{\Delta,s}\left[\frac{1}{2}\ln
s-\frac{s}{2}(\Delta-\beta_{i})^{2}\right]+C$
$\displaystyle=\operatorname{E}_{\sigma}\left[\sum_{c=1}^{M}-\frac{\sigma}{2}(y_{c}-x_{c}\beta_{i})^{2}\right]+\operatorname{E}_{\Delta,s}\left[-\frac{s}{2}(\Delta-\beta_{i})^{2}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\sigma}\left[-\frac{\sigma}{2}\sum_{c=1}^{M}{y_{c}^{2}}+\sigma\beta_{i}\sum_{c=1}^{M}{y_{c}\mathbf{x}_{c}}-\frac{\sigma\beta_{i}^{2}\sum_{c=1}^{M}\mathbf{x}^{2}_{c}}{2}\right]+\operatorname{E}_{\Delta,s}\left[-\frac{s\Delta^{2}}{2}+s\Delta\beta_{i}-\frac{s\beta_{i}^{2}}{2}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\sigma}\left[\sigma\beta_{i}\sum_{c=1}^{M}{y_{c}\mathbf{x}_{c}}-\frac{\sigma\beta_{i}^{2}\sum_{c=1}^{M}\mathbf{x}^{2}_{c}}{2}\right]+\operatorname{E}_{\Delta,s}\left[s\Delta\beta_{i}-\frac{s\beta_{i}^{2}}{2}\right]+\text{const}$
$\displaystyle=\underbrace{(\operatorname{E}_{\sigma}\left[\sigma\right]\sum_{c=1}^{M}{\mathbf{x}_{c}y_{c}}+\operatorname{E}_{s}\left[s\right]\operatorname{E}_{\Delta}\left[\Delta\right])}_{\text{Mean}}\beta_{i}-\frac{1}{2}\underbrace{(\operatorname{E}_{\sigma}\left[\sigma\right]\sum_{c=1}^{M}\mathbf{x}^{2}_{c}+\operatorname{E}_{s}\left[s\right])}_{\text{Variance
}}\beta_{i}^{2}+\text{const}$

By completing the square over $\beta$, we can derive the parameters of the
Gaussian distribution. There are still the expectations of the variables with
respect to their variational distribution, which we consider to be constant,
that we have to replace. The variational distribution of the terms
$\operatorname{E}_{\sigma}\left[\sigma\right]$ and
$\operatorname{E}_{s}\left[s\right]$ is a Gamma distribution and the expected
value for a Gamma distribution is
$\operatorname{E}\left[\operatorname{Gamma}(a,b)\right]=\frac{a}{b}$.
Therefore, the expect values are
$\operatorname{E}_{\sigma}\left[\sigma\right]=\frac{a_{n}}{b_{n}}$ and
$\operatorname{E}_{s}\left[s\right]=\frac{c_{n}}{d_{n}}$ respectively. The
expectation for $\Delta$ is with respect to a normal distribution, the
expected value of a normal distribution is the mean, which results into the
expected values of $\operatorname{E}_{\Delta}\left[\Delta\right]=\Delta$. This
leads to the following update equations for $\beta_{i}$:

$\displaystyle\lambda_{\beta_{i}}$
$\displaystyle=\frac{a_{n}}{b_{n}}\sum_{c=1}^{M}\mathbf{x}^{2}_{c}+\frac{c_{n}}{d_{n}},$
$\displaystyle\beta_{i}$
$\displaystyle=\frac{\frac{a_{n}}{b_{n}}\sum_{c}^{M}{\mathbf{x}_{c}y_{c}}+\frac{c_{n}}{d_{n}}\Delta}{\lambda_{\beta_{i}}}.$

The update equations shows that the client individual weights are dependent on
their individual data and the hierarchical prior.

###  2.3 Update factor q($\Delta$)

Next, we estimate the update equations for the factor q($\Delta$) to derive
the variational posterior distribution. The procedure is the same as before,
we select the terms that depend on $\Delta$ because the other terms get
absorbed into the constant term. We multiply everything out and rearrange
until we derive the same form as for the previous factor, where we recognize
the terms of the resulting normal distribution.

$\displaystyle\ln q_{\Delta}^{*}(\Delta)$
$\displaystyle=\operatorname{E}_{\beta_{i},s}\left[\ln
p(\beta_{i}\mid\Delta,s)+\ln
p(\Delta\mid\mu_{0},\mathbf{w})\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},s}\left[\sum_{i=1}^{C}\frac{1}{2}\ln
s-\frac{s}{2}(\Delta-\beta_{i})^{2}\right]+\operatorname{E}_{\mathbf{w}}\left[\frac{1}{2}\ln\mathbf{w}-\frac{\mathbf{w}}{2}(\Delta-\mu_{0})^{2}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},s}\left[\sum_{i=1}^{C}-\frac{s}{2}(\Delta-\beta_{i})^{2}\right]+\operatorname{E}_{\mathbf{w}}\left[-\frac{\mathbf{w}}{2}(\Delta-\mu_{0})^{2}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},s}\left[-\frac{Cs\Delta^{2}}{2}+s\Delta\sum_{i=1}^{C}{\beta_{i}}\right]+\operatorname{E}_{\mathbf{w}}\left[\mathbf{w}\Delta-\frac{w\Delta^{2}}{2}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},s}\left[-\frac{1}{2}(Cs+w)\Delta^{2}\right]+\operatorname{E}_{\mathbf{w},s,\beta_{i}}\left[\Delta(s\sum_{i}^{C}{\beta_{i}}+\mathbf{w}\mu_{0})\right]+\text{const}$

Similar to before, we have to fill in the expected values
$\operatorname{E}_{\mathbf{w}}\left[\mathbf{w}\right]=\frac{e_{n}}{\mathbf{f}_{n}}$,
and $\operatorname{E}_{\beta_{i}}\left[\beta_{i}\right]=\beta_{i}$, and the
expected value of $s$ has already been stated in the update of the previous
factor. This leads to the following update equations for $\Delta$:

$\displaystyle\lambda_{\Delta}$
$\displaystyle=C\frac{c_{n}}{d_{n}}+\frac{e_{n}}{\mathbf{f}_{n}},$
$\displaystyle\Delta$
$\displaystyle=\frac{\frac{c_{n}}{d_{n}}\sum_{i}^{I}{\beta_{i}}+\frac{e_{n}}{\mathbf{f}_{n}}\mu_{0}}{\lambda_{\Delta}}.$

We indicated that $\mathbf{f_{n}}$ is a vector with dimensions entries whereas
$e_{n}$ is a single value. Because each weight will have its own Gamma prior,
but the parameter $e$ will be the same for all $\mathbf{f_{n}}$ assuming they
have the same prior value. We will see this in the derivation of the update
for $\mathbf{w}$. The update equations show that the hierarchical prior is
similar to a mean weight, where it consists of the sum of all weights which
are scaled by the individual Gamma distributions of $\mathbf{w}$.

###  2.4 Update factor q($\sigma$)

Now, we derive the update equation for the precision of the noise $\sigma$. We
proceed as in the previous updates and since we use a conjugate prior, we
expect a Gamma distribution as the result.

$\displaystyle\ln q_{\sigma}^{*}(\sigma)$
$\displaystyle=\operatorname{E}_{\beta_{i}}\left[\ln
p(\mathbf{Y}\mid\mathbf{x}\beta_{i},\sigma)+\ln p(\sigma\mid
a,b)\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i}}\left[\sum_{i=1}^{C}\sum_{c=1}^{M}\frac{1}{2}\ln\sigma-\frac{\sigma}{2}(y_{i,c}-\mathbf{x}\beta_{i})^{2}\right]+\left[(a-1)\ln\sigma-b\sigma\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i}}\left[-\frac{\sigma}{2}\sum_{i=1}^{C}\sum_{c=1}^{M}{(y_{i,c}-\mathbf{x}\beta_{i})^{2}}\right]+(a-1)\ln\sigma-b\sigma+\frac{\sum_{i=1}^{C}\sum_{c=1}^{M}{1}}{2}\ln\sigma+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i}}\left[-\frac{\sigma}{2}\sum_{i=1}^{C}\sum_{c=1}^{M}{y_{i,c}^{2}-2y_{i,c}\mathbf{x}\beta_{i}+\mathbf{x}^{2}\beta_{i}^{2}}\right]+((a-1)+\frac{\sum_{i=1}^{C}\sum_{c=1}^{M}{1}}{2})\ln\sigma-b\sigma+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i}}\left[-\underbrace{((\sum_{i=1}^{C}\mathbf{x}\beta_{i}^{2}-2\sum_{i=1}^{C}\sum_{c=1}^{M}y_{i,c}\mathbf{x}^{2}\beta_{i}+\sum_{i=1}^{C}\sum_{c=1}^{M}{y_{i,c}^{2})\frac{1}{2}+b)}}_{\text{Rate
parameter
b}}\sigma\right]+\underbrace{((a-1)+\frac{\sum_{i=1}^{C}\sum_{c=1}^{M}{1}}{2})}_{\text{Shape
parameter a}}\ln\sigma+\text{const}$

We recognize the result as a $\operatorname{log}\operatorname{Gamma}$
distribution. This time we need to estimate the expected value of the
$\beta^{2}$, note that
$\operatorname{E}\left[\mathbf{X}^{2}\right]=\operatorname{E}\left[\mathbf{X}\right]+Var(\mathbf{X})$.
This leads to
$\operatorname{E}_{\beta_{i}}\left[\beta_{i}^{2}\right]=\beta_{i}+\lambda_{\beta_{i}}$
for the expectation of $\beta_{i}^{2}$. If we substitute the expectations and
rearrange the terms, we receive the following update equations:

$\displaystyle N$ $\displaystyle=\sum_{i=1}^{C}\sum_{c=1}^{M}{1},$
$\displaystyle a_{n}$ $\displaystyle=a+\frac{N}{2},$ $\displaystyle b_{n}$
$\displaystyle=b+\frac{1}{2}\left(\sum_{i=1}^{C}\sum_{c=1}^{M}{(y_{i,c}-\mathbf{x}_{i,c}\beta_{i})^{2}}+\sum_{i}^{C}\mathbf{x}^{T}_{i}\lambda_{\beta_{i}}\mathbf{x}_{i}\right).$

Considering the update equations, we notice that it represents the noise of
measuring the target with the prediction error term and the second term is the
sum of the standard errors. If we assume that the fitting error on the data is
high compared to the number of the sample, the expected value
$\frac{a_{n}}{b_{n}}$ of the precision would be low, so the variance is high.
Of course, the opposite would be true too if we make no prediction error but
have high uncertainty in the estimated clients’ weights, we expect the
measures to be less reliable.

###  2.5 Update factor q(s)

The result for factor $s$ which regulates the influence of the hierarchical
prior on the client individual weights will also be a Gamma distribution.

$\displaystyle\ln q_{s}^{*}(s)$
$\displaystyle=\operatorname{E}_{\beta_{i},\Delta}\left[\ln
p(\beta_{i}\mid\Delta,s)+\ln p(s\mid c,d)\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},\Delta}\left[\sum_{i=1}^{C}\sum_{d=1}^{D}\frac{1}{2}\ln
s-\frac{s}{2}(\beta_{i,d}-\Delta_{d})^{2}\right]+\left[(c-1)\ln
s-ds\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},\Delta}\left[\sum_{i=1}^{C}\sum_{d=1}^{D}-\frac{s}{2}(\beta_{i,d}-\Delta_{d})^{2}\right]+((c-1)+\frac{\sum_{i=1}^{C}\sum_{d=1}^{D}{1}}{2})\ln
s-ds+\text{const}$
$\displaystyle=\operatorname{E}_{\beta_{i},\Delta}\left[-((\sum_{i=1}^{C}\sum_{d=1}^{D}\beta_{i,d}^{2}-2\sum_{i=1}^{C}\sum_{d=1}^{D}\beta_{i,d}\Delta_{d}+\sum_{d=1}^{D}{\Delta_{d}^{2}})\frac{1}{2}+d)s\right]+((c-1)+\frac{CD}{2})\ln
s+\text{const}$

The result contains the squared expected values for $\beta$ and $\Delta$. The
expected value for $\Delta$ is derived similar to the one of $\beta$ in the
previous derivation which results in
$\operatorname{E}_{\Delta}\left[\Delta^{2}\right]=\Delta+\lambda_{\Delta}$.
After substituting and rearranging the terms we derive the update equations:

$\displaystyle c_{n}$ $\displaystyle=c+\frac{DC}{2},$ $\displaystyle d_{n}$
$\displaystyle=d+\frac{1}{2}\left(\sum_{i=1}^{C}\sum_{d=1}^{D}{(\beta_{i,d}-\Delta_{d})^{2}}+Spur(\lambda_{\Delta})+\sum_{i}^{C}Spur(\lambda_{\beta_{i}})\right).$

The update equations show that the variance depends on the squared differences
between the individual weights and the hierarchical prior, the variance of the
clients’ data, and the variance of the estimated weights. Similarly, if the
difference between the hierarchical prior and the clients’ weights is large,
the hierarchical prior will have less influence on the individual clients’
weights. However, if we imagine this term to be zero, then the variance of the
clients’ data and variance of the estimated weights still contributes. This
means, if there is a high variance in the clients’ data or variance between
the estimated client individual weights, the influence of that prior will be
reduced.

###  2.6 Update factor $q(\mathbf{w})$

Finally, we derive the update equation for $\mathbf{w}$. For simplicity, we
show the derivation for one of the factors.

$\displaystyle\ln q_{\mathbf{w}_{d}}^{*}(\mathbf{w}_{d})$
$\displaystyle=\operatorname{E}_{\Delta}\left[\ln p(\Delta_{d}\mid
0,\mathbf{w}_{d})+\ln p(\mathbf{w}_{d}\mid e,f)\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\Delta}\left[\frac{1}{2}\ln\mathbf{w}_{d}-\frac{w}{2}(\Delta_{d})^{2}\right]+\left[(e-1)\ln\mathbf{w}_{d}-f\mathbf{w}_{d}\right]+\text{const}$
$\displaystyle=\operatorname{E}_{\Delta}\left[-\frac{\mathbf{w}_{d}}{2}(\Delta_{d})^{2}\right]+((e-1)+\frac{1}{2})\ln\mathbf{w}_{d}-f\mathbf{w}_{d}+\text{const}$
$\displaystyle=\operatorname{E}_{\Delta}\left[-\mathbf{w}_{d}(\frac{1}{2}(-\Delta_{d})^{2}+f)\right]+((e-1)+\frac{1}{2})\ln\mathbf{w}_{d}+\text{const}$

The result shows that this is too a log Gamma distribution, where we have to
substitute the expected values by the momentums of the distribution. To keep
the notation simple, we show the update equation for $\mathbf{f}$ as a vector.

$\displaystyle e_{n}$ $\displaystyle=e+\frac{1}{2}$
$\displaystyle\mathbf{f}_{n}$
$\displaystyle=f+\frac{1}{2}(\Delta^{2}+trace(\lambda_{\Delta}))$

If we inspect the update equation for $\mathbf{w}$, we notice how it penalizes
the weights for each dimension. The first term is a quadratic penalty for
being away from zero and the second term adds a penalty which describes the
variance in the individual weights. If a particular weight varies among the
clients, this is stronger penalized than when it shows less variation among
the clients.

###  2.7 Variational lower bound

After we derived all update equations, we know that an iterative updating of
the variational factors will maximize the variational lower bound. In order to
determine when the algorithm is converged, we estimate the value of the
variational lower bound. This allows us to keep track of the optimization
process and provides an approximation of the evidence $\ln p(\mathbf{x})$.
Therefore, we plug our model joint distribution and variational distribution
into the negative $KL$ divergence from Equation 1. This time, however, we use
the continues definition and because the evidence lower bound was defined as
the negative divergence, the factors in the log fraction are switched. The
variational lower bound $\mathcal{L}(q)$ is then given by:

$\displaystyle\mathcal{L}(q)$ $\displaystyle=\int\int\int\int\int
q(\mathbf{\beta},\Delta,\sigma,s,\mathbf{w})\ln\left\\{\frac{p(\mathbf{\beta},\Delta,\sigma,s,\mathbf{w}\mid\mathbf{D})}{q(\beta,\Delta,\sigma,s,\mathbf{w})}\right\\}d\mathbf{\beta}\;d\Delta\;d\sigma\;ds\;d\mathbf{w}$
$\displaystyle=\;\mathbb{E}_{\mathbf{\beta},\Delta,\sigma,s,\mathbf{w}}\Big{[}\ln
p(\mathbf{\beta},\Delta,\sigma,s,\mathbf{w}\mid\mathbf{D})\Big{]}-\mathbb{E}_{\beta,\Delta,\sigma,s,\mathbf{w}}\Big{[}q(\beta,\Delta,\sigma,s,\mathbf{w})\Big{]}$
$\displaystyle=\;\mathbb{E}_{\mathbf{\beta},\sigma}\Big{[}\ln
p(\mathbf{Y}\mid\mathbf{x}\beta,\sigma)\Big{]}+\mathbb{E}_{\beta,\Delta,s}\Big{[}\ln
p(\beta\mid\Delta,s)\Big{]}+\mathbb{E}_{\mathbf{w}}\Big{[}\ln p(\Delta\mid
0,\mathbf{w})\Big{]}+\mathbb{E}_{\sigma}\Big{[}\ln p(\sigma)\Big{]}$
$\displaystyle\quad+\mathbb{E}_{s}\Big{[}\ln
p(s)\Big{]}+\mathbb{E}_{\mathbf{w}}\Big{[}\ln
p(\mathbf{w})\Big{]}-\mathbb{E}_{\beta}\Big{[}\ln
q(\beta)\Big{]}-\mathbb{E}_{\Delta}\Big{[}\ln
q(\Delta)\Big{]}-\mathbb{E}_{\sigma}\Big{[}\ln
q(\sigma)\Big{]}-\mathbb{E}_{\mathbf{w}}\Big{[}\ln q(\mathbf{w})\Big{]}$

The terms that involve expectations of the variational distributions
$\operatorname{log}q(\cdot)$ are the entropies $\mathbb{H}(\cdot)$ of that
distribution. The various terms are given in the following:

$\displaystyle\mathbb{E}_{\mathbf{\beta},\sigma}\Big{[}\ln
p(\mathbf{Y}\mid\mathbf{x}\beta,\sigma)\Big{]}$
$\displaystyle=\frac{N}{2}(\psi(an)-\ln
bn)-\frac{\sigma}{2}\left(\sum^{C}(\mathbf{x}_{i}\beta_{i}-\mathbf{y}_{i})^{T}(\mathbf{x}_{i}\beta_{i}-\mathbf{y}_{i})+\sum^{C}\mathbf{x}^{T}_{i}\lambda_{\beta_{i}}\mathbf{x}_{i}\right),$
$\displaystyle\mathbb{E}_{\mathbf{\beta},\Delta,\sigma}\Big{[}\ln
p(\beta\mid\Delta,s)\Big{]}$ $\displaystyle=\frac{C}{2}(\psi(c_{n})-\ln
d_{n})-\frac{s}{2}\left(\sum^{C}(\beta_{i}-\Delta)^{T}(\beta_{i}-\Delta)+Spur(\lambda_{\Delta})+\sum^{C}Spur(\lambda_{\beta_{i}})\right),$
$\displaystyle\mathbb{E}_{\Delta,\sigma}\Big{[}\ln p(\Delta\mid
0,\mathbf{w})\Big{]}$ $\displaystyle=\sum^{D}\frac{D}{2}\psi(e_{n})-\ln
f_{n}-\sum^{D}\frac{\mathbf{w}_{d}}{2}(\Delta^{2}_{d}+\lambda_{\Delta_{d,d}}),$
$\displaystyle\mathbb{E}_{\sigma}\Big{[}\ln p(\sigma)\Big{]}$
$\displaystyle=(a_{0}-1)(\psi(a_{n})-\ln bn)-b_{0}\sigma,$
$\displaystyle\mathbb{E}_{s}\Big{[}\ln p(s)\Big{]}$
$\displaystyle=(c_{0}-1)(\psi(c_{n})-\ln d_{n})-d_{0}s,$
$\displaystyle\mathbb{E}_{\mathbf{w}}\Big{[}\ln p(\mathbf{w})\Big{]}$
$\displaystyle=\sum^{D}(e_{0}-1)(\psi(e_{n})-\ln\mathbf{f}_{n_{d}})-e_{0}\mathbf{w}_{d},$
$\displaystyle\mathbb{E}_{\mathbf{\beta}}\Big{[}\ln q(\mathbf{\beta})\Big{]}$
$\displaystyle=\frac{1}{2}\sum^{C}\ln det(\lambda_{\beta_{i}}),$
$\displaystyle\mathbb{E}_{\Delta}\Big{[}\ln q(\Delta)\Big{]}$
$\displaystyle=\frac{1}{2}\ln det(\lambda_{\Delta}),$
$\displaystyle\mathbb{E}_{\sigma}\Big{[}\ln q(\sigma)\Big{]}$
$\displaystyle=a_{n}-\ln b_{n}+\ln\Gamma(a_{n})+(1-a_{n})\psi(a_{n}),$
$\displaystyle\mathbb{E}_{s}\Big{[}\ln q(s)\Big{]}$ $\displaystyle=c_{n}-\ln
d_{n}+\ln\Gamma(c_{n})+(1-c_{n})\psi(c_{n}),$
$\displaystyle\mathbb{E}_{\mathbf{w}}\Big{[}\ln q(\mathbf{w})\Big{]}$
$\displaystyle=\sum^{D}e_{n}-\ln\mathbf{f}_{n_{d}}+\ln\Gamma(e_{n})+(1-e_{n})\psi(e_{n}),$

where $\psi(\cdot)$ is the digamma function. During the optimization, the
bound is maximized and the optimization stopped when it reaches a plateau
$|\mathcal{L}(q_{n})-\mathcal{L}(q_{n+1})|<\epsilon$.

###  2.8 Predictive density

After we obtained the approximated posterior distribution, we might want to
make predictions for new target variables $y_{*}$ based on new observations
$\mathbf{x}_{*}$. Therefore, we require the posterior predictive distribution,
which allows us to estimate confidence intervals. The posterior predictive
distribution for the new target variable of the newly observed sample is
calculated by marginalizing the distribution of the target variable given the
observation and parameters over the posterior distribution of the parameters.
We substitute the posterior distribution for the variational posterior.

$\displaystyle p(y_{*}|\mathbf{x}_{*},\mathbf{D})$ $\displaystyle=\int\int
p(y_{*}\mid\mathbf{x}_{*},\beta,\sigma)p(\beta,\sigma\mid\mathbf{D})\;d\beta\;d\sigma$
$\displaystyle\approx\int\int
p(y_{*}\mid\mathbf{x}_{*},\beta,\sigma)q(\beta,\sigma)\;d\beta\;d\sigma$
$\displaystyle=\int\int\mathcal{N}\left(y_{*}\mid\mathbf{x}_{*}\beta,\sigma^{-1}\right)\mathcal{N}\left(\beta\mid\beta_{i},\lambda_{\beta_{i}}\right)\operatorname{Gamma}(\sigma\mid
a_{n},b_{n})\;d\beta\;d\sigma$
$\displaystyle=\int\mathcal{N}\left(y_{*}\mid\mathbf{x}_{*}\beta_{i},\sigma^{-1}+\mathbf{x}_{*}^{T}\lambda_{\beta}\mathbf{x}_{*}\right)\operatorname{Gamma}(\sigma\mid
a_{n},b_{n})\;d\sigma$
$\displaystyle=St\left(y_{*},\frac{a_{n}}{b_{n}}+\mathbf{x}_{*}^{T}\lambda_{\beta_{i}}\mathbf{x}_{*},2a_{n}\right)$

To obtain the result, we used results for the convolution of normal and Gamma
distributions [2, 10]. First, we convoluted the two normal distributions to
integrate out $\beta$. The marginal distribution of the resulting normal with
the Gamma distribution results into a Student’s t distribution with mean
$\mathbf{x}_{*}\beta_{i}$, precision
$\frac{a_{n}}{b_{n}}+\mathbf{x}_{*}^{T}\lambda_{\beta}\mathbf{x}_{*}$, and
$2a_{n}$ degrees of freedom. The result shows that the predicted uncertainty
is the sum of the noise $\sigma^{-1}$ and the variance in client individual
weights $\beta_{i}$. The number of degrees of freedom is approximately the
number of observed samples. After we have observed 30 samples we could
approximate this with a normal distribution. Interestingly, these samples do
not have to come from an individual since it is the number of the overall
observed samples. We can also make a prediction for clients that we have never
observed before by using the hierarchical prior $\Delta$ and the variance $s$.

##  3 Data Analysis

For a demonstration of the algorithm, we chose a freely available dataset from
the machine learning repository
(http://archive.ics.uci.edu/ml/datasets/Turkiye+Student+Evaluation/). We chose
the Turkiye Student Evaluation Data Set [5], which is a data set that consists
of the evaluation score of courses from the Gazi University in Ankara
(Turkey). To rate the individual courses, 28-course specific questions were
inquired and 5 additional attributes. An overview of the questions and items
is shown in Table 1

Item | Description
---|---
instr | Instructor’s identifier
class | Course code (descriptor)
repeat | Number of times the student is taking this course
attendance | Code of the level of attendance
difficulty | Level of difficulty of the course as perceived by the student; values taken from
Q1 | The semester course content, teaching method and evaluation system were provided at the start.
Q2 | The course aims and objectives were clearly stated at the beginning of the period.
Q3 | The course was worth the amount of credit assigned to it.
Q4 | The course was taught according to the syllabus announced on the first day of class.
Q5 | The class discussions, homework assignments, applications and studies were satisfactory.
Q6 | The textbook and other courses resources were sufficient and up to date.
Q7 | The course allowed field work, applications, laboratory, discussion and other studies.
Q8 | The quizzes, assignments, projects and exams contributed to helping the learning.
Q9 | I greatly enjoyed the class and was eager to actively participate during the lectures.
Q10 | My initial expectations about the course were met at the end of the period or year.
Q11 | The course was relevant and beneficial to my professional development.
Q12 | The course helped me look at life and the world with a new perspective.
Q13 | The Instructor’s knowledge was relevant and up to date.
Q14 | The Instructor came prepared for classes.
Q15 | The Instructor taught in accordance with the announced lesson plan.
Q16 | The Instructor was committed to the course and was understandable.
Q17 | The Instructor arrived on time for classes.
Q18 | The Instructor has a smooth and easy to follow delivery/speech.
Q19 | The Instructor made effective use of class hours.
Q20 | The Instructor explained the course and was eager to be helpful to students.
Q21 | The Instructor demonstrated a positive approach to students.
Q22 | The Instructor was open and respectful of the views of students about the course.
Q23 | The Instructor encouraged participation in the course.
Q24 | The Instructor gave relevant homework assignments/projects, and helped/guided students.
Q25 | The Instructor responded to questions about the course inside and outside of the course.
Q26 | The Instructor’s evaluation system (midterm, assignments, etc.) effectively measured the course objectives.
Q27 | The Instructor provided solutions to exams and discussed them with students.
Q28 | The Instructor treated all students in a right and objective manner.
Q1-Q28 are all Likert-type, meaning that the values are taken from 1,2,3,4,5
Table 1: Items and questions in the data set

The data has a hierarchical structure in a way that each class is measured
multiple times. However, we can also see that each instructor has multiple
classes, but for modeling this structure we would have to extend the model.
Therefore, we neglect this variable for now. The question of interest is, what
are the characteristics of a difficult class, and can we predict the
difficulty given the questions and items. Since the data in the data set is
ordinal, it is clearly not perfectly suited for analysis with the developed
algorithm, but it should help in demonstrating the algorithm. We start the
analysis with a linear regression model and the results are shown in Table 2.

Item | Estimate | Std. Error | t value |  Pr($>$$|$t$|$) | Item | Estimate | Std. Error | t value |  Pr($>$$|$t$|$)
---|---|---|---|---|---|---|---|---|---
nb.repeat | 0.858 | 0.024 | 35.32 | <0.000 | Q14 | 0.082 | 0.048 | 1.72 | 0.086
attendance | 0.453 | 0.011 | 39.58 | <0.000 | Q15 | 0.020 | 0.044 | 0.47 | 0.642
Q1 | 0.065 | 0.027 | 2.36 | 0.018 | Q16 | -0.151 | 0.040 | -3.74 | <0.001
Q2 | 0.035 | 0.035 | 0.99 | 0.320 | Q17 | 0.201 | 0.035 | 5.81 | <0.000
Q3 | -0.020 | 0.031 | -0.65 | 0.517 | Q18 | -0.074 | 0.039 | -1.89 | 0.059
Q4 | 0.026 | 0.033 | 0.80 | 0.424 | Q19 | -0.012 | 0.040 | -0.30 | 0.763
Q5 | 0.065 | 0.037 | 1.73 | 0.083 | Q20 | 0.001 | 0.042 | 0.02 | 0.987
Q6 | -0.019 | 0.034 | -0.56 | 0.575 | Q21 | 0.024 | 0.046 | 0.52 | 0.605
Q7 | 0.020 | 0.038 | 0.52 | 0.603 | Q22 | 0.048 | 0.046 | 1.03 | 0.302
Q8 | 0.064 | 0.036 | 1.80 | 0.073 | Q23 | -0.033 | 0.043 | -0.75 | 0.453
Q9 | -0.037 | 0.030 | -1.23 | 0.219 | Q24 | 0.053 | 0.039 | 1.34 | 0.181
Q10 | -0.091 | 0.041 | -2.21 | 0.027 | Q25 | 0.087 | 0.042 | 2.11 | 0.035
Q11 | -0.028 | 0.031 | -0.91 | 0.363 | Q26 | -0.084 | 0.036 | -2.32 | 0.020
Q12 | -0.026 | 0.030 | -0.85 | 0.396 | Q27 | 0.032 | 0.032 | 1.02 | 0.308
Q13 | 0.013 | 0.042 | 0.30 | 0.767 | Q28 | 0.003 | 0.036 | 0.09 | 0.927
Table 2: Result of the linear regression to explain the difficulty of a class

The regression results show that both items, the number of a repeated taking
of the class and attendance indicate a class with higher difficulty.
Apparently, students are more inclined to visit the lectures, when the content
is difficult. Regarding the questions, the 1st, 10th, 16th, 17th, 25th, and
26th show a significant relationship to the course difficulty. The first
question describes if the course topic and content were presented adequately
at the start of the course, potentially instructors of difficult classes tend
to state the goals and evaluation method of difficult classes more clearly in
the beginning. The 10th item is negatively related to the difficulty, which
states if the expectancies of the course result were met. Supposedly, students
are more pleased with the outcome of the class if its perceived difficulty is
lower. The 16th question is negatively associated with the difficulty, which
asks if the instructor was committed to the course and easily understandable.
Enthusiastic and clearly understandable instructors are associated with easier
overall classes. Interestingly, the 17th question inquires if the instructor
arrived in time and is positively associated with the class difficulty. This
suggests that if the instructor arrives early to the class, its content might
be more difficult. Similarly, the 25th question indicates that if the
instructor responded to questions inside and outside the class, the course was
more difficult. Question 26 indicates that the evaluations methods of courses
with higher difficulty were perceived as effective to measure the course
objective.

This provides some overview of the data. We could fit the variational model to
the data to estimate the weights and confidence intervals. However, keep in
mind that variational inference only provides an approximation and will
underestimate the variance in the posterior distribution. To visualize this,
we train the model on the data and use a Gibbs sampler for the estimation of
the same model. The weight estimates and 95% confidence intervals are
visualized in Figure 2

Figure 2: Comparison of the estimated weights and confidence interval of the
variational approximation and a Gibbs sampler

The variational model also provides a ranking of the importance of the
features, which we compare to the feature ranking of the Lasso regression[12].
For the estimation of the ranking for the Lasso regression, we add the class
variable and rank the features according to its first appearance in the lambda
path. The hierarchical model, however, does not have a variance parameter
$\mathbf{w}_{i}$ for each class. Therefore, we use the estimated weights for
their ranking. Because we used two different estimates for ranking of the
features for the variational model, we present the results in separately in
Table 3.

| Lasso regression | Variational model
---|---|---
1 | attendance | attendance
2 | nb.repeat | nb.repeat
3 | Q17 | Q17
4 | Q9 | Q16
5 | Q11 | Q26
6 | Q12 | Q9
7 | Q10 | Q22
8 | Q16 | Q1
9 | Q1 | Q25
10 | Q22 | Q18
11 | Q25 | Q5
12 | Q26 | Q2
13 | Q18 | Q10
14 | Q4 | Q11
15 | Q2 | Q3
16 | Q5 | Q14
17 | Q14 | Q4
18 | Q3 | Q21
19 | Q21 | Q12
20 | Q24 | Q27
21 | Q7 | Q7
22 | Q27 | Q20
23 | Q13 | Q24
24 | Q19 | Q28
25 | Q23 | Q13
26 | Q28 | Q6
27 | Q20 | Q19
28 | Q6 | Q23
29 | Q8 | Q15
30 | Q15 | Q8
Comparison of the ranking for the questions

| Lasso Regression | Variational model
---|---|---
1 | class 7 | class 7
2 | class 13 | class 8
3 | class 2 | class 2
4 | class 1 | class 11
5 | class 6 | class 4
6 | class 11 | class 6
7 | class 10 | class 3
8 | class 5 | class 12
9 | class 3 | class 9
10 | class 9 | class 5
11 | class 4 | class 10
12 | class 12 | class 1
13 | class 8 | class 13
Comparison of the ranking for the classes

Table 3: Comparison of estimated feature importance for Lasso regression and
the variational model

The importance ranking of the questions is similar in both models. However,
the ranking for the classes, where we use the absolute value of the weights
for the variational model appears to differ. Although both models agree on the
importance of class 7, the Lasso regression ranks class 13 as second most
important whereas the variational model ranks this last. Utilizing the
individual classes weights for ranking might not be a reliable method to rank
the importance.

As a next objective, we are interested in the prediction performance of the
variational model in comparison to similar models. Therefore, we use 10 fold
cross where we train the model on 9 folds and predict the remaining fold.
Before we estimate the mean square error of that fold, we round the
predictions. For the estimation of the cross-validated error, we use linear
regression, lasso regression, ridge regression, the variational hierarchical
Bayesian regression, and a linear mixed-effects model with an individual
intercept for each class.

| Regression | Lasso | Variational model | LME model | Ridge
---|---|---|---|---|---
Mean square error | 1.45 | 1.43 | 1.45 | 1.45 | 1.54
Table 4: Cross validated mean square error of the various models

The results of the cross-validated error in Table 4 show that the Lasso
regression has the lowest prediction error followed by the linear regression.
The error of the hierarchical regression is higher but still lower than the
prediction error of the ridge regression. Its prediction error is comparable
to the error of the linear mixed effect model. Apparently, the setup of the
analysis and data are not favoring the hierarchical models, however, we hope
that this provides some intuition into the models use, capabilities, and
limitations.

##  4 Concluding Remarks

We derived a variational approximation for a hierarchical Bayesian regression
model and demonstrated its application on real-world data set. The main
benefit of this model might be that it can be estimated considerably faster
than the same model using Gibbs sampling. Therefore, we can see its use for
application where the model has to be estimated quickly, repeatedly, or on a
bigger dataset. Furthermore, we hope that his document could shed some light
and raise interest in Bayesian modeling, variational approximation, and its
potential use in applications and data analysis.

##  5 Acknowledgements

We want to thank the machine learning repository and the publisher of the data
because they made this analysis possible.

## References

  * [1] MJ Matthew J Beal.  Variational algorithms for approximate Bayesian inference.  PhD Thesis, (May):1–281, 2003.
  * [2] Christopher M. Bishop.  Pattern Recognition and Machine Learning (Information Science and Statistics).  Springer-Verlag, Berlin, Heidelberg, 2006.
  * [3] Jan Drugowitsch.  Variational Bayesian inference for linear and logistic regression.  2013\.
  * [4] S Geman and D Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  IEEE Transactions on Pattern Analysis and Machine Intelligence, (6):721–741.
  * [5] N. Gunduz and E. Fokoue.  Uci machine learning repository, 2013.
  * [6] W. K. Hastings.  Monte carlo sampling methods using markov chains and their applications.  Biometrika, 57(1):97–109, 1970.
  * [7] Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul.  An introduction to variational methods for graphical models.  Mach. Learn., 37(2):183–233, November 1999.
  * [8] S. Kullback and R. A. Leibler.  On information and sufficiency.  Ann. Math. Statist., 22(1):79–86, 1951.
  * [9] MacKay.  Baysian Interpolation.  MIT Press Journal, 447:415–447, 1992.
  * [10] Kevin P. Murphy.  Machine learning : a probabilistic perspective.  MIT Press, 1 edition, August 2013.
  * [11] Christian P. Robert and George Casella.  Monte Carlo Statistical Methods (Springer Texts in Statistics).  Springer-Verlag, Berlin, Heidelberg, 2005.
  * [12] Robert Tibshirani.  Regression Selection and Shrinkage via the Lasso, 1996.
  * [13] Michael E. Tipping.  The Relevance Vector Machine.  In Advances in Neural Information Processing Systems (NIPS’ 2000), number 1, pages 652–658, 2000.
  * [14] David Wipf and Srikantan Nagarajan.  New view of automatic relevance determination.  In Advances in Neural Information Processing Systems 20, pages 1625–1632.


* [14] L. Onsager, Il Nuovo Cimento (1943-1954) 6, 279 (Mar 1949), 10.1007/BF02780991.
  * [15] A. A. Migdal, Turbulence as statistics of vortex cells, in The First Landau Institute Summer School, 1993: Selected Proceedings, ed. V. Mineev (Gordon and Breach, 1993), pp. 178–204.  arXiv:hep-th/9306152v2 [hep-th].
  * [16] G. L. Eyink, A. Gupta and T. A. Zaki, Journal of Fluid Mechanics 901, A2 (2020), 10.1017/jfm.2020.491.
  * [17] V. Yakhot and V. Zakharov, Physica D: Nonlinear Phenomena 64, 379 (1993), https://doi.org/10.1016/0167-2789(93)90050-B.
  * [18] V. Aji and N. Greenfeld, Physical Review Letters 86, 1007–1010 (Feb 2001), 10.1103/physrevlett.86.1007.
  * [19] A. Zamolodchikov and A. Zamolodchikov, Nuclear Physics B 477, 577–605 (Oct 1996), 10.1016/0550-3213(96)00351-3.
  * [20] I. Klebanov, String theory in two dimensions (1991), arXiv:hep-th/9108019v2.
  * [21] H. Dorn and H.-J. Otto, Nuclear Physics B 429, 375–388 (Oct 1994), 10.1016/0550-3213(94)00352-1.
  * [22] G. E. Volovik and K. Zhang, Phys. Rev. Research 2, 023263 (Jun 2020), 10.1103/PhysRevResearch.2.023263.
  * [23] K. Zhang, Phys. Rev. Research 2, 043356 (Dec 2020), 10.1103/PhysRevResearch.2.043356.
  * [24] K. P. Iyer, K. R. Sreenivasan and P. K. Yeung, Phys. Rev. X 9, 041006 (Oct 2019), 10.1103/PhysRevX.9.041006.
  * [25] K. P. Iyer, S. S. Bharadwaj and K. R. Sreenivasan, Area rule for circulation and minimal surfaces in three-dimensional turbulence (2020), arXiv:2007.06723 [physics.flu-dyn].
  * [26] T. Matsuzawa and W. Irvine, Realization of confined turbulence through multiple vortex ring collisions (03/12/2019), https://www.quantamagazine.org/an-unexpected-twist-lights-up-the-secrets-of-turbulence-20200903/, "Talk at the Flatiron Conference Universality Turbulence Across Vast Scales".


* Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378.
  * Golub and van Loan, (1996) Golub, G. H. and van Loan, C. F. (1996).  Matrix Computations.  John Hopkins University Press, 3 edition.
  * Hahn et al., (2018) Hahn, P. R., Martin, R., and Walker, S. G. (2018).  On recursive bayesian predictive distributions.  Journal of the American Statistical Association, 113(523):1085–1093.
  * Hall, (1987) Hall, P. (1987).  On Kullback-Leibler loss and density estimation.  The Annals of Statistics, 15(4):1491–1519.
  * Hjort and Jones, (1996) Hjort, N. L. and Jones, M. C. (1996).  Locally parametric nonparametric density estimation.  The Annals of Statistics, pages 1619–1647.
  * J. Ross and Markwick, (2019) J. Ross, G. and Markwick, D. (2019).  dirichletprocess: Build Dirichlet Process Objects for Bayesian Modelling.  R package version 0.3.1.
  * Jara et al., (2011) Jara, A., Hanson, T., Quintana, F., Müller, P., and Rosner, G. (2011).  DPpackage: Bayesian semi- and nonparametric modeling in R.  Journal of Statistical Software, 40(5):1–30.
  * Keith et al., (2010) Keith, M., Jameson, A., Van Straten, W., Bailes, M., Johnston, S., Kramer, M., Possenti, A., Bates, S., Bhat, N., Burgay, M., et al. (2010).  The High Time Resolution Universe Pulsar Survey I, System configuration and initial discoveries.  Monthly Notices of the Royal Astronomical Society, 409(2):619–627.
  * Lavine, (1992) Lavine, M. (1992).  Some aspects of Polya tree distributions for statistical modelling.  Annals of Statistics, 20(3):1222–1235.
  * Lavine, (1994) Lavine, M. (1994).  More aspects of Polya tree distributions for statistical modelling.  The Annals of Statistics, 22(3):1161–1176.
  * Loader, (2006) Loader, C. (2006).  Local regression and likelihood.  Springer Science & Business Media.
  * Loader, (1996) Loader, C. R. (1996).  Local likelihood density estimation.  The Annals of Statistics, 24(4):1602–1618.
  * Loftsgaarden and Quesenberry, (1965) Loftsgaarden, D. O. and Quesenberry, C. P. (1965).  A nonparametric estimate of a multivariate density function.  The Annals of Mathematical Statistics, 36(3):1049–1051.
  * Lorimer and Kramer, (2012) Lorimer, D. R. and Kramer, M. (2012).  Handbook of Pulsar Astronomy.
  * Lyon, (2016) Lyon, R. J. (2016).  Why are pulsars hard to find? PhD thesis, The University of Manchester (United Kingdom).
  * Ma and Li, (2019) Ma, H. and Li, J. (2019).  A true O$(n\log n)$ algorithm for the all-k-nearest-neighbors problem.  In International Conference on Combinatorial Optimization and Applications, pages 362–374. Springer.
  * Mack and Rosenblatt, (1979) Mack, Y. and Rosenblatt, M. (1979).  Multivariate k-nearest neighbor density estimates.  Journal of Multivariate Analysis, 9(1):1–15.
  * Mildenberger and Weinert, (2012) Mildenberger, T. and Weinert, H. (2012).  The benchden package: Benchmark densities for nonparametric density estimation.  Journal of Statistical Software, 46(14):1–14.
  * Miller and Dunson, (2019) Miller, J. W. and Dunson, D. B. (2019).  Robust Bayesian inference via coarsening.  Journal of the American Statistical Association, 114(527):1113–1125.
  * Newton, (2002) Newton, M. A. (2002).  On a nonparametric recursive estimator of the mixing distribution.  Sankhyā: The Indian Journal of Statistics, Series A, 64(2):306–322.
  * Newton and Zhang, (1999) Newton, M. A. and Zhang, Y. (1999).  A recursive algorithm for nonparametric analysis with missing data.  Biometrika, 86(1):15–26.
  * Ni et al., (2020) Ni, Y., Ji, Y., and Müller, P. (2020).  Consensus monte carlo for random subsets using shared anchors.  Journal of Computational and Graphical Statistics, 29(4):1–12.
  * Pólya, (1920) Pólya, G. (1920).  Über den zentralen grenzwertsatz der wahrscheinlichkeitsrechnung und das momentenproblem.  Mathematische Zeitschrift, 8(3-4):171–181.
  * R Core Team, (2018) R Core Team (2018).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rivoirard and Rousseau, (2012) Rivoirard, V. and Rousseau, J. (2012).  Bernstein–von Mises theorem for linear functionals of the density.  The Annals of Statistics, 40(3):1489–1523.
  * Rousseau and Mengersen, (2011) Rousseau, J. and Mengersen, K. (2011).  Asymptotic behaviour of the posterior distribution in overfitted mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(5):689–710.
  * Scott, (2015) Scott, D. W. (2015).  Multivariate Density Estimation: Theory, Practice, and Visualization.  John Wiley & Sons.
  * Sheather and Jones, (1991) Sheather, S. J. and Jones, M. C. (1991).  A reliable data-based bandwidth selection method for kernel density estimation.  Journal of the Royal Statistical Society: Series B (Methodological), 53(3):683–690.
  * Song et al., (2020) Song, H., Wang, Y., and Dunson, D. B. (2020).  Distributed Bayesian clustering using finite mixture of mixtures.  arXiv preprint, page arXiv:2003.13936.
  * Song, (2000) Song, P. X.-K. (2000).  Multivariate dispersion models generated from Gaussian copula.  Scandinavian Journal of Statistics, 27(2):305–320.
  * Terrell and Scott, (1992) Terrell, G. R. and Scott, D. W. (1992).  Variable kernel density estimation.  The Annals of Statistics, pages 1236–1265.
  * Teschl, (2009) Teschl, G. (2009).  Mathematical methods in quantum mechanics.  Graduate Studies in Mathematics, 99:106.
  * Tsybakov, (2004) Tsybakov, A. B. (2004).  Introduction to nonparametric estimation, 2009.  URL https://doi. org/10.1007/b13794. Revised and extended from the.
  * Vaidya, (1986) Vaidya, P. M. (1986).  An optimal algorithm for the all-nearest-neighbors problem.  In 27th Annual Symposium on Foundations of Computer Science, pages 117–122.
  * Wand and Jones, (1994) Wand, M. P. and Jones, M. C. (1994).  Multivariate plug-in bandwidth selection.  Computational Statistics, 9(2):97–116.
  * Wang and Dunson, (2011) Wang, L. and Dunson, D. B. (2011).  Fast Bayesian inference in Dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 20(1):196–216.
  * West, (1992) West, M. (1992).  Hyperparameter estimation in Dirichlet process mixture models.  Duke University ISDS Discussion Paper# 92-A03.
  * Wong and Ma, (2010) Wong, W. H. and Ma, L. (2010).  Optional Polya tree and Bayesian inference.  The Annals of Statistics, 38(3):1433–1459.
  * Zhang et al., (2014) Zhang, X., Nott, D. J., Yau, C., and Jasra, A. (2014).  A sequential algorithm for fast fitting of Dirichlet process mixture models.  Journal of Computational and Graphical Statistics, 23(4):1143–1162.


# COVID-19 Detection using Transfer Learning with Convolutional Neural Network

Pramit Dutta1, Tanny Roy2 and Nafisa Anjum3 Department of Electronics and
Telecommunication Engineering
Chittagong University of Engineering and Technology
Chittagong, Bangladesh
<EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

The Novel Coronavirus disease 2019 (COVID-19) is a fatal infectious disease,
first recognized in December 2019 in Wuhan, Hubei, China, and has gone on an
epidemic situation. Under these circumstances, it became more important to
detect COVID-19 in infected people. Nowadays, the testing kits are gradually
lessening in number compared to the number of infected population. Under
recent prevailing conditions, the diagnosis of lung disease by analyzing chest
CT (Computed Tomography) images has become an important tool for both
diagnosis and prophecy of COVID-19 patients. In this study, a Transfer
learning strategy (CNN) for detecting COVID-19 infection from CT images has
been proposed. In the proposed model, a multilayer Convolutional neural
network (CNN) with Transfer learning model Inception V3 has been designed.
Similar to CNN, it uses convolution and pooling to extract features, but this
transfer learning model contains weights of dataset Imagenet. Thus it can
detect features very effectively which gives it an upper hand for achieving
better accuracy.

###### Index Terms:

convolutional neural network, covid-19, ct-scan, f1 score, inception v3,
precision, recall

##  I Introduction

The COVID-19(SARS-CoV-2) was originally detected in China at the end of
December 2019. COVID-19 became an epidemic all over the world [1]. Statiscally
on November 9, 2020 , 508.6 million people were reported as COVID positive in
about 216 countries, the number of death crossed 1.2 million whereas more than
35.85 million have recovered [2]. This is a fatal respiratory disease, caused
by an acute respiratory syndrome of coronavirus or SARS-CoV-2 virus. According
to these statistics, the death rate has been increasing at an alarming rate.
Now as medical bio-molecular process, RT-PCR (Reverse Transcription Polymeras
Chain Reaction) is one of the standard and effective diagnostic method for
diagnosis of COVID-19. But this process is not only time-consuming but also
not effecient as RT-PCR is not widely accessible in different countries [1].
On the other hand, CT images are preferable for analysis of tissues of lung
muscles as X-RAY cannot effectively distinguish soft tissues in muscles [1].

In order to make the process of diagnosing the virus more effective, it has
become imperative to develop a time effecient method for the process. With the
rapid development of Artificial Intelligence, the Convolutional Neural
Networks (CNNs) have been justified as highly effective for the diagnosis of
many medical image processing. As such, it becomes easier to detect COVID-19
in Computed Tomography (CT) image. In the proposed model, a multi-layer
Convolutional neural network(CNN) has been used. Here, DTL(Deep Transfer
Learning) trains the image data by pre-trained model ,to specify, the
InceptionV3 model which contains weights from Imagenet.

The paper has been categorised into sections with section I introducing the
work done here. Section II contains the works related to this field and
Section III follows this with the mechanism for the working of proposed model.
Data-set, classifier model, image augmentation and hyper-parameter tuning are
the subsections included in section III. Section IV contains the result
obtained from the designed model. The training process, model evaluation as
well as comparison data depicting the effectiveness of proposed technique with
respect to existing techniques have been mentioned as subsections of section
IV. Section V concludes the paper and is followed by section VI that pertains
to the future work in this field.

Figure 1: Workflow Diagram of the Proposed Model

##  II Related Work

Throughout recent times, researchers have proposed various methods to detect
COVID-19 positive patients and to isolate them as soon as possible. One of
these works developed a DenseNet201 based deep transfer learning (DTL) to
classify the patients as COVID infected or not i.e. COVID-19 (+) or COVID-19
(-) [1]. Another one of these works has implemented a chest CT image
segmentation using a deep convolutional neural network (DCNN), a feature
variation block that adjusts the global parameters of the features for
segmenting COVID-19 infection was introduced here [3]. Another paper for
detecting COVID-19 has compared various deep learning model performances and
developed a model for X-ray, Ultrasound and CT images [4]. Another paper has
introduced a deep transfer learning technique and a top-2 smooth loss function
with cost sensitive attributes was also utilized for the classification of
COVID-19 infected patients [5]. In this paper, an image fusion machine
learning model was built by using a deep belief network model [6]. This paper
discusses on conventional two-dimensional ultrasound and color Doppler
ultrasound to observe the characteristics of lesions in non-critical COVID-19
patients [7]. This paper presents a novel, fully automatic approach for
classifying vessels from chest CT images into arteries and veins [8]. Another
paper discusses different CNN architectures and evaluates the influence of
data-set scale and spatial image context on performance [9].

##  III Methodology

The system diagram (Fig.1) shows how the whole classification process works.
Firstly, the training image was processed and augmented. Then, these images
were used to train the neural network built with a combination of transfer
learning and DNN model. After this step, the model was evaluated in the
validation set.

###  III-A Dataset

The COVID-19 CT scan images dataset was collected from kaggle [10]. In this
dataset, there were two classes for covid positive and covid negative. For
training, there were 279 images for covid positive and also 279 images for
covid negative. For validation, the dataset contained 70 images for covid
positive and the same amount of images for covid negative. It was 80% to 20%
split in the training set and validation set respectively.

###  III-B Classifier Model

In this paper, a deep learning-based neural network model has been proposed in
which there are 2 parts. The first part of this network contains a transfer
learning model named as Inception V3 and the second part is a customized deep
neural network (DNN) layer. This Inception V3 model contains weights of
imagenet dataset. The COVID-19 is a newly originated disease. As a result the
dataset of covid chest image is not large enough. So, transfer learning model
has been used to detect this disease. In this case, Inception V3 has been used
which has such filter and convolutional layer that can easily detect features.
That gave us a upper hand in the detection of COVID-19. The full structure of
this model was not used here rather only a small portion with a DNN was used
to extract the feature from that portion.

####  III-B1 Inception V3

The classification model was developed with the assistance of modified pre-
trained Inception-v3 transfer learning model. The Inception-v3 model could
factorize the convolution layer which could reduce the number of parameters
without even affecting the accuracy. Again, it concatenated the max-pooling
and convolutional layer which make feature downsizing more effective. The
model had the advantage to extract output from any particular concatenation
node. They were named as mixed layer and there are 11 mixed layer in total. In
order to make the model more effective for the experiment, the general
structure of the model was modified by using only 4 of them as shown in Fig.2.
Otherwise, it could cause overfitting due to our small dataset.

####  III-B2 DNN Model

To extract the output the last few layers of the model were replaced with a
DNN which used four customized layers whereas flatten used for transforming
mixed4 layer output into one-dimensional array, dense layer had 1024 neurons,
the next layer was used for 20% dropout. Finally, dense_1 layer with 1 neuron
with a sigmoid activation function was used.

Figure 2: Neural Network Model Structure

###  III-C Image Augmentation

For image detecting purposes it’s an efficient technique to learn more
features from image datasets from different dimensions. In this case, the
dataset was very small due to insufficient data. So, image augmentation
increased the diversity and size of the dataset helping the model to extract
unique features. For this purpose, operations like zoom, horizontal flip,
shear, shifting were performed on this dataset.

###  III-D Hyperparameter Tuning

In this model, the RMSprop optimizer was used with a 0.00003 learning rate and
batch size of 32. It was a binary classification problem. So, the loss
function was binary cross-entropy and metrics were precision and recall.

##  IV Result

###  IV-A Training Process

This model was trained in Google Colab using TensorFlow package. We used 1.66
GB RAM and 1.40 GB GPU. The model was trained for 31 epoches[callback stop]
taking about 176 seconds time(5 seconds per epoch).

###  IV-B Model Evaluation

The model performance was evaluated after every epoch by our validation
dataset. In the case of evaluation, four-parameter such as loss, accuracy,
precision, and recall were considered.

TABLE I: Loss & Accuracy Epoch |  Training Loss |  Validation Loss |  Training Accuracy |  Validation Accuracy
---|---|---|---|---
After Epoch 1 | 1.4701 | 0.6272 | 52.51% | 68.57%
After Epoch 6 | 0.5719 | 0.7031 | 67.92% | 58.87%
After Epoch 11 | 0.4897 | 0.4438 | 76.16% | 77.14%
After Epoch 16 | 0.4370 | 0.6054 | 80.65% | 70.00%
After Epoch 21 | 0.2887 | 0.4549 | 86.56% | 83.57%
After Epoch 26 | 0.3123 | 0.4597 | 86.38% | 81.43%
After Epoch 31 | 0.2073 | 0.4432 | 91.40% | 84.29%

After the training process, loss was 0.2073 and 0.4432 for training and
validation following (TABLE I). It can be also be observed from the graph
(Fig.3) how training and validation loss are decreasing over epochs. In this
case there no major dump(most of the cases between 0.9 to 0.4 for validation
set) can be observed during the training process.

Figure 3: Classifiaction Loss vs Epoch

On the other hand, after the training process, accuracy was 91.40% for
training and 84.29% for validation((TABLE I)). Though Fig.4 illustrates some
fluctuation during the training, they are small in size (most of the time it
fluctuates between 60% to 75% for validation set).In a word, our training
process was smooth and robust.

Figure 4: Classifiaction Accuracy vs Epoch

The model was also evaluated by precision and recall metrics in every epoch
and F1 Score for the final epoch was calculated which also gave a proper
evaluation on basis of both precision and recall.

TABLE II: Precision & Recall Factor |  After Epoch 1 |  After Epoch 6 |  After Epoch 11 |  After Epoch 16 |  After Epoch 21 |  After Epoch 26 |  After Epoch 31
---|---|---|---|---|---|---|---
Training Precision | 0.5263 | 0.6678 | 0.7664 | 0.8087 | 0.8778 | 0.8638 | 0.9140
Validation Precision | 0.6711 | 1.0000 | 0.8519 | 0.8500 | 0.7901 | 0.8548 | 0.8636
Training Recall | 0.5018 | 0.7133 | 0.7527 | 0.8029 | 0.8495 | 0.8638 | 0.9068
Validation Recall | 0.7286 | 0.1714 | 0.6571 | 0.4857 | 0.9143 | 0.7571 | 0.8143

In the case of precision, it was 0.9140 and 0.8636 for training and validation
following (TABLE II). On the other hand training recall was 0.9068 and
validation recall was 0.8143 (TABLE II). We know,

$F\textsubscript{1}Score=2\times\left(\frac{Precision\times
Recall}{Precision+Recall}\right)$ (1)

From this equation (1) and result from (TABLE II), 0.9104 was found as
training F1 Score and 0.8382 as validation F1 Score.

###  IV-C Comparison

MICHAEL. J. HORRY et al [4] discussed on different transfer learning model
performance on different types of data. In this proposed model the model
performance was evaluated on CT scan data. In the following table, the model
performance was compared with other models build for this classification task
and used CT scan images for the training and evaluation process.

TABLE III: Comparing Model Performance Model | Accuracy | Training Comment
---|---|---
VGG-16 | 79% |  Converged.Overfitting evident from 5 epochs.
VGG-19 | 78% |  Converged. Overfitting evident after 20 epochs.
Xception | 70% |  Did not converge. Overfitting evident immediately
InceptionResnet | 63% |  Did not converge. Overfitting evident after 1st epoch.
InceptionV3 | 71% |  Did not converge. Overfitting evident after 2 epochs.
NasNetLarge | 64% |  Did not converge. Overfitting evident immediately.
Densenet121 | 75% |  Converged. Overfitting evident after 8 epochs.
ResNet50V2 | 66% |  Converged poorly. Overfitting evident immediately.
Proposed Model | 84% |  Converged well.Overfitting not evident after 31 epochs.

From this table, it can be seen that the previously used VGG-16 model achieved
the highest accuracy among all the models working with a CT scan dataset while
the inception V3 model had 71% accuracy. So, our proposed model with 84%
accuracy outperformed all the models proposed in this classification task.

##  V Conclusion

In this paper, a deep learning concept of transfer learning to detect COVID-19
was proposed. This model shows that computer vision has the power to bring
about radical changes in the analysis of radiological images.Hence a time
efficient solution can be designed for identifying and isolating infected
patients. With a small dataset, the proposed model results in an exceptional
result with a validation accuracy of 84% comparing to 71% in the InceptionV3
model [3]. This model also outperforms all other models worked with CT Scan
data (previously 79%) (TABLE III). It’s because of the combination of
InceptionV3 model and the customized deep neural network(DNN) model that
enables this model to outperform all the model that have been proposed based
on CT scan images. Moreover, day by day the size of such dataset will increase
and so the model will be more accurate and robust.

##  VI Future Work

The future scope of the paper is to test and improve the model performance for
other lung diseases as well.

## References

  * [1] A. Jaiswal, N. Gianchandani, D. Singh, V. Kumar & M. Kaur, “Classification of the COVID-19 infected patients using DenseNet201 based deep transfer learning,“Journal of Biomolecular Structure and Dynamics, Jul 2020
  * [2] https://www.worldometers.info/coronavirus/
  * [3] Qingsen Yan et al, “COVID-19 Chest CT Image Segmentation – A Deep Convolutional Neural Network Solution“,arXiv:2004.10987v2 [eess.IV], April 2020
  * [4] MICHAEL. J. HORRY et al, “COVID-19 Detection through Transfer Learning using Multimodal Imaging Data,“ IEEE Access, Auguest 2020
  * [5] Y.Pathaka, P.K.Shuklab, A.Tiwaric, S.Stalind, S.Singhe & P.K.Shuklaf, “Deep Transfer Learning Based Classification Model for COVID-19 Disease,“ Elsevier Masson France,MAY 2020.
  * [6] M. Kaur, D. Singh, “Fusion of medical images using deep belief networks,“ Springer Science+Business Media, LLC, part of Springer Nature 2019, November 2019
  * [7] Yi Huang et al, “A preliminary study on the ultrasonic manifestations of peripulmonary lesions of non-critical novel coronavirus pneumonia (COVID-19)“SSRN, February, 2020
  * [8] Pietro Nardelli et al, “Pulmonary Artery-Vein Classification in CT Images Using Deep Learning,“ IEEE Transactions on Medical Imaging, in 2018.
  * [9] Hoo-Chang Shin et al, “Deep Convolutional Neural Networks for Computer-Aided Detection: CNN Architectures, Dataset Characteristics and Transfer Learning,“ IEEE TRANSACTIONS ON MEDICAL IMAGING, MAY 2016
  * [10] https://www.kaggle.com/andrewmvd/covid19-ct-scans