text:
\section{Introduction}
Cyber security can be broadly defined as the set of tasks and procedures required to defend computers and individuals from malicious attacks. Its origin can be traced back to 1971, a period where the Internet as we know it today was not even born. Amongst the computer science community it is widely accepted that it all started with Bob Thomas and his harmless experimental computer program known as the Creeper. This program was designed to move through the Advanced Research Projects Agency Network (ARPANET) leaving the following message: "I'm the creeper: catch me if you can". Inspired by Bob Thomas' Creeper, Roy Tomlinson created an enhanced version, allowing the Creeper to self-replicate, therefore coding the first computer worm. Later on, he would also design the Reaper which can be considered the first antivirus system, since it was designed to move across the ARPANET and delete the Creeper.

Despite a harmless origin, some years later the world would find out that network breaches and malicious activity were more dangerous than expected and cyber threats became a serious matter. Nowadays, cyber security is considered a major concern that affects people, organisations and governments equally, due not only to the growth of computer networks and Internet usage but also to the fact that cyber attacks are more sophisticated and frequent than ever. These cyber attacks represent a complex new challenge that demands more innovative solutions, hence, it requires a multi-disciplinary effort in order to be well-prepared and protected against such threats. Some of the disciplines involved in this task include computer science, computer and network architecture and statistics \citep{adams+heard}.

In this review we are mainly interested in the statistical approaches to cyber security problems and we centre our attention on how the discovery of cyber threats has been tackled as an anomaly detection problem. In particular, we discuss three different classes of anomaly detection problems within cyber security research and the typical types of data used in each one of them. The first of these problems is about volume-traffic anomaly detection, the second one about network anomaly detection and, finally, the third one is about malware detection and classification. For each of these problems we gently describe and discuss some of the most common Bayesian models that have been used to address them. For other approaches to anomaly detection and cyber security, such as classical statistics, machine learning and data mining approaches the reader could direct their attention to recent reviews by \cite{Buc+Guv}, \cite{Chandola}, \cite{Gupta+Gao} and \cite{adams+heard}.

Traditionally, cyber security threat detection systems have been built around signature-based methods; in this approach, large datasets of signatures of known malicious threats are developed and the network is constantly monitored to find appearances of such signatures. These systems have been proved effective for known threats but can be slow or ineffective when dealing with new ones, with mutations of known ones or with time-evolving threats. These are some of the reasons why we need to consider alternatives to signature-based methods. In order to do so, statistics offers a wide range of options for cyber security problems; these include both classical and Bayesian approaches that, in general, can be built on either parametric or nonparametric assumptions. Statistical anomaly detection methods usually build a model of normal behaviour to be considered as a benchmark, so that departures from this behaviour might be an indication that an anomaly has occurred.

Cyber security research from a mathematical and statistical point of view has proved to be an interesting and complex challenge that has led to an increasing interest in recent years. There are various reviews and reports \citep[see {\it e.g.}][]{Willinger+paxson, Catlett08,meza+juan,Sandia}, that outline some of the key areas, problems and challenges the mathematical community faces. It was early remarked  \citep{Willinger+paxson} how the constant changes in time and sites made the Internet such a difficult object to understand. Since then all authors have agreed that, due to the continuously exponential growth of the Internet and computer networks, there is a need for statistical models able to scale well to high-volume datasets of real time heteroscedastic and non-stationary data, which represents a significant theoretical and computational challenge. Moreover, as pointed out by \cite{Catlett08}, the mathematical models used should be able to effectively distinguish between harmless anomalies and malicious threats.

The need to design on-line detection methods able to handle high-volumes of data is not the only challenge discussed in these reviews. For example, \cite{Catlett08} also discussed the role mathematics play, by allowing us to understand computer networks, the Internet and malware behaviour in providing predictive awareness for secure systems. On a separate note, the author also remarked the need to advance the state of the art in graph theory and large-scale simulation to understand the spreading process of malicious code. \cite{meza+juan} further emphasised the importance of having access to reliable data. The lack of it is (mainly) due to privacy and confidentiality reasons and has made researchers study the best way to anonymise or sanitise the data. Some methods, as discussed in \cite{Bishop}, include using synthetic data, extracting the data from sources with no privacy constraints or a proper sanitising process. For example, in user-systems any characteristic that can be associated to an individual should be suitably changed (e.g. IP address or user name). More complex anonymization processes have also been developed, for example, in \cite{Tang_Wu} the authors describe a process based on subnet clustering where three parts of a whole IP address are anonymized by different methods. Fortunately, as we see in the following sections, there are some publicly available datasets that can be used for research purposes.

All of the other challenges just described can be well-grouped into three general cyber security research areas \citep{Sandia}. The first area deals with the modelling of large-scale dynamic networks like the modern Internet or any current computer network, that cannot longer be well-modelled with the classical graph theory formulation. Hence, there is a need to develop more sophisticated network mathematical formulations and new statistical techniques for comparing them. The reader could refer to \cite{olding} for a review on classical graph theory methods applied to modern network data. Discovering cyber threats is the second cyber security research area. As already established, cyber attacks are more sophisticated and frequent than ever, hence, the need for models capable of detecting malicious activity along with their variations, complicated multi-stage attacks and, if possible, the source of the cyber attack \citep{Sandia}. This is the are we explore in more depth in this review. Finally, the last research area is related to network dynamics and cyber attacks. This area is mainly dedicated to understanding the spreading characteristics of the malicious code through a computer network, before and after it has been detected and protections have been released. Particularly interesting problems consist in determining the potential limit of the infection and the interplay of the malicious spreading and the protection processes.

The remainder of the paper is organised as follows: in Sections 2, 3 and 4 we describe and explore, respectively, volume-traffic anomaly detection, network anomaly detection and malware detection and classification. In each of these sections we provide a gentle description of the kind of data used and how Bayesian models have been used to address these problems. In Section 5 we provide an insight into alternative cyber threat anomaly detection procedures. In Section 6 we describe some of emerging challenges. Lastly, Section 7 presents final points of discussion.

\section{Volume-traffic anomaly detection}

To begin developing statistical methods for computer network data, it is useful to have a high-level description of a computer network. The Open Systems Interconnect (OSI) is a widely-used conceptual set of rules for computer systems to be able to communicate with one another. The correct and reliable transmission of information is achieved through the joint work of seven sequentially connected layers, each one with its own purpose. As such, malicious activity could be targeted to any of the layers in order to destabilise the communication process between computer systems. For the purposes of this section we restrict our attention to the third layer of the OSI-model: the network layer. This layer is in charge of structuring and managing a multi-edge network including addressing, routing and traffic control \citep{hall2000internet}. The data is transmitted by breaking it down into pieces called packets that contain the user data (or payload) and the control information which provides data for delivering the payload, e.g. source and destination network addresses, error detection codes and segment information.

