Detection of epidemic outbreaks with Persistent Causal-Chain Dynamic Bayesian Networks

[0001] 1. Field of the Invention

[0002] The present invention relates to detection of epidemic outbreaks in a large population of hosts. In particular, the present invention relates to detection of an epidemic outbreak by making use of a particular type of Dynamic Bayesian Network which is defined herein as a Persistent Causal-Chain Dynamic Bayesian Network (PCCDBN).

[0003] 2. Description of the Prior Art

[0004] A Bayesian Network is a type of Directed Acyclic Graph which is used to show causal relationships between random variables. A Directed Graph contains a set of nodes and a set of arrows describing a path from node to node. A node X is said to be a parent of another node Y if there is an arrow from X to Y. A Directed Acyclic Graph is a special type of Directed Graph in which there is no directed path from a node back to itself.

[0005] In a Bayesian Network, the nodes represent random variables and the arrows represent the fact that the child node may be conditionally dependent on the parent. The random variables may take on a discrete value such as true or false or a continuous value such as one of the real numbers.

[0006] In a Bayesian Network, for any node X the conditional probability that the random variable X takes on any particular value, given the value of all of X's parents is specified. Given these conditional probabilities, it is then possible to calculate the probability that an event occurred given other events occurring. This process is known as inference.

[0007] In a Bayesian Network, nodes are either hidden or observed. A hidden node, as opposed to an observed node, is a node whose value is not known. A node may be explicitly represented in the model but still hidden due to a lack of observable information. A random variable may also be hidden due to a lack of known conditional independencies, and thus not be represented by a node in the Bayesian network.

[0008] A Dynamic Bayesian Network (DBN) extends the concept of the Bayesian Network into the time dimension. In a Dynamic Bayesian Network, the Bayesian Network is repeated throughout time. Just as arrows create causal connections between nodes in a Bayesian Network, allows are used in a Dynamic Bayesian Network to causally link the network at one time instance to the network at the next time instance.

[0009] When a Dynamic Bayesian Network uses discrete random variables, it may be computationally intractable to solve for the conditional probabilities among nodes. Standard approaches using exact inference models require exponential time due to the large number of cross-temporal dependencies that exist between nodes. Such an exponential approach is intractable and cannot be solved in real time. Other inference models have been developed to solve the network in less than exponential time, however these approaches introduce approximations and cannot give exact solutions.

[0010] Thus, it would be beneficial to have an exact inference model that could solve a type of discrete Dynamic Bayesian Network in less than exponential time by taking advantage of certain properties of the network.

[0021] In the following description, various aspects of the present invention will be described. For purposes of explanation, specific configurations and details are set forth in order to provide a thorough understanding of the present invention. However, it will be apparent to one skilled in the art that the present invention may be practiced without the specific details presented herein. Furthermore, well-known features may be omitted or simplified in order not to obscure the present invention. Various examples are given throughout this description. These are merely descriptions of specific embodiments of the invention. The scope of the invention is not limited to the examples given.

[0022] Unless specifically stated otherwise, as apparent from the following discussions, it is appreciated that throughout the specification, discussions utilizing terms such as “processing,” “computing,” “calculating,” “determining,” or the like, refer to the action and/or processes of a processor, computer or computing system, or similar electronic computing device, that manipulates and/or transforms data represented as physical, such as electronic, quantities within the computing system's registers and/or memories into other data similarly represented as physical quantities within the computing system's memories, registers or other such information storage, transmission or display devices.

[0023] The processes and displays presented herein are not inherently related to any particular computer, communication device or other apparatus. The desired structure for a variety of these systems will appear from the description below. In addition, embodiments of the present invention are not described with reference to any particular programming language, machine code, etc. It will be appreciated that a variety of programming languages, machine codes, etc. may be used to implement the teachings of the invention as described herein.

[0024] Embodiments of the present invention may include apparatus for performing the operations herein. This apparatus may be specially constructed for the desired purposes, or it may comprise a general purpose computer selectively activated or reconfigured by a computer program stored in the computer. Such a computer program may be stored in a computer readable storage medium, such as, but is not limited to, any type of disk including floppy disks, optical disks, magnetic-optical disks, read-only memories, compact disc read-only memories, random access memories, electrically programmable read-only memories, electrically erasable and programmable read only memories, magnetic or optical cards, or any other type of media suitable for storing electronic instructions, and capable of being coupled to a computer system bus.

