The ranking function takes w(di, cl**) as an input.**
Various forms of ranking functions can be constructed. For example, [Maeno (2009)] studied a simple form in eq.(5) where the function u(x) returns 1 if the real variable x **is positive, or 0 otherwise.**

$$\begin{array}{ll}s(d_{i})&\propto&\sum_{l=0}^{C-1}u(w(d_{i},c_{l}))\\ &=&\sum_{l=0}^{C-1}B(d_{i}\cap c_{l}\neq\phi).\end{array}\tag{5}$$

The i-th most suspicious log is given by dσ(i) **where**
σ(i) is calculated by eq.(6). Suspiciousness s(dσ(i)**) is** always larger than s(dσ(i
′)) for any **i < i**′.

$$\sigma(i)=\arg\max_{m\neq\sigma(n)\ \mbox{for}\ \forall n<i}s(d_{m})\ (1\leq i\leq D).\tag{6}$$

The computational burden of the method remains light as the number of nodes and surveillance logs increases. The method is expected to work generally for clustered networks but moderately even if the network topological and stochastic mechanism to generate the surveillance logs is not understood well. The method works without the knowledge about the hub-and-spoke model; the parametric form with rjk and fj **in Section 3.**
The result, however, can not be very accurate because of the heuristic nature. A statistical inference method which requires heavy computational burden, but outputs more accurate results is presented next.

## 4.2 Statistical Inference Method

The statistical inference method employs the maximal likelihood estimation to infer the topology of the network, and applies an anomaly detection technique to retrieve the suspicious surveillance logs which are not likely to realize without the covert nodes. The maximal likelihood estimation is a basic statistical method used for fitting a statistical model to data and for providing estimates for the model's parameters. The anomaly detection refers to detecting patterns in a given dataset that do not conform to an established normal behavior.

A single symbol θ represent both of the parameters rjk and fj for the nodes in O. θ **is the target**
variable, the value of which needs to be inferred from the surveillance log dataset. The logarithmic likelihood function [Hastie (2001)] is defined by eq.(7). The quantity p({di}|θ) denote the probability where the surveillance log dataset {di} **realizes under a given** θ.

$$L(\mathbf{\theta})=\log(p(\{d_{i}\}|\mathbf{\theta})).$$
L(θ) = log(p({di}|θ)). (7)
The individual surveillance logs are assumed to be independent. eq.(7) becomes eq.(8).

$$\begin{array}{r c l}{{L(\mathbf{\theta})}}&{{=}}&{{\log(\prod_{i=0}^{D-1}p(d_{i}|\mathbf{\theta}))}}\end{array}$$
$$\begin{array}{r l r l}{={}}&{{}\sum_{i=0}^{D-1}\log(p(d_{i}|\mathbf{\theta})).}&{}&{{}}&{{}{\mathrm{(8)}}}\end{array}$$
4
The quantity qi|jk **in eq.(9) is the probability where**
the presence or absence of the node nk **as a responder to** the stimulating node nj **coincides with the surveillance**
log di.

$$q_{i|jk}=\left\{\begin{array}{ll}r_{jk}&\mbox{if$\mathbf{d}_{ik}=1$for given$i$and$j$}\\ 1-r_{jk}&\mbox{otherwise}\end{array}\right..\tag{9}$$

Eq.(9) is equivalent to eq.(10) since the value of dik is either 0 or 1.

$$(10)$$
$$d_{i k})(1-r_{j k}).$$

qi|jk = dikrjk + (1 − dik)(1 − rjk). **(10)**
The probability p({di}|θ**) in eq.(8) is expressed by**
eq.(11).

$$p(d_{i}|\mathbf{\theta})=\sum_{j=0}^{N-1}\mathbf{d}_{ij}f_{j}\prod_{0\leq k<N\ \wedge\ k\neq j}q_{i|jk}.\tag{11}$$

The logarithmic likelihood function takes an explicit formula in eq.(12). The case k = j **in multiplication**
(Qk
) is included since d 2 ik = dik **always holds.**

$$L(\mathbf{\theta})=\sum_{i=0}^{D-1}\log(\sum_{j=0}^{N-1}\mathbf{d}_{ij}f_{j}\prod_{k=0}^{N-1}\{1-\mathbf{d}_{ik}$$ $$+(2\mathbf{d}_{ik}-1)r_{jk}\}).\tag{12}$$
$$(13)$$

The maximal likelihood estimator θˆ **is obtained by**
solving eq.(13). It gives the values of the parameters rjk and fj . A pair of nodes nj and nk **for which** rjk > 0 possesses a link between them.

$${\hat{\boldsymbol{\theta}}}=\arg\operatorname*{max}_{\boldsymbol{\theta}}L({\boldsymbol{\theta}}).$$

L(θ). **(13)**
A simple incremental optimization technique; the hill climbing method (or the method of steepest descent)
is employed to solve eq.(13). Non-deterministic methods such as simulated annealing [Hastie (2001)] can be employed to strengthen the search ability and to avoid sub-optimal solutions. These methods search more optimal parameter values around the present values and update them as in eq.(14) until the values converge.

$$\left\{\begin{array}{l l}{r_{j k}\to r_{j k}+\Delta r_{j k}}&{(0\leq j,k<N).}\\ {f_{j}\to f_{j}+\Delta f_{j}}&{}\end{array}\right.$$
$$(14)$$
$\left(7\right)$. 
The change in the logarithmic likelihood function can be calculated as a product of the derivatives (differential coefficients with regard to r and f**) and the**
amount of the updates in eq.(15). The update ∆rnm and ∆fn **should be in the direction of the steepest ascent**
in the landscape of the logarithmic likelihood function.

$$\Delta L(\mathbf{\theta})=\sum_{n,m=0}^{N-1}\frac{\partial L(\mathbf{\theta})}{\partial r_{n m}}\Delta r_{n m}+\sum_{n=0}^{N-1}\frac{\partial L(\mathbf{\theta})}{\partial f_{n}}\Delta f_{n}.\quad(15)$$