The packet rate which is defined as the number of packets per time unit moving across the network is one of the most common volume-traffic characteristics used for analysing a network's traffic. Their constant surveillance is useful for the detection of certain cyber threats that create changes in the network's normal traffic behaviour, such as distributed denial of service (DDoS) attacks which are intended to saturate the victim's network with traffic. Volume-traffic datasets can be obtained upon request from Los Angeles Network Data Exchange and Repository (LANDER) project. Another free network flow dataset is described in  \cite{akent-2015-enterprise-data}. The downloadable file "flows.txt.gz" presents network flow events from 58 consecutive days within Los Alamos National Laboratory's corporate internal computer network; each event is characterised by 9 variables: time, duration, source computer, source port, destination computer, destination port, protocol, packet count and byte count. The first three events included in the file are reported, as an illustration, in Table~\ref{FLOWS}.
\begin{table}[ht]\footnotesize{
\centering
\begin{tabular}{ccccccccc}
  \hline
time & duration & source comp. & source port & dest. comp& dest. port & prot. & packet count & byte\\
  \hline
1&0&C1065&389&C3799&N10451&6&10&5323 \\
1&0&C1423&N1136&C1707&N1&6&5&847  \\
1&0&C1423&N1142&C1707&N1&6&5&847  \\
   \hline
\end{tabular}
\caption{Extract form the network flow events (Los Alamos National Library).}
\label{FLOWS}}
\end{table}

In this dataset we can identify two different kind of variables. It is clear that we have access to volume-traffic characteristics such as the packet or byte count which can be used for volume-traffic anomaly detection purposes. The second set of variables characterise each event by providing the source and destination computer, the ports and the protocol used which allow us to perform a more refined analysis. This data can be used if we are trying to detect a rogue user or if there is a particularly sensitive and important computer that should not be involved in more than a certain amount of events or should not be the destination node for a certain kind of source computers. The above piece of information also allow us to develop multi-channel detectors by splitting the network traffic characteristics into separate bins represented by the source or destination. Moreover, a dynamic analysis of the network can be performed by taking into account the time at which each event took place.

\subsection{Bayesian approaches to volume-traffic anomaly detection}

Volume-traffic anomaly detection is concerned with detecting cyber attacks that produce changes in traffic measures such as the packet rate. The main goal is to detect as fast as possible changes in the normal behaviour. Once a change has been detected an alarm needs to be sent off so that the system can be checked and then decide whether there has been an attack or not (false alarm). It is important to remark that false alarms could yield important interruptions in the computer network, so there is a need to find the true change by seeking a low false positive rate as well. This yields a tradeoff between the detection delay and the false alarm rate that we need to consider for the detection procedure. The methods used to analyse these kind of data are mainly based on the statistical theory of change-point analysis.

\subsubsection{Change-point analysis}

The main objective of change-point analysis is the accurate detection of changes in a process or system that occur at unknown moments in time. In a single change-point setting we assume that there is a sequence of random variables $\{X_n\}_{n\geq1}$ (not necessarily independent nor identically distributed (iid)) such that $X_n\sim f(X_n|X^{(n-1)}),$ where $X^{(n-1)}=(X_1,...,X_{n-1})$ and the common probability density function (pdf) $f$ is known as the pre-change density. Then, at an unknown time $\nu$, something unusual occurs and from the time $\nu+1$ onwards $X_n\sim g(X_n|X^{(n-1)}).$ In this setting $\nu$ is known as the change-point and the pdf $g\neq f$ is called the post-change density. It is important to remark that theoretically, the densities $f$ and $g$ might depend on $n$ and $\nu$ as well, that is, $f(X_n|X^{(n-1)})=f_n(X_n|X^{(n-1)})$ and $g(X_n|X^{(n-1)})=g_{n,\nu}(X_n|X^{(n-1)}).$ Allowing these densities to depend on $n$ and $\nu$ might help us to more realistically explain time-evolving data found whilst doing cyber security research. In practice the post-change distribution $g_{n,\nu}(\cdot)$ might only be known up to some unknown parameter or vector of parameters $\theta$. Therefore, in some applications the problem can be reduced to detecting changes in mean, changes in variance or changes in both.

In order to deal with change-point detection problems, one could either follow a non-\linebreak sequential approach where the objective is to detect the changes in a fixed set of observations, or a sequential approach where the goal is to detect changes as new data arrives. Since both of these approaches have been tackled from a classical and a Bayesian perspective, the choice will certainly depend on the type of problem at hand and the objective of the analysis. From a cyber security point of view there is a need for constant surveillance of the computer network, therefore it is important to have fast on-line detection procedures. That is why in this review we only provide an insight into the sequential change-point analysis theory (for a complete review the reader can refer to \cite{Polunchenko_2011}) and how it has been applied to volume-traffic anomaly detection problems.

\subsubsection{Sequential change-point analysis}
The objective of a sequential change-point analysis is to decide after each new observation whether the observations' common pdf is currently $f$ or not. One of the main challenges is that the change should be detected with as few observations as possible past the true change-point, which must be balanced against the risk of false alarms. In order to address this, we must reach a tradeoff between the loss associated with the detection delay and the loss associated with false alarms \citep{polun+tartak}. A good sequential detection procedure must minimise the average detection delay (ADD) subject to a constraint on the false alarm rate (FAR).

There have been several approaches to analyse the tradeoff between early detection and number of false alarms. For the purposes of the review we focus our attention on the minimax formulation \citep{Lorden}, where $\nu$ is unknown (but not random) and the Bayesian formulation, where $\nu$ is considered a random variable with an appropriate prior distribution assigned. From a statistical perspective, at each step we need to test the hypothesis $H_k:\nu=k\geq0$ vs $H_{\infty}:\nu=\infty.$ Once a detection statistic, based on the likelihood ratio given by $\Lambda_n^k=\frac{p(X^{(n)}|H_k)}{p(X^{(n)}|H_{\infty})},$ has been chosen we supply it to an appropriate sequential detection procedure i.e. a stopping time $T$ with respect to the natural filtration $\mathscr{F}_n=\sigma(X^{(n)})$. For example, in the minimax setting the cumulative sum control chart (or simply CUSUM) procedure \citep{Page} is widely used, whereas in the Bayesian approach the most commonly used one is the Shiryaev-Roberts (SR) procedure \citep[][]{shiryaev_63,roberts}. Once a realisation of the detection procedure takes place at, say, time $T$, we have a false alarm if $T\leq\nu$ otherwise, the detection delay is given by the random variable $T-\nu$, which is commonly measured through the average detection delay given by $ADD(T)=\mathbb{E}(T-\nu|T>\nu)$.

A Bayesian optimal detection procedure will be the one that minimises the ADD subject to the probability of false alarm, $\mathbb{P}(T\leq \nu)$, does not exceed a fixed threshold $\alpha\in (0,1)$. It is important to remark that in both the minimax and Bayesian formulation the detection procedure is only applied once, hence, what happens after the alarm is not important. This certainly represents an issue for
surveillance applications, such as cyber security, where there is a need to repeatedly apply the detection procedure. After each false alarm or an actual attack we need to start monitoring the system again, hence a specification of a renewal mechanism is required. For example, assuming an homogeneous process, we can start from scratch after every alarm, yielding a multi-cyclic model \citep{tartak1}. In this scenario it is very likely that if an attack occurs it will happen a long time after the surveillance began, situation that has led to interesting theoretical results. For example, \cite{pollak} proved that under some assumptions the multi-cyclic SR detection procedure is optimal with respect to the stationary average delay to detection (SADD) which can be thought as the limiting case of the ADD.