[0025] An embodiment of a Dynamic Bayesian Network with a time horizon of length M is shown in FIG. 1. In an embodiment of the invention this network may represent a single host under attack. A host may be any computing device that provides a service. Thus, a host may be a personal computer, a file server, a cell phone, a data base system, a web server or the like. An attack may be a malicious attempt to take over or interrupt the normal operation of the host. Thus, an attack may be a virus, a worm, a Trojan horse, a denial of service or the like. In another embodiment of the invention, the network may represent a catastrophic system failure such as a power outage, an assembly line stoppage, or the like.

[0026] In the DBN of FIG. 1, the root cause node A_i may represent a binary random variable that indicates whether an attack has taken place against the host at time instance i. The node I_i may represent a binary random variable that indicates whether the host has been infected at time instance i. The node O_i may represent a binary random variable that indicates the presence or absence of some observable symptom in the host at time instance i. Within a time instance i, an arrow connects A_i to I_i and an allow connects I_i to O_i. This is known as an instantaneous causal chain. The term “instantaneous causal chain” is defined below. Thus, there may be a chain of conditional probabilities linking an attack occurring on the host, an infection occurring on the host given the host is attacked and an observation occurring of a symptom of the host given the host is infected. Such a network may be used to discover whether a host is under attack or whether a host is infected using observations of symptoms within a host. In other words, observable node O_i may be used to determine information about hidden nodes I_i and A_i.

[0027] Between time instances i and i+1, an allow connects A_i to A_i+1 and an arrow connects I_i to I_i+1. It is a reasonable assumption that a host that is under attack will remain under attack. Similarly, it is a reasonable assumption that a host that is infected will remain infected. Thus, for all i, if A_i is true then A_i+1 is also true and if I_i is true then I_i+1 is also true. In other words, the conditional probability of A_i+1 occurring given that A_i has occurred or I_i+1 occurring given that I_i has occurred is 100%. This quality is known as persistence. The term “persistence” is defined below. If a Dynamic Bayesian Network has instantaneous causal chains and exhibits the quality of persistence it is known as a Persistent Causal Chain Dynamic Bayesian Network (PCCDBN).

[0028] To more rigorously define a Persistent Causal Chain Bayesian Network, an embodiment of a more generalized Dynamic Bayesian Network with a causal chain of length N and a time horizon of length M is shown in FIG. 2. A Dynamic Bayesian Network is a PCCDBN if it meets the following criteria:

• • [0029] 1) Each time slice j contains an instantaneous causal chain:
    X₁^j→X₂^j→ . . . →X_N^j→O_j
    • [0030] where X_i^j is a binary hidden random variable for i, j and O_j is a single binary observed random variable.
  • [0031] 2) Each state variable X_i^j exhibits persistence in its on state (the on state is defined as “1” and the off state is defined as “0”):
    P(X_i^j=1|X_i^j−1=1,X_i−1^j=x)=1
    • [0032] For all x ε {0,1} and for 2<j≦M,2<i≦N, and
      P(X₁^j=1|X₁^j−1=1)=1 for all j>1
  • [0033] 3) The conditional distribution of each variable given its parent set is independent of time:
    P(X_i^j=1|X_i^j−1=0,X_i−1^j=1)=P_i,
    P(X_i^j=1|X_i^j−1=0,X_i−1^j=0)= P_i,
    • [0034] with boundary conditions
      P(X₁¹=1)=P(X₁^j=1|X₁^j−1=0)=P₁,
      P(X_i¹=1|X_i−1¹=1)=P_i,
      P(X_i¹=1|X_i−1¹=0)= P_i.

[0035] In Dynamic Bayesian Networks such as the network of FIG. 1, the standard method for performing inference is to group the hidden variables into a single state variable S={A,I}. Doing such converts the DBN into a standard Hidden Markov Model (HMM). Each variable A and I have two possible states 1 and 0 which may also be represented as true and false or on and off. Thus, each variable A and I have a state space n of size two. Such an HMM has a state space of size n_S=n_A×n_I=2×2=4. Because the variables are binary and thus have two possible states it is mote convenient to define n_S as equal to 2^(number of state variables in S).

[0036] An embodiment of a Hidden Markov Model representation of a Dynamic Bayesian Network is shown in FIG. 3. Inference in this model may be performed in O(n_S·M) using belief propagation and/or the Viterbi algorithm. In order to convert the DBN of FIG. 2 into a Hidden Markov Model, state variable S may be defined as {X₁, X₂, . . . , X_N}. Thus, n_S would have a state space of 2^N. As stated above, the size of the state space of n_S grows exponentially with respect to the number of nodes in the causal chain, N. Thus, when a Hidden Markov Model of the network is used, the inference calculation grows exponentially as well and can be performed in O(2^N·M). Such exponential growth can quickly become computationally intractable.