In network security, change-point detection theory provides a natural framework for \linebreak volume-traffic anomaly detection. Still, there are some important considerations we need to contemplate. \cite{tartak1} argues that the behaviour of both pre- and post-attack traffic is poorly understood, as result, neither the pre- nor post-change distributions are known. Therefore, it is suggested that we can not rely on the likelihood ratio at all and we might think of using score-based statistics instead. The statistics proposed are suitable modifications of the CUSUM and SR statistics that include a score sensitive function. The way in which this score function is selected will depend on the type of change one is trying to detect. For example, for a change in mean, \citet{tartak06a,tartak06b} proposed a linear memoryless score function. Another important observation due to \citet{tartak06a,tartak06b} is that, in certain conditions, splitting packets in bins and considering multichannel detectors helps localise and detect attacks more quickly. This multichannel setting can be thought as a generalisation of the classical change-point detection problem, where we assume that an $n$-dimensional stochastic process is going to be observed simultaneously and at a random time only one of the entries will change its behaviour. This setting might be useful when dealing for example with DoS attacks where it has been observed that an increase number of packets of certain size occurs during the attack.

\subsection{Case study: ICMP reflector attack (LANDER project)}
The Internet Control Message Protocol (ICMP) reflector attack is a Denial of Service (DoS) attack that lasted 240 units of time. In this section we describe how the multi-cyclic Bayesian sequential $N(\mu,a\mu)-$to$-N(\theta,a\theta)$ change-point model \citep{polun+tartak} was used in order to detect it. In this model, the first Gaussian component describes the normal packet rate's behaviour and the second one the behaviour during the attack. Although the volume data is a time series of counts, the authors found that the Gaussian assumption was more realistic that Poisson process models used in other approaches to traffic data\citep[see {\it e.g.}][]{Cao2003, Kara_Molle}. A reproduction of the trace using the point estimates of the true data is showed in Figure~\ref{sim_packet_rate}. Due to the abrupt change in the trace, any good change-point detection procedure should be able to detect it. Hence, and in order to test the model and its capabilities, the authors transformed the original process for it to behave as a $N(\mu,a\mu)-$to$-N(\theta,a\theta)$. A reproduction of the transformed trace is illustrated in Figure~\ref{sim_packet_rate_2}.
\begin{figure}[h!]
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=3in]{sim_packet_rate.png}
  \caption{Reproduction of the original trace.}
  \label{sim_packet_rate}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=3in]{sim_packet_rate_2.png}
  \caption{Reproduction of the transformed trace.}
  \label{sim_packet_rate_2}
\end{subfigure}
\caption{Reproduction of ICMP reflector attack.}
\label{fig:ICMP}
\end{figure}

In order to find the threshold used as the benchmark for raising an alarm, the first step is to make sure that average run length (ARL) to false alarm rate is at a desired level $\gamma$. Setting this $\gamma=1000$ allows to use the asymptotic analysis described by the authors and found the threshold for the SR procedure at 731.3. Then the multi-cyclic procedure can be applied, so that each time an alarm is raised the process starts afresh. The results of this procedure are illustrated in Figure~\ref{SR}. It can be appreciated that the true attack was detected 23 seconds after it started and that the change-point model also raised a false alarm. It is clear that this model works reasonably well since the attack was detected (even though at first sight it was not visible) with only one false alarm. However, it has to be noted that the procedure was used on synthetic data. In the real trace, this procedure detected the true change-point but it also raised three false alarms.
\begin{figure}[h!]
  \begin{center}
    \includegraphics[width=3in]{SR.png}
    \caption{Multi-cyclic SR procedure on the packet rate illustrated in Figure~\ref{sim_packet_rate_2}.}
\label{SR}
  \end{center}
\end{figure}
\color{black}
\vspace{-20pt}
\section{Network anomaly detection}

Monitoring volume-traffic data is one important method of cyber security anomaly detection. However, there are other variables we can consider for network analysis. For example, monitoring the features that characterise each packet such as its length, the version, the header length, the priority, or characteristics of the connections between computers such as the source and destination computer or protocol can be used and need different anomaly detection methods. Most of the procedure used will require a graph theory formulation and classification models.

As established before, statistical anomaly detection models usually build a model of normal behaviour that works as a benchmark. Therefore, network anomaly detection models also allow us to analyse the normal behaviour of a computer network. In particular, studying the connections between the components of a computer network from a graph theory perspective has proved to be helpful for threats such as compromised credentials and intrusion detections. In research some of the most commonly used datasets for modelling computer network behaviour belong to Los Alamos National Laboratory (LANL) (\cite{hagberg-2014-credential}, \cite{akent-2015-enterprise-data} and \cite{turcotte+kent}). These datasets are mainly comprised of network and computer events collected from LANL enterprise network. For example, the User-Authentication Associations in Time dataset \citep{hagberg-2014-credential} encompasses 9 months of successful event authentications for a total of 708,304,516 connections. As an illustration, the first four events are shown in Table~\ref{auth_events}.
\begin{table}[ht]
\centering
\begin{tabular}{ccc}
  \hline
time & user & computer\\
  \hline
1&U1&C1\\
1&U1&C2  \\
2&U2&C3  \\
3&U3&C4  \\
   \hline
\end{tabular}
\caption{User-computer authentications associations in time.}
\label{auth_events}
\end{table}

This kind of datasets allow us to view the computer network as a bipartite graph with users and computers as nodes and the connections as edges. Some of the key tasks for anomaly detection consist in creating a set of normal pattern connections amongst users and computers and grouping together similar users.

\subsection{Bayesian approaches to network anomaly detection}
We now turn our attention to Bayesian network anomaly detection in the Bayesian framework. The basic idea is that we can characterise normal pattern connections within a computer network, either by the constant surveillance of packet features or the pattern connections between the nodes of the computer network. The detection of anomalous activity occurs when a connection can not be grouped with normal behaviour clusters. In the following subsections we explore some of the Bayesian approaches that have been used for network anomaly detection.

\subsubsection{Bayesian clustering}
The task of creating groups of normal behaviour can be achieved by cluster analysis. Distance-based algorithms are some of the most common ways to create such clusters; however, each one of these algorithms require the specification of what similarity between groups looks like through the choice of a distance measure (e.g. Euclidean distance or Manhattan distance). This implies that the clustering process will rely on both the algorithm and the distance (or similarity measure) choice. Besides using distances between data points, clustering can be performed by considering how dense parts of the data space are. In other words, clusters represent areas densely populated by data points, separated by regions of low density.

Bayesian clustering models have also been used for cyber security research. In particular, \cite{MeHe16} used a 2-step procedure for inferring cluster configurations of users with similar connection patterns and at the same time modelling new connections across the network. The first step uses a Bayesian agglomerative clustering algorithm with the choice of the multiplicative change in the posterior probability as a similarity measure. This algorithm yields an initial cluster configuration of users with similar connection behaviours, which is then used in a Bayesian Cox proportional hazards model \citep{Cox} with time-dependent covariates for the identification of new edges within the computer network. In this case, the chosen covariates for a connection between the $i$-th user and the $j$-th computer are the overall unique number of authentications over time for the computer and the restriction of these authentications to the user's cluster. This 2-step procedure requires a Markov Chain Monte Carlo (MCMC) algorithm for the joint update of the initial clusters and the coefficient parameters.

Working along this line, \cite{metelli2019} presented a Bayesian model for new edge
prediction and anomaly detection using a Bayesian Cox regression model like the one previously introduced in \cite{MeHe16}. However, in the most recent approach the authors used a more robust set of covariates and the initial cluster configuration was obtained through the spectral biclustering algorithm \citep{Dhillon}. The covariates used can be grouped into two different classes: the first group is comprised of the unique number of authentication over time for each client (time-varying out-degree), the unique number of authentication over time to each computer (time-varying in-degree) and two indicator functions respectively specifying if the last connection and the last two connections made by the client were new. The second set of covariates represent what the authors described as the notion of attraction between clients and servers. For their construction both hard-threshold and soft-threshold clustering models were used in a latent feature space.

\subsubsection{Bayesian networks}
A Bayesian network (BN) is a directed acyclic graph in which the nodes are the variables and the vertices represent the direct influences among the variables and their parent nodes \citep{Pearl_BN}. These influences are measured through the conditional probabilities and hence the model is completely characterised by them. For cyber security research purposes, BN's have been mainly used for intrusion detection. It is commonly argued that BN's yield robust models able to capture more realistic scenarios \citep{Frigault}. This is mainly due to the fact that BN's can directly model the combined effects of the vulnerabilities, contrary to other models where individual vulnerabilities are measured and then aggregated. In \cite{Kruegel} it was also discussed that using BN's might reduce the number of false positives that other anomaly detection models usually face.

One of the first approaches to intrusion detection through BN's can be found in \cite{valdes}, where the authors developed the eBayes TCP model in order to analyse temporally contiguous bursts of traffic at periodic intervals. In this model, the root node is the (unobserved) session class and the child nodes different (observed) variables, such as number of unique ports, service distribution and event intensity.     At each interval the idea is to know if an attack is taking place and the session class is assumed to propagate as a discrete Markov chain through the intervals. A similar approach for network packet traces can be found in \cite{shelton}, where a time continuous Markov chain is used instead. Other approaches for intrusion detection can be found in \cite{Kruegel}, where BN's are used to classify events as \textit{normal} or \textit{anomalous}, and in \cite{pauwels}, where the authors made an extension of dynamic Bayesian networks in order to model and detect anomalies on log files within the context of Business Processes, which are series of structured activities in order to perform a task.

Finally, another interesting use that researchers have given to BN's is the modelling of \textit{attack graphs} that represent how different network vulnerabilities could be combined in order to breach the network's security. This has been especially useful in order to measure and assess the risk associated to these vulnerabilities and therefore, for risk management \citep[see {\it e.g.}][]{Dantu, Frigault, Dewri}.

\subsubsection{Topic modelling and LDA}
Topic models are probabilistic models that have their origins in latent semantic indexing (LSI) modelling \citep{Deerwester90}. LSI uses the singular value decomposition of a \linebreak term-document association matrix in order to create a space where related terms and documents are placed near one another. \cite{hofmann}, using this idea, developed the probabilistic latent semantic analysis (PLSA) which can be considered the first topic model. His approach uses a generative latent class model to perform a probabilistic mixture decomposition. Other models, such as the latent Dirichlet allocation (LDA) model \citep{blei+lda}, have also been introduced with the task of discovering the topics that occur in a set of documents. However, they have been widely used in other fields where there is a need of unsupervised clustering.

The topics produced by these models are clusters of similar words which allow examining a set of documents. As discussed in \cite{blei+lda} two of the usual assumptions made are that the words in a document are exchangeable (also known as the \enquote{bag-of-words} assumption) and so are the documents. These assumptions allows us to exploit de Finetti's representation theorem for exchangeable random variables \citep{hewitt}. In the LDA model, the basic idea is that every document in the corpus can be represented as a random mixture over a known and fixed number $k$ of latent topics, which are characterised by a distribution over words. It is further assumed that the word probabilities are characterised by an unknown but fixed matrix $\beta$ of dimensions $k\times M$ that needs to be estimated, where $M$ is the size of the vocabulary. In practice $M$ is usually large, thereby creating issues related to sparsity and with the prediction of new documents. In order to address this, \cite{blei+lda} also developed a fuller Bayesian approach, usually called the smoothed LDA, by allowing $\beta$ to be random and a Dirichlet prior distribution assigned to each row $\beta_i$.

The LDA model's original setup includes a corpus $X$ with $N$ documents $\textbf{w}_1$,...,$\textbf{w}_N$, each one having $P_n$ words, $w_{1,1},...,w_{n,P_n}$. The length of each document can be sampled from a Poisson distribution or from a more realistic document-length distribution. For each of the words in the $n$-th document we first select a random topic $z_i$ from which a random word will be assigned. The mathematical representation of the generative process is:
\begin{align*}
\theta_n&\sim \text{Dirichlet}(\alpha)&n\in\{1,...,N\}\\
z_i&\sim \text{Multinomial}(\theta)&i\in\{1,...,P_n\}\\
w_{i}&\sim\text{Multinomial}(\beta_{z_i}).
\end{align*}

The LDA model has been shown to be a valid and useful technique for network anomaly detection and for modelling computer network behaviour. The reader can refer to \cite{HePaSk16} for an example on how the LDA model can be used to analyse computer network connection traffic data to determine the number of users present. In their approach, each document is represented by the day's authentication records, different users are the topics and the destination computers play the role of the words. In this scenario each of the entries of $\theta_n$ indicates how active was the respective user in the $n$-th day. As discussed by the authors this procedure could play an important role for detecting misuse of credentials.

In matters of anomaly detection, \cite{CaChLiFu16} used the LDA model to analyse features obtained from the packet headers captured using the \textit{tcpdump} software. In their approach, the documents are represented by the tcpdump traffic obtained within a time slot and the words are the unique packet's network features. The LDA model is used on free attack traffic data in order to learn its feature patterns and new traffic data is then compared against it. The authors proposed using the likelihood of a new document as the detector of anomalous activity. A similar procedure can be found in \cite{CrCa11}, where the LDA model and the dynamic LDA (dLDA) \citep{dlda} are considered to analyse Ethernet packets. In their approach the data is also divided into fixed time intervals and the words are any event of interest observed across 45 well-known ports used for network topic modelling. In their results the dLDA proved to be a better choice for modelling this kind of data due to the dLDA's ability to analyse time-dependent documents by letting the weights over the topics to change in time.