[0037] The temporal dependence that is due strictly to persistence may lead to mole efficient inference. This is because persistence leads to a simpler inference model than general dependence would normally require. For example, persistence reduces the number of allowed configurations of {A₁, A₂, . . . , A_M} and {I₁, I₂, . . . , I_M} from 2^M to M+1. Instead of having M variables each with 2 possible states (2^M configurations), persistence allows for only M+1 configurations because of the assumption that once a variable is “on” it stays on. Thus, if M is 3 there are only the following 4 configurations: {0, 0, 0}, {0, 0, 1}, {0, 1, 1}, and {1, 1, 1}. The first configuration represents the variable never turning on, whereas the last configuration represents the variable turning on immediately.

[0038] Perhaps the simplest way to exploit persistence is to rewrite the model in FIG. 3 to an embodiment of a Persistent Causal Chain Dynamic Bayesian Network as shown in FIG. 4. Two new variables are defined in FIG. 4, à and Ĩ, denoting the time at which A and I become true, respectively. Thus, the effective size of the “state space” in this model is O(M) (where “state space” refers to the number of states of variables à and Ĩ). Standard exact inference models such as belief propagation may be performed on this model in O(M²·N). Reorganizing the model in this way removes the exponential dependence on N, but the calculation is now quadratic in M. In typical DBN applications it may be beneficial to let M grow indefinitely. A quadratic dependence on M may make the calculations quickly become intractable which may not be an ideal solution.

[0039] There may be more structure in the inference calculation above than is captured in the model of FIG. 4. To illustrate this structure, a sketch detailing an embodiment of how to calculate P(O) for the model of FIG. 4 is shown in FIG. 5. In calculating P(O), it may be necessary to marginalize out all the combinations of hidden states, so this calculation can serve as a crucible of sorts for any inference queries that one would want to ask about the network.

[0040] FIG. 5 presents this calculation in a notation that mirrors the notation that will be used later in a more generalized proof. Here P_I is used to denote the probability off I firing (turning on, being set to true, changing to 1) given that A has fired, i e. the true-positive probability. P_I is used to denote the probability of I firing given that A has not fired, i.e. the false-positive probability. Both probabilities assume that I has not filed previously. For conciseness of notation, P_Ī and P_Ī are also defined to be the complements of P_I and P_I, respectively. Similarly, P_A and P_A are used to denote the probability of A filing or not, respectively. Because A does not have any parent nodes, there is no need for P_Ā and P_Ā notation. In the general notation below, P_ī, etc. will be used to denote the conditional probabilities for variable X_i given variable X_i−1.

[0041] The possible configurations of A's and I's through time are indexed by the variables j_A and j_I, respectively. Because of persistence, there are only O(M) values of each of these in general. In the present example j_A ε {0,1,2,3} corresponds to {AAA,ĀAA, AAA, AAA}, respectively, and similarly for j_I. For this notation, AAA means that A fires immediately (configuration 0) and AAA means that A does not fire (configuration 3). These variables have a natural semantics: j_I denotes the latest time slice in which I fails to fire. Finally,

[0042] $\sum _O^{j_I}$
denotes the probability of all observed nodes given that the I variables are in configuration j_I.

[0043] $\sum _O^{j_l}=\prod _{i=1}^{M}P\left(O_i=o_i|I_i=\delta \left(i>j_l\right)\right),\mathrm{where}\delta \left(x\right)=1\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}x\mathrm{is}\mathrm{true}.$

[0044] An important notational point is evident here: in FIG. 5 many P_Ī, etc., variables are raised to a power. These powers are easily confused with superscripts that are also frequently used in this application. Thus, care has been taken to ensure that only P's will be raised to powers and no P's will have superscript indices.

[0045] To better clarify the notation, row j_A equals 1 is explained below. As explained above, j_A equals 1 denotes the configuration in which A does not fire in time instance 0 and instead files (and stays true) at time instance 1. For every configuration of j_A there are four configurations of j_I. Namely, I fires at time instance 0, 1, 2 or not at all. As with A firing, if I fires at time instance x, it does not fire before time instance x and remains true for all time instances after x. The probability for A not firing at time instance 0 and then firing at time instance 1 is equal to P_ĀP_A. Thus, this portion of P(O) in which j_A equals 1 will have P_ĀP_A as a multiplicative factor for all configurations of j_I within the row. For the configuration ĀAAIII, I fires before A, thus this portion of P(O) is equal to the multiplicative factor P_ĀP_A multiplied by P_I (probability of I firing before A) times