\subsubsection{Poisson factorisation}
Topic models are not the only probabilistic models used for cyber security research that have been originally designed for other purposes. Poisson factorisation (PF) models, which are widely used for recommender systems in machine learning \citep[see {\it e.g.}][]{Gopalan2014}, have also been used for network anomaly detection. PF is a probabilistic model of users and items that was proposed as an alternative to the classical probabilistic matrix factorisation (PMF) \citep{Sala_PMF}.

The main assumption is that the data can be represented in a matrix, in the recommender system the rows are the clients and the columns the number of items. Each entry of this matrix is assumed to be the rating given by a certain user to a particular item, and these are modelled using the dot-product of latent factors for both the users and the items. PMF assumes that each entry is normally distributed and Gaussian priors are assigned on the latent factors for both the users and movies. This will theoretically imply that the ratings could become negative which is something we would not desire to have. In order to address this issue, PF assumes that both the users and the latent factors are non-negative and so a Poisson distribution for the entries and gamma distributions for the latent factors are used instead. These assumptions make PF more applicable to real datasets like The Netflix Prize dataset \citep{netflix}, where we have a set of users and the rankings for each movie in the catalogue.

With respect to cyber security research, \cite{TuMoHeMc16} considered PF models for peer-based user analysis which provides a better understanding of the individuals by learning their peer's behaviour. The basic idea is that computer users with similar roles within an organisation will have similar patterns of behaviour. This type of analysis can be particularly important for quickly detecting rogue users. The behaviour of a new user  can be compared to their peers and anomalies detected. The data used is the recorded authentication events (Table~\ref{auth_events}) as briefly explained in Section 2. The model is completely specified by letting $Y_{ui}$ be the number of times that user $u$ authenticates on machine $i$ and where it is assumed that $Y_{ui} \sim\mbox{Poisson}(\theta_u \beta_i)$, where $\theta_u$ for $u=1,\dots,U$ and $\beta_i$ for $i=1,\dots,C$ are $k$-dimensional vectors of positive values. The model is interpreted as having $k$ latent features characterising the users (such as job title, department, etc.) with $\theta_u$ representing their scores for the $u$-th user and $k$ latent features characterising each computer (such as the number of daily processes, the type of computer, etc.), with $\beta_i$ representing their scores for the $i$-th computer.

A certain feature might have a high score for all machines within one department and low scores otherwise. If a user had a high-score on that feature then they are likely to have many authentication events on machines in that department (perhaps, representing that they work in that department). If a user had a low-score on that feature then they are likely to have a very low number of authentication events. In general, the mean number of authentication events for a user on a machine is the sum over products of many features which allows similarities between users and computers to be learnt from the data. The specification of the model is completed by assuming that
$\theta_{u,j}\stackrel{i.i.d.}{\sim}\mbox{Ga}(a, \zeta_u)$, $\beta_{u,j}\stackrel{i.i.d.}{\sim}\mbox{Ga}(b, \eta_i)$, $\zeta_u\sim\mbox{Ga}(a',b')$ and $\eta_i\sim\mbox{Ga}(c',d')$. The model is fitted to a training sample and anomalies can be detected by comparing predictions from this model to observed values from a testing sample.

\subsubsection{Dirichlet process}
So far, we have only described Bayesian models working under parametric assumptions. Although these models have been proved effective in detecting network anomalies, there is still a missing methodology that we need to consider: the nonparametric realm. Bayesian nonparametric models have become increasingly popular and appealing in research areas such as finance, biology and machine learning, amongst others, because they posses a flexibility that can hardly be achieved by parametric models. This flexibility is explained through the fact that these models in contrast to parametric ones, assume that the parameter space is an infinite-dimensional object, something particularly useful as more data becomes available. For cyber security research, where huge amounts of time evolving data are available almost instantly, Bayesian nonparametrics models should provide a more suitable modelling techniques and provide interesting insights and solutions to the problem at hand.

Historically, it is widely accepted that Bayesian nonparametrics had its beginnings with the introduction of the Dirichlet Process (DP) by \cite{ferguson73} and since then, the DP has played a vital role in Bayesian nonparametrics and its applications \citep[see {\it e.g.}][]{hjort2010}. The DP works as a prior on the space of probability distributions and just as the Dirichlet distribution we will preserve a nice conjugacy property. More precisely, if $\{X_i\}_{i=1}^n\stackrel{\text{i.i.d}}{\sim} P$ and $P\sim DP(\alpha)$ then $P|X_1,...,X_n\sim \text{DP}\left(\alpha+\sum_{i=1}^n\delta_{X_i}\right)$. In this setting, $\alpha$ is a finite measure defined on the same space as $P$ and it is commonly expressed as $\alpha=\theta P_0$, where $\theta$ is the total mass and $P_0$ a probability measure. Using this structure we obtain a nice expression for the predictive distribution
\begin{eqnarray*}
X_{n+1}|X_1,...,X_n&\sim&\begin{cases}
\delta_{X_j^*} \hspace{.5cm}\text{with probability proportional to } \hspace{.1cm}n_j\\
P_0\hspace{.7cm}\text{with probability proportional to } \hspace{.1cm}\theta,
\end{cases}
\end{eqnarray*}
where $\{X_j^*\}_j$ is the set of distinct values observed in $X_1,...,X_n$ and $n_j$ their frequencies. Hence, we will be observing either a previously seen value or a completely new one, a characteristic which is really useful in clustering applications.

Exploiting the structure of the posterior and predictive distribution of the DP, \cite{HeRD16} developed a Bayesian nonparametric approach to intrusion detection by assuming a DP-based model for each message recipient on a set of computers and the directed connections amongst them that represent the node set $V$ and the set of edges $E$ respectively in a directed graph $(V,E)$, which can be thought as a set of objects connected together where the source node and destination node can be identified for each connection (the direction matters). The first step of their anomaly detection procedure is to obtain the predictive $p$-value for the event $x_{n+1}$ defined as $ p_{n+1}=\sum_{x\in V: \theta_x^*\leq \theta_{x_{n+1}}^*} \frac{\theta_x^*}{\theta^*}$, where $\theta_x^*=\theta P_0(x)+\sum_{i=1}^n\delta_{X_i}(x)$ and $\theta^*=\theta+n$. These $p$-values quantify the level of surprise of a new connection. Since the goal is to detect anomalies in each source computer, the $m$ $p$-values observed in the edge $(x,y)$ are reduced to a single score using Tippett's method \citep{tippett}. Then a single score for each node $x$ is obtained using Fisher's method \citep{fisher}. Finally, the computers are ranked through these scores, and compromised ones should have higher ranks.