[0046] $\sum _O^0$
(probability of all observed nodes given that I is in configuration 0; I fires immediately in time instance 0). This process may be repeated for all j_A and j_I using similar reasoning.

[0047] FIG. 6 shows an embodiment of how linear-time inference can be performed on the network represented in FIG. 4. All the terms necessary for j_A=0 may be calculated in O(M) time. Once these terms are calculated, all the remaining terms may be calculated by making M constant time moves. This may be accomplished by partitioning all the terms for each value of j_A into the σ_I^jA, σ_I^jA, and {circumflex over (σ)}_I^jA terms shown in FIG. 5. For a given configuration, a term may be partitioned into σ_I^jA if I fires at the same time or after A fires. A term may be partitioned into σ_I^jA if I fires before A files. Lastly, a term may be partitioned into {circumflex over (σ)}_I^jA if I never fires. Recursion relations may then be derived relating σ_I^j to σ_I^j+1, σ_I to σ_I^j+1, and {circumflex over (σ)}_I^j to {circumflex over (σ)}_I^j+1. The details of these recursion relations are given below. While calculating the various σ_I^jA, σ_I^jA, and {circumflex over (σ)}_I^jA terms the

[0048] $\sum _I^{j_A}$
terms may also be calculated since they are just the sum of the various σ_I^jA, σ_I^jA, and {circumflex over (σ)}_I^jA terms. If there were more than two variables in the causal chain (for example, in the network of FIG. 2), then this calculation may be recursed by associating

[0049] $\sum _I^{j_A}$
with

[0050] $\sum _O^{j_I}$
and P_I with P_A, etc., N times to give a total complexity of O(N M).

[0051] It is useful to return to the network of FIG. 2 in order to prove this result more generally. For uniformity of notation, P_I≡P_I. M denotes the number of time steps being considered by the model. There is assumed to be an instantaneous causal chain of length N hidden variables, each terminating with a single observation O. The phrase “X_i^j has filed” is used to indicate that X_i^j=1.

[0052] A common shorthand notation X_i:j^j:k is used to denote the set of variables:

[0053] $X_{i\text{:}l}^{j\text{:}k}=\left\{\begin{array}{cccc}{X}_i^j& X_i^{j+1}& \cdots & X_i^k\\ X_{i+1}^j& X_{i+1}^{j+1}& \cdots & X_{i+1}^k\\ ⋮& ⋮& ⋮& ⋮\\ X_l^j& X_l^{j+1}& \cdots & X_l^k\end{array}\right\}$
The notation X_i may be used to refer to an arbitrary member of X_i^1:M, or X^j may be used to refer to an arbitrary member of X_1:N^j.

[0054] One version of the problem of calculating the value of one node from other nodes (known as the filtering problem) is to calculate the probability distribution of the unobserved state variables: P(X_1:N^1:M|O_1:M=o_1:M).

[0055] This quantity may be derived from the evidence (observed nodes) P(O_1:M=o_1:M) in a straightforward manner using Bayes' Rule and the Bayesian network joint factorization. Thus, calculating the evidence will be concentrated upon. Notational shorthand P(O_1:M) may be used to denote P(O_1:M=o_1:M) and similar abbreviations may be used for other probability statements.

[0056] The evidence can be expanded out as follows:

[0057] $\begin{array}{cc}P\left(O_{1\text{:}M}\right)=\sum _{X_{1\text{:}N}^{1\text{:}M}}P\left(O_{1\text{:}M}|X_{1\text{:}N}^{1\text{:}M}\right)·P\left(X_{1\text{:}N}^{1\text{:}M}\right)& \left(\mathrm{EQ}1\right)\end{array}$

[0058] The persistence assumption assigns zero probability to any configuration of X_i^1:M such that X_i^j=0 when X_i^j′=1 for j>j′. Thus, there are only M+1 non-zero probability configurations of X_i^1:M for each i, as illustrated above.

[0059] Let {j_i:0≦j_i≦M} index the sequence of X_i^1:M in which variable X_i^j, is the last (highest-time) variable to be in the off state, unless j_i=0 in which case it indexes the sequence in which all X_i are in the on state. As an example, if M=3, then j_i={0, 1, 2, 3} indexes the states {111, 011, 001, 000}, respectively. All configurations not indexed by j_i have zero probability due to the persistence assumption (for example state 010). To simplify notation, j_i is used to denote the event that X_i^1:M is in the sequence indexed by j_i. Equation 1 can be decomposed according to the network structure as follows:

[0060] $\begin{array}{cc}P\left(O_{1\text{:}M}\right)=\sum _{j_1=0}^{M}P\left(j_1\right)\sum _{j_2=0}^{M}P\left(j_2|j_1\right)·\sum _{j_3=0}^{M}P\left(j_3|j_2\right)\dots \sum _{j_N=0}^{M}P\left(j_N|j_{N-1}\right)P\left(O_{1\text{:}M}|j_N\right)& \left(\mathrm{EQ}2\right)\end{array}$

[0061] Recall that P_i is the probability that variable X_i will fire for the first time given that its causal parent has fired, and P_i is the probability that X_i will fire for the first time given that its causal parent has not fired:
P_i≡P(X_i^j=1|X_i^j−1=0,X_i−1^j=1),  (EQ 3)
and
P_i≡P(X_i^j=1|X_i^j−1=0,X_i−1^j=0).  (EQ 4)

[0062] Let P_ī and P_ī denote the complements of 1−P_i and 1− P_i, respectively.

[0063] Σ_k^I. can be defined recursively to denote the j_k sum from Equation 2 such that j_k−1=L:

[0064] $\begin{array}{cc}\sum _k^L\equiv \sum _{j_k}P\left(j_k|j_{k-1}=L\right)·\sum _{k+1}^{j_k}& \left(\mathrm{EQ}5\right)\end{array}$
with boundary condition

[0065] $\sum _{N+1}^L\equiv P\left(O_{1\text{:}M}|j_N=L\right)$
The simplified example of FIG. 5 shows

[0066] $\sum _I^{j_A}$
for all j_A and with

[0067] $\begin{array}{cc}\sum _{N+1}^L\equiv P\left(O_{1\text{:}M}|j_N=L\right)& \end{array}$

[0068] Using this notation, Equation 2 can be rewritten as:

[0069] $\begin{array}{cc}\sum _{N+1}^{j_I}=\sum _O^{j_I}.& \left(\mathrm{EQ}6\right)\end{array}$

[0070] To clarify this notation, the above equations will be explained with reference to FIG. 5. In this case, N is equal to 2 since there are two hidden nodes A and I. M is equal to 3 and thus there are 4 configurations of each hidden node. Plugging into Equation 2,

[0071] $P\left(O_{1\text{:}3}\right)=\sum _{j_{1=0}}^{3}P\left(j_1\right)\sum _{j_{2=0}}^{3}P\left(j_2|j_1\right)·P\left(O_{1\text{:}3}|j_2\right).$
As explained above, j_i indicates the state in which X_i^1:M is in the sequence indexed by j_i. In FIG. 5 j_I is associated with j_A which indexes A and j₂ is associated with j₁ which indexes I. Thus, the equation can be rewritten as

[0072] $P\left(O_{1\text{:}3}\right)=\sum _{j_{A=0}}^{3}P\left(j_A\right)\sum _{j_{I=0}}^{3}P\left(j_I|j_A\right)P\left(O_{1\text{:}3}|j_I\right).$
Considering only j_A equals 0:
P(j_A=0)*[P(j₁=0|j_A=0)*P(O_1:3|j₁=0)+P(j₁=1|j_A=0)*P(O_1:3|j₁=1)+P(j₁=2|j_A=0)*P(O_1:3|j₁=2)+P(j₁=3|j_A=0)*P(O_1:3|j₁=3)]

[0073] If P_Ī, P_Ī, P_I, P_I, P_A, and P_A are plugged into the equation above, the resultant equation will match the first row (j_A=0) in FIG. 5.

[0074] Each Σ_k^L can be partitioned as follows:
Σ_k^L= σ_k^L+σ_k^L+{circumflex over (σ)}_k^L  (EQ. 7)
σ_k^L contains all the terms in the sum such that X_k fires before X_k−1:

[0075] $\begin{array}{cc}{\stackrel{\_}{\sigma }}_k^L=\sum _{j_k<L}{\stackrel{\_}{P}}_{\stackrel{\_}{k}}^{L-j_k-1}{\stackrel{\_}{P}}_k\sum _{k+1}^{j_k},& \left(\mathrm{EQ}8\right)\end{array}$
σ_k^L contains all the terms in the sum such that X_k fires at the same time as or after X_k−1:

[0076] $\begin{array}{cc}{\sigma }_k^L=\sum _{L\le j_k<M}{\stackrel{\_}{P}}_{\stackrel{\_}{k}}^L{P}_{\stackrel{\_}{k}}^{j_k-L}{P}_k·\sum _{k+1}^{j_k},& \left(\mathrm{EQ}9\right)\end{array}$
and {circumflex over (σ)}_k^L contains all the terms in the sum such that X_k never fires:

[0077] $\begin{array}{cc}{\hat{\sigma }}_k^L=\sum _{L\le j_k<M}{\stackrel{\_}{P}}_{\stackrel{\_}{k}}^L{P}_{\stackrel{\_}{k}}^{M-L}·\sum _{k+1}^M& \left(\mathrm{EQ}10\right)\end{array}$

[0078] This partitioning matches the above partitioning in FIG. 5 into σ_I^jA, σ_I^jA, and {circumflex over (σ)}_I^jA, and the terms here have similar definitions to the terms above.

[0079] In order to calculate Equation 6 in time O(M N), it is necessary to pre-compute σ₂^L, σ₂^L, and {circumflex over (σ)}₂^L for all values of L in O(M N).

[0080] As a boundary condition for the recursion, assume that

[0081] $\sum _{N+1}^i$
has been calculated for all 0≦i≦M. It will be shown below how this is accomplished. It is also necessary to pre-calculate and cache P _lⁱ for 0≦l≦N and 0≦i≦M, which can be done recursively in O(M N) time and space.

[0082] Given these initial quantities, one can calculate σ_Nⁱ for 0≦i≦M in O(M) time using the following recursion:

[0083] $\begin{array}{cc}{\stackrel{\_}{\sigma }}_l^{i+1}={\stackrel{\_}{\sigma }}_l^i+{\stackrel{\_}{P}}_{\stackrel{\_}{l}}^i-{\stackrel{\_}{P}}_l-\sum _{l+1}^i& \left(\mathrm{EQ}11\right)\end{array}$
with boundary condition σ_N⁰=0.

[0084] One can calculate σ_Nⁱ for 0≦i≦M in O(M) time using the following recursion:

[0085] $\begin{array}{cc}{\sigma }_l^{i-1}=\frac{{\sigma }_l^i}{{\stackrel{\_}{P}}_{\stackrel{\_}{l}}}{P}_{\stackrel{\_}{l}}+{\stackrel{\_}{P}}_{\stackrel{\_}{l}}^{i-1}-{\stackrel{\_}{P}}_l-\sum _{l+1}^{i-1}& \left(\mathrm{EQ}12\right)\end{array}$
with boundary condition σ_N^M=0.

[0086] One can calculate {circumflex over (σ)}_Nⁱ for 0≦i≦M in O(M) time using the following recursion:

[0087] $\begin{array}{cc}{\hat{\sigma }}_l^{i-1}=\frac{{\hat{\sigma }}_l^i}{{\stackrel{\_}{P}}_{\stackrel{\_}{l}}}{P}_{\stackrel{\_}{l}}& \left(\mathrm{EQ}13\right)\end{array}$
with boundary condition

[0088] ${\hat{\sigma }}_N^M={\stackrel{\_}{P}}_{\stackrel{\_}{N}}^M-\sum _{N+1}^M.$

[0089] Once σ_Nⁱ, σ_Nⁱ, and {circumflex over (σ)}_Nⁱ are calculated, one can calculate all Σ_Nⁱ for 0≦i≦M in O(M) time using Equation 7. Once Σ_N^0:M is calculated, Equation 5 can be used to calculate Σ_N−1^0:M in time O(M), and so on N times to get all values of Σ_1:N^0:M. Thus, the entire calculation takes O(M N) time.

[0090] In order to calculate Σ_k^L in O(M) time it is necessary to compute Σ_N+1ⁱ (the probability of the observations for a given configuration i of X_N^1:M) for all 0≦i≦M in O(M) time as well. Recall that

[0091] $\sum _{N+1}^l\equiv P\left(O_{1\text{:}M}❘j_N=i\right).$
Since the parents of each O_j is given, for each i, this calculation is simply the product of the observations:

[0092] $\begin{array}{cc}P\left(O_{1\text{:}M}❘j_N=i\right)=\prod _{k=1}^{M}P\left(X_N^k,j_N=i\right)& \left(\mathrm{EQ}14\right)\end{array}$

[0093] Using existing notation, the following are defined:
P_N+1=P(O_k|X_N^k=1),  (EQ 15)
P_N+1=P(O_k|X_N^k=0),  (EQ 16)
and with P _N+1=1−P_N+1 and P _N+1=1− P_N+1. Equation 14 can then be written as:

[0094] $\begin{array}{cc}P\left(O_{1\text{:}M}❘j_N=i\right)=\prod _{k=1}^M{\varphi }_k^{\delta \left(k\le i\right)}-{\stackrel{\_}{\varphi }}_k^{\delta \left(k>i\right)}\mathrm{where}{\varphi }_k=\{\begin{array}{ccc}{P}_{N+1}& \mathrm{if}& O_i=1\\ P_{\stackrel{\_}{N+1}}& \mathrm{otherwise}& \end{array}\mathrm{and}{\stackrel{\_}{\varphi }}_k=\{\begin{array}{ccc}{\stackrel{\_}{P}}_{N+1}& \mathrm{if}& O_l=1\\ {\stackrel{\_}{P}}_{\stackrel{\_}{N+1}}& \mathrm{otherwise}& \end{array}& \left(\mathrm{EQ}17\right)\end{array}$

[0095] Equation 17 can be calculated for all 0≦i≦M in O(M) time via the recursion relation:

[0096] $\begin{array}{cc}\sum _{N+1}^0=\prod _{k=1}^M{\varphi }_i,\mathrm{and}& \left(\mathrm{EQ}18\right)\\ \sum _{N+1}^{i+1}=\frac{\sum _{N+1}^i}{{\varphi }_i}{\stackrel{\_}{\varphi }}_i& \left(\mathrm{EQ}19\right)\end{array}$

[0097] FIG. 1 shows an embodiment of a PCCDBN that may be used for the detection of outbreaks in a population of N_p hosts (the A nodes are shared among the N_p hosts while each host has its own unique set of I nodes and O nodes). In an embodiment of the invention, the model may be used to detect a worm outbreak in a collection of networked devices. In other embodiments of the invention, a nearly identical model may be used to detect any contagious spread in any population as long as each member of the population is being monitored through time (by contrast, many disease surveillance applications for example sample a random subset of the population over time).

[0098] For detection with this model, at each time slice the quantity may be calculated:
ξ=log P(A_M=1|O_1:M)−log P(A_M=0:O_1:M),  (EQ 20)
which may have a threshold to determine if a global attack on the system is taking place.

[0099] This model differs significantly from prior art approaches. Most, if not all, time-series techniques for this task operate according to the model where some statistic(s) from the population is(are) sampled from the population at regular intervals, either by intentional periodic surveillance or by incidental measurements. Under this schema, the models that come closest to PCCDBN are HMM models. However, this model is novel for this task in that it may maintain a time-series of observations of each individual in the population and may use the model to merge those time-series in a principled manner. This added information qualitatively improves the ability to detect an outbreak.

[0100] Besides being novel in the way it monitors the population as well as being tractable for inference, the model in FIG. 1 may be very efficiently parallelized and distributed. Each local host H_i may maintain its own observation history and calculate its own signal ξ_i. The aggregate signal for the entire population may be calculated by summing over all local signals (ignoring the prior): ξ=Σ_iξ_i.

[0101] In a distributed setting, this may be implemented over a network by having hosts send messages containing their local ξ_i values. Those ξ_i values may be aggregated with other messages that the host has received to allow for exponential growth in the number of hosts that can participate in a given message, although the messages still can maintain a constant size. This technique may be prone to over-counting; although that may be mitigated by allowing messages to expire after a few hops and by ensuring that no host contributes to a message with the same originating host-time pair more than once.

[0102] The empirical setup for this model includes a system of noisy local detectors (LDs), which are noisy observable detectors sitting on the end-hosts. The LDs scan outgoing packets and count the number of UDP packets to unique destination IPs, itemized by source and destination port. The LDs may also scan memory used, CPU usage, and the like. PCCDBN inference was tested on five weeks of network traffic data for 37 hosts within a corporate enterprise network. In experiments where more than 37 hosts were used, the data was bootstrapped by uniformly-at-random selecting a starting-point for each virtual host at one of the original 37 host's data stream (original host again chosen uniformly-at-random).

[0103] FIG. 7 shows a graph detailing noisy observations of outgoing packets by the LDs aggregated over all machines for all five weeks. The LD was operated at a low threshold of 4 connections per 50 second intervals (shown as the vertical line in the figure). Thus, if the count is 4 or greater in time slice i, O_i is set to 1.

[0104] This setup is designed to catch very slow worms that try to escape detection by sliding well under the background traffic, as seen by the relatively large histogram mass to the right of the LD threshold in FIG. 7. The lines in this plot labeled W1, W2, . . . , W5 show where some recent worm attacks would appear using this heuristic; the point being that this heuristic is very good for detection of fast worms, but for sneaky worms with slow propagation rate, this detector is very noisy, having a false positive rate on the order of thousands per week.

[0105] The simplicity of the LDs made simulating worm behavior on top of real network traffic trivial. A worm was modeled with two parameters: the spread rate S which indicated how many hosts the worm attempted to spread to per unit time, and the address density with which the worm was able to spread. For example, if the worm is randomly choosing addresses from the entire 32-bit address space and a network has 37 hosts, then the address density will be 37/2³². All results shown here use an S of 1 connection per 20 seconds, and set the address density to 1/1000.

[0106] Simulations superimposed artificial worm traffic on the real enterprise traffic by releasing worms, letting the proportion of infected hosts grow to 1 in the population, then repeating 20 times. To test for false positive rates the enterprise traffic with no worm traffic superimposed was used. The results were averaged over all runs.

[0107] Very similar experiments were performed, where a Naive Bayes approximation to PCCDBNs was used for efficiency and was compared to several standard algorithms including a modified HMM and a Cumulative Sum (CUSUM) detector. Both the Naive-DBN and the modified HMM out-performed the rest of the tested models. The HMM model actually performed slightly better than the Naive-DBN model; however, the Naive Bayes DBN approximation performed qualitatively similarly and was much faster to run, so it was used it as the baseline in the experiments presented here.

[0108] The local detector false-positive rate P(O_i=1|I_i=0) was estimated online using a Bayesian update assuming that all traffic was clean. After an outbreak, the parameters were reset to avoid having outbreak data inflate the assumed false positive rates. Each LD was allowed to maintain its own false positive and true positive rates independent of the other LDs. The true-positive rate P(O_i=1|I_i=1) was fixed to a constant 0.95. The same method for estimating parameters was used for the Naive-DBN and the PCCDBN.

[0109] Regarding the state variables, it is assumed that:
P(I_i=1|I_i−1=0, A_i=1)=1/N_p  (EQ 22)
P(I_i=1|I_i−1=0, A_i=0)=0,  (EQ 23)
and that
P(A_i=1|A_i−1=0)=0.001,  (EQ 24)
however, detection may be measured by sweeping through all detection thresholds, so the results are expected to be insensitive to the exact values of the priors if they are non-zero.

[0110] In the experiments here, three separate models are compared: the distributed Naive-DBN approximation and the PCCDBN model in both centralized form (C-PCCDBN), where all hosts report to a single detector every 10 seconds, and a distributed model (D-PCCDBN) where each host sent out 1 message every 10 seconds to a random other host. Messages were aggregated and given a time-to-live of 4 hops.

[0111] The number of hosts was varied over the set N_p ε {50, 100, 200, 400, 800, 1600}. The results ate shown in FIG. 8

[0112] These results are striking for at least two reasons. First, they show that an embodiment of D-PCCDBN possesses qualitatively different detection ability than either of the other two models (and all the models tested previously for that matter). Namely, D-PCCDBN made a detection when roughly a fixed absolute value of hosts was infected, as opposed to when a fixed percentage was infected. The second striking fact from FIG. 8 was that an embodiment of C-PCCDBN did not display the same qualitative behavior that D-PCCDBN displayed, and performed much worse. This is in line with previous results which showed that the distributed embodiment of Naive-DBN slightly out-performed the centralized version. However, the difference there was relatively minor, and the general trend for scaling up was qualitatively similar.

[0113] Further results are shown in FIG. 9. The goal of epidemic outbreak detection is to generate a large enough detection signal before a significant fraction of the population gets infected The graph of FIG. 9 shows an infection spreading through a population of 50 computer hosts using data from the corporation enterprise over 5 weeks. At time 50,000 a worm is introduced in the system and begins to spread. The intermediate spikes in the detection signal during periods of non-infection are due to diurnal cycles of activity, causing the signal to increase during the high-traffic workday periods. As can be seen from the graph, the detection signal spikes significantly after the infection begins.

[0114] Finally, PCCDBNs have been tested on the task of detecting a worm outbreak in a computer network and it has been shown that when the calculation is distributed in a particular manner, the detection results are qualitatively better than other cutting edge techniques, allowing detection (at less than 1 false positive per week) when a fixed absolute number of hosts are infected rather than when a fixed percentage is infected. The qualitative difference does not seem to hold when PCCDBNs are used in a centralized manner, reaffirming the previous findings that distributed detection can actually perform better than centralized detection, despite the fact that the centralized detector has available all information; whereas no single distributed detector has that information.