Working along this line,  \cite{Sanna_Heard} examined a joint model of a sequence of computer network links $\{(x_i,y_i)\}_{i=1}^n$, with $x_i$ and $y_i$ representing the source and destination computer respectively, based on the Pitman-Yor process (PY) \citep{Perman}. The PY process, also known as the two-parameter Poisson-Dirichlet process, requires two parameters which are usually denoted by $\sigma\in[0,1)$ (the discount parameter) and $\theta>-\sigma$ (the strength parameter). An appealing characteristic of the PY process is the closed form of the predictive distribution given by:
\begin{eqnarray*}
X_{n+1}|X_1,...,X_n&\sim&\begin{cases}
\delta_{X_j^*} \hspace{.5cm}\text{with probability proportional to } \hspace{.1cm}(n_j-\sigma)\\
P_0\hspace{.7cm}\text{with probability proportional to } \hspace{.1cm}\theta+k\sigma,
\end{cases}
\end{eqnarray*}
where $k$ is the number of different observations and $P_0$, $n_j$ and $X_j^*$ are defined just as for the DP. We can immediately notice that if $\sigma\rightarrow0$ we recover the DP, thus, the PY process can be thought as a generalisation of the DP. Another interesting remark is that the probability of a new value depends on $k$, so the more unique observations we have, the more likely it will be to obtain new samples from $P_0$. This is certainly useful for applications where power-law behaviour is observed, something the DP can not achieve. Finally, as $\sigma\rightarrow1$, unique observations $X_j$'s having small frequencies $n_j$ are unlikely to be sampled.

In \cite{Sanna_Heard} the joint modelling of $p(x,y)$ is achieved through the decomposition $p(x|y)p(y)$ by assuming the sequence of destination nodes, $\{y_i\}$, to be exchangeable with a hierarchical PY distribution and conditioned on the destination node the sequence of source nodes connecting to that destination, $\{x_i|y\}$, are also exchangeable with a hierarchical PY distribution with parameters depending on $y$. The mathematical representation of the model is:
\begin{align*}
x|y&\sim F_{x|y_i}\\
y_i&\sim G\\
F_{x|y}&\sim \text{PY}(\alpha_y,\theta_y,F_0)\\
G&\sim\text{PY}(\alpha,\theta,G_0)
\end{align*}
As for the detection procedure, it follows the same reasoning as in \cite{HeRD16}, that is, for each source computer one needs to obtain the predictive $p$-values and combine them into a single score. Besides the use of the PY process rather than the DP, there are two other interesting results found in this approach. The first one is that the authors do not restrict their attention to the use of $p$-values and they explore the use of mid $p$-values (the reader can direct the attention to \cite{lancaster} and \cite{delanchy_heard} for an insight on mid $p$-values and why in some problems they might be preferred over $p$-values). Finally, they also explored Pearson \citep{Pearsons} and Stouffer's \citep{stouffer} $p$-value combiners.

\subsection{Case study. LANL user-authentication dataset}
The user-authentication dataset by Los Angeles National Laboratory \citep{kent-2015-cyberdata1} is a comprehensive summary of 58 days of traffic of the enterprise's network, which also contains a set of 48,079 red team compromise events resulting from a simulated intrusion attack. In total, there are four source computers compromised (C17693, C18025, C19932, and C22409). The task is to identify the unusual activity occurring and flag these computers. These connections between source and destination nodes can be modelled as a bipartite graph and hence, we could use the model proposed by \cite{HeRD16} or \cite{Sanna_Heard} in order to identify the compromised computers.
As described in Section 3.1.5, these models allow us to quantify the surprise of each new connection through the $p$-values or mid $p$-values, which are then aggregated for each source computer to have a unique score. These scores are then used to rank the computers and, in the case of the compromised ones, a high anomaly score should be assigned. In Table~\ref{case_study} we present the results that \cite{HeRD16} obtained with the DP $p$-values aggregated using Fisher's method \citep{fisher} and the best results of \cite{Sanna_Heard} obtained with the PY mid $p$-values aggregated using Pearson's method \citep{Pearsons}.
\begin{table}[ht]
\centering
\begin{tabular}{ccc}
  \hline
Source Computer & Rank with DP $p$-values & Rank with PY mid $p$-values\\
  \hline
C17693&5&1\\
C18025&94&74  \\
C19932&5347&2754  \\
C22409&7172&6984 \\
   \hline
\end{tabular}
\caption{Rank of the compromised source computers.}
\label{case_study}
\end{table}
\vspace{-20pt}
\section{Malware detection and classification}

The malware detection and classification problem is the third group of cyber threat investigation problems that we consider in this review. A malware is defined as a software specifically designed to disrupt, damage or gain access to a computer system. Nowadays there are many types of malware, such as spyware, adware, ransomware, amongst others, including several variations of them. That is why the fast detection of unknown malware is one of the biggest concerns of cyber security. However, accurate detection is not the only task required when dealing with malicious software. Malware have to be classified into families for a better understanding of how they infect computers, their threat level and therefore, how to be protected against them. Correct classification of new malware into known families may also speed-up the process of reverse-engineering to fix computer systems that were infected. In order to have good explanatory and predictive models for the malware detection and classification problem, researchers have mainly used the content of the malware in two ways, either by using the hexadecimal representation of the binary code or the dynamic trace.

\subsection{Hexadecimal representation and $n$-grams}
In computer science the byte is the basic unit of information for storage and processing and it is most commonly represented by a sequence of 8 binary digits (bits). Every instruction given to a computer (malicious or not) can be broken down into sequences of bytes, which form the instruction's binary code. These sequences can be expressed in a more compact way through their hexadecimal representation, where each byte is written as a combination of two elements of the set $\{0,1,...,9,A,B,...,F\}$. For malware detection and classification the hexadecimal representation of the binary code is used in conjunction to what is known in literature as $n$-grams.

An $n$-gram is usually defined as a contiguous sequence of $n$ elements; these elements, depending on the field, can be letters, numbers, words, etc. In computer science and for malware detection and classification purposes, these elements are the hexadecimal representation of each byte. For example, considering $n=4$, for the code extract given by $00 \hspace{.1cm} 00  \hspace{.1cm} 1C  \hspace{.1cm} 40  \hspace{.1cm} 2A  \hspace{.1cm} 28$, we take all the possible sequences of 4 contiguous bytes to create the set of $4$-grams, i.e.\linebreak $\{00 00 1C 40, 00 1C 40 2A, 1C 40 2A 28\}$. The elements of the set of all the different $n$-grams are assumed to completely characterise benign and malicious code through their presence and absence. Hence, the set of $n$-grams needs to be modelled for the prediction of new observations, contrary to the classical $n$-gram analysis where these structures are used as a probabilistic model for the prediction of the next item in a sequence. In practice, it has been observed that the set of unique $n$-grams is huge even for small datasets (easily reaching the order of millions), so it becomes computationally unfeasible to use all this information and a feature selection procedure is required as a first step.

To the best of our knowledge, there is no Bayesian approach to malware detection using the hexadecimal representation of the binary code. However, we believe this line of work should be appealing to the Bayesian community as well due to assumptions made and the challenge it represents. First of all, the data can be represented by a $N\times M$ binary matrix where the rows are either malicious and benign code. From a Bayesian perspective one could be interested in modelling this binary matrix using a conjugate beta-Bernoulli hierarchical model or by choosing more complex models like a logistic regression for malware detection. It should also be clear that this is a high-dimensional problem in both rows and columns, so Bayesian nonparametric models could be definitely helpful for analysing this sort of data. Moreover, the assumptions made can be easily generalised for the classification into known families problem.

\subsection{Dynamic traces as Markovian structures}
The second approach to malware detection and classification is to use the dynamic traces of the malware, which is basically the set of instructions executed in order to infect the system. Following the ideas of \cite{Storlie_2014}, many authors have assumed that these traces have a Markovian structure and in that way the interest relies in modelling and analysing the probability transition matrix. Since there are hundreds of commonly used instructions and thousand of them overall, modelling the one-to-one transition is not feasible. However, there are some instructions that perform the same or similar task so creating groups of similar instructions is a reasonable first step. \cite{Storlie_2014} developed four different categorisations with $8, 56, 86 $ and $122$ groups of similar instructions. In practice the most most widely used categorisation is the one with 8 classes, which include amongst others: math, memory, stack, other.

The mathematical framework is fully specified by letting $c$ to be the number of instruction categories previously chosen (e.g. 8), then a dynamic trace is defined as a sequence $\{x_1,...,x_n\}$ with $x_i\in\{1,...,c\}$ that are modelled as a Markov chain (MC). Hence, we let $\bold{Z}_i$ be the transition counts matrix for the $i$-th program, $\bold{P}_i$ to be the probability transition matrix and $B_i$ an indicator if the program is malicious or not. \cite{Storlie_2014} proposed that the entries of the estimated $\bold{P}_i$, denoted by $\hat{\bold{P}}_i$, should be used as predictors to classify a program as malicious or not through a logistic spline regression model. In practice, the actual predictors used to model $B_i$ are: $\text{logit}(\hat{P}_{i,1,1})$, $\text{logit}(\hat{P}_{i,1,2})$, ..., $\text{logit}(\hat{P}_{i,c,c-1}),\text{logit}(\hat{P}_{i,c,c})$. Finally, a symmetric prior Dirichlet distribution with parameter $\nu$ is used as a prior for each row of $\bold{P}_i$.

Working directly on this approach, \cite{Kao+Reich} proposed a flexible Bayesian nonparametric approach to modelling the probability transition matrices $\{\bold{P}_i\}_i$ by using a mixture of Dirichlet processes (MDP) \citep{antoniak}, that is,
\begin{eqnarray*}
\bold{P}_i|\bold{Q}_i,\sigma&\stackrel{\text{i.i.d.}}{\sim}& MD(\sigma\bold{Q}_i)\\
\bold{Q}_i|G&\stackrel{\text{i.i.d.}}{\sim}& G\\
G&\sim& DP(\theta,P_0),
\end{eqnarray*}
where $P_0$ follows a matrix Dirichlet (MD) distribution centered in some constant matrix $\hat{P}$, $\sigma>0$ controls the variance of $\bold{P}_i$ and the matrix $\bold{Q}_i$ is the shape parameter. This MD distribution implies that each row of $\bold{P}_i$ independently follows a $c$-dimensional Dirichlet distribution with concentration parameter equal to the corresponding row of $\sigma\bold{Q}_i$. The model is completely specified by letting $B_i$ be the indicator random variable of maliciousness just like before. A new program $i^*$ is classified as malicious if $\mathbb{P}(B_{i^*}=1|\bold{Z_{i^*}},\mathscr{Z})$ exceeds a predefined threshold, with $\mathscr{Z}$ being the collection of all observed counts matrices. Moreover, if the program is malicious it can be further classified into a cluster with existing programs that share common features.

A different approach for the modelling of dynamic traces was developed by \cite{Bolton+Heard}. Their approach followed the same assumptions, i.e. that dynamic traces, specified by the prior clustering of common instructions, have a Markovian structure. However, they further assumed that this structure changes over time with recurrent regimes of transition patterns. Hence, each dynamic trace can be modelled as a MC with a time varying transition probability matrix $\bold{P}(t)$. In order to detect the regime changes three change-point models were described. The basic idea is that there are $k\geq0$ change-points that partition the dynamic trace and where within each segment the trace follows a homogeneous MC. The methods vary in the way the probability transition matrices are defined within each segment. The first one changes the whole matrix in each segment, the second one only allows some of the rows to change and finally the regime switching method allows the change in rows not only to be forward but also to go back to a vector of probabilities that governed the Markov chain in earlier segments. Finally, the authors proposed a classification procedure based on a similarity measure of the vectors of change-points and their regimes, that obtains the minimum value for the two samples of the proportions of instructions which occur within regimes shared by both traces. A high level of similarity requires that a large number of observations in both traces are drawn from common regimes.

\subsection{Case study: Microsoft malware classification challenge}
The Microsoft Malware Classification Challenge dataset \citep{MMCC} contains almost 11,000 malware representing a mix of 9 known families of malicious threats. For each file the hexadecimal representation of the binary code can be obtained and therefore, we could obtain the $n$-grams that characterise each family. For example, using a subset of 590 malware we are able to extract 826 $4$-grams appearing at least once in each family (Figure~\ref{malware_class}). It is important to note that other feature extraction processes could be considered.
\begin{figure}[h!]
  \begin{center}
    \includegraphics[width=3in]{malware_class.png}
    \caption{Graphical representation of the binary matrix comprised of 590 malware divided in 9 families and 826 binary features}
\label{malware_class}
  \end{center}
\end{figure}

Using this binary matrix as the training set, supervised learning algorithms such as naive Bayes, multinomial-logistic regression and decision tree (with its adaptive, gradient and extreme gradient boosted versions), can be used in order to classify new malware. The classification performance of the above classifiers is showed in Table~\ref{accuracy}.
\begin{table}[ht]
\centering
\begin{tabular}{cccccc}
  \hline
naive Bayes & Decision Tree & ABoost & GBoost & XGBoost & Mult. Log. Reg.\\
  \hline
77.38&80.15&94.05&91.67&91.27&84.92\\
   \hline
\end{tabular}
\caption{Classification accuracy (in$\%$) comparison.}
\label{accuracy}
\end{table}
\vspace{-20pt}
\section{Alternative cyber threat anomaly detection approaches}

It is imperative to stress that the models and the kind of data described in the previous sections are far from being exhaustive. They represent the ones that we have found to be the most frequently used for the general class of cyber threat anomaly detection problems presented here. However, there are other kind of cyber security related problems that have been tackled from a Bayesian perspective and that could be appealing for future research purposes.

For example, \cite{melissa_turcotte} used computer event logs to identify misuse of credentials within a computer network. In their approach these logs are treated as an aggregated multivariate data stream comprised of the client computer $x$, the server computer $y$ and the type of event $e$. These features are modelled independently for each user credential and the probability of an observed triplet $(x_t,y_t,e_t)$ at time $t$ is modelled using the conditional probabilities, that is, $\mathbb{P}(x_t,y_t,e_t)=\mathbb{P}(x_t)\mathbb{P}(y_t|x_t)\mathbb{P}(e_t|x_t,y_t)$. For each of these components an appropriate multinomial-Dirichlet based model is used, and in order to detect anomalies, the predictive distributions are used to obtain the $p$-value that can be compared against a predefined threshold. A second example for detecting compromised credentials can be found in \cite{price_w}, where the authors proposed a users' activity anomaly detection approach by analysing the amount of user activity on a given day and the times where these activities were realised. A seasonal behaviour model is developed by first constructing a model to measure the user's activity in a certain period and then using the events registered, a change-point density estimation model is used to estimate the times at which the events occurred. In this setting, users working at hours that differ from their normal schedule are considered to be anomalous activities.

Other interesting approaches are aimed to obtain a better understanding of a computer network's behaviour. An example of this can be found in \cite{price_w_2}, where the authors main goal is to detect automated events that can be viewed as polling behaviour from an opening event originated by a user. It is discussed that achieving this should yield an improvement in the statistical model used and enhance its anomaly detection capabilities. In this approach a change-point model is used in each edge of the computer network in order to separate human behaviour from automated events. This methodology works as an alternative to the one presented in \cite{6975589}, where the discrete Fourier transform is used.

Being able to detect automated events from human activity is not the only approach undertaken for a better understanding of a computer network. Nowadays, there are computer networks that contain a vast number of nodes and thereby, a large number of connections amongst them. In practice, temporal independence for the nodes is usually assumed to have mathematical tractability; however, this is a strong assumption and a deeper understanding of the dependance amongst these nodes is required. A recent approach by \cite{price_w_3} aims to detect and understand the interactions between computer nodes in order to detect correlated traffic patterns to reduce false positives when performing anomaly detection. A test based on higher criticism \citep{donoho2004} is used to detect this dependence.

\section{New and emerging challenges}
Some of the challenges the mathematical community face when dealing with cyber security problems have already been described in the Introduction. These include challenges related to the data itself, like the privacy and ethical issues, and to the models and how there is a need to handle large volumes of non-homogeneuous data. However, we would like to describe two of the main challenges that still need to be fully considered in Bayesian cyber security research.
\subsection{Robustness}
As already established, in cyber security an anomalous activity might be a sign that an attack is occurring, so it is imperative to detect it as fast as possible while keeping a low false positive rate. Clearly, the performance of most of the anomaly detection models heavily rely on the data used in the training step. Ideally, this data should be as reliable as possible and among other things, it should be noise-free. However, in real case scenarios, and especially nowadays with larger and more complex networks, delivering noise-free data might not be an easy task to achieve. Therefore, designing robust models able to deal with noisy data and to capture more complex (and realistic) scenarios is also a crucial task.

From a non-Bayesian perspective, there are already some robust anomaly detection models that have been used in cyber security. For example, \cite{Eskin} described a probabilistic approach to anomaly detection without training on normal data in order to detect intrusions in UNIX system call traces; \cite{paffenroth} developed a robust principal component analysis (RPCA) assuming noisy and missing data on network packets; \cite{RSVM} used robust support vector machines (RSVM) to study intrusion detection over noisy data. Finally, it is worth mentioning that deep learning models have also been used in order to provide more robust models, like the deep autoencoders which are a class of unsupervised neural networks \citep{BermanDS}.

From a Bayesian perspective, robust models applied to cyber security data are still little explored. Some of the already known robust models include Bayesian networks which are able to model more realistic scenarios by analysing the combined effect of the vulnerabilities. For volume-type traffic data, there is also a need of more robust models, since this type of data usually contains outliers, deriving from normal activities. In this direction, but not directly applied to cyber data, \cite{Knoblauch} developed a robust Bayesian online change-point detection algorithm that achieved promising results when dealing with outliers in well-log data and analysing noisy measurements of nitrogen oxide levels.

\subsection{Scalability}
Bayesian models have become appealing due to their rich theoretical background and ability to model complex data. Bayesian inference relies on the posterior distribution of the parameters which, in most of the cases, will not be known in a closed form. In order to obtain samples from the posterior, one could use approximations to the unknown integral or to the posterior, with the intent of minimising the discrepancy that forms the basis of variational Bayes. Alternatively, in some instances an MCMC scheme could be designed to produce correlated samples from the posterior. Although theoretically efficient, MCMC schemes often face several computational issues due to a slow mixing and slow convergence, which usually get worse when dealing with high-dimensional data. Possible solutions include, parallel MCMC, approximate MCMC, C-Bayes and Hybrid algorithms. A thorough review on theoretical and practical aspects of scalable Bayesian models can be found in \cite{scalingBay}.

From a Bayesian perspective, some scalable approaches have been designed for modelling and detecting anomalies in cyber security applications. Examples of these include: \cite{Clausen}, where a Markov-modulated Poisson process embedded in a fast and scalable Bayesian framework was used in the modelling for network flow data; \cite{scalableDN}, where a novel class of Bayesian dynamic models was introduced and applied to Internet traffic and, according to the authors, the sequential analysis is fast, scalable and efficient; \cite{efficientBAG} explored two methods for scalable inference on Bayesian attack graphs; other models like the one described in \cite{HeRD16} and \cite{Sanna_Heard} are fully parallelisable and suitable for platforms designed for Big Data analysis like Hadoop.

\section{Final remarks}

Cyber security research from a mathematical and statistical point of view is challenging due to the inherent complexity of the problems and the nature of the data. We believe that in order to be well-prepared against the current cyber threats, Bayesian statistics offers a wide range of flexible models that might be the key for a deeper understanding of the generative process at the basis of malicious attacks and, at the same time, for us to have predictive models able to handle large volumes of time-evolving data. That is why in this review we have presented the statistical approach to cyber security anomaly detection methods, making particular emphasis on Bayesian models. However, as remarked in Section 5, the methodologies described in this review are far from being exhaustive. In a highly connected world with cyber threats being more dangerous than ever there is a need for a thorough understanding on the computer networks' behaviour. That is why as the interest in cyber security keeps increasing we are able to find (in a frequent basis) new models that work directly along the line of some of the ones we have presented here, Moreover, alternative approaches to the ones described in this review have been considered and proved useful for both network modelling and anomaly detection.

We would also like to point out that, although there has been an actual increase in cyber security research from a Bayesian point of view, to the best of our knowledge, there are some areas that have not been as widely explored as others. Most of the work we have encountered corresponds to either volume-traffic or network anomaly detection. Malware related problems, like detection and classification, are still open areas of research that need to be deeply developed.

As a final comment, we would like the reader to note that, although it was not mentioned directly in each of the sections of the review paper, anomaly detection models for cyber security research require the analysis of high-volumes of data. No matter if it is for volume-traffic analysis, network modelling or malware detection and classification, all of them require handling and learning from datasets that are usually very large. This definitely plays a vital role in cyber security research, since we have always to keep in mind that whilst developing statistical models for this kind of problems, there is a need for algorithms able to scale well, to be parallelised and, in preference, able to perform in a sequential procedure as new data is observed.

\section*{Acknowledgements}
We would like to thank Professor Nick Heard for his useful suggestions and insightful feedback on the first version of this paper. His comments allowed us to provide a better organised and clearer review. 	We would also like to thank the Associate Editor and the anonymous reviewer for their thoughtful suggestions and comments to improving the review.

\bibliographystyle{hapalike}

redpajama_set_name: RedPajamaArXiv