## 5.3 Network failure

In this section we consider a more serious application of importance sampling. Though still a toy example, where we can find an exact solution, the example illustrates well the type of application where we want to approximate a small probability using a Monte Carlo average. Importance sampling can then increase the probability of the rare event and as a result make the variance of the Monte Carlo average smaller.

We will consider the following network consisting of ten nodes and with some of the nodes connected.

The network could be a computer network with ten computers. The different connections (edges) may “fail” independently with probability \(p\), and we ask the question: what is the probability that node 1 and node 10 are disconnected?

We can answer this question by computing an integral of an indicator function, that is, by computing the sum \[\mu = \sum_{x \in \{0,1\}^{18}} 1_B(x) f_p(x)\] where \(f_p(x)\) is the point probability of \(x\), with \(x\) representing which of the 18 edges in the graph that fail, and \(B\) represents the set of edges where node 1 and node 10 are disconnected. By simulating edge failures we can approximate the sum as a Monte Carlo average.

The network of nodes can be represented as a graph adjacency matrix \(A\) such that \(A_{ij} = 1\) if and only if there is an edge between \(i\) and \(j\) (and \(A_{ij} = 0\) otherwise).

```
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 0 1 0 1 1 0 0 0 0 0
## [2,] 1 0 1 0 0 1 0 0 0 0
## [3,] 0 1 0 0 0 1 1 1 0 1
## [4,] 1 0 0 0 1 0 0 1 0 0
## [5,] 1 0 0 1 0 0 0 1 1 0
## [6,] 0 1 1 0 0 0 1 0 0 1
## [7,] 0 0 1 0 0 1 0 0 0 1
## [8,] 0 0 1 1 1 0 0 0 1 0
## [9,] 0 0 0 0 1 0 0 1 0 1
## [10,] 0 0 1 0 0 1 1 0 1 0
```

To compute the probability that 1 and 10 are disconnected by Monte Carlo integration, we need to sample (sub)graphs by randomly removing some of the edges. This is implemented using the upper triangular part of the (symmetric) adjacency matrix.

```
sim_net <- function(Aup, p) {
ones <- which(Aup == 1)
Aup[ones] <- sample(
c(0, 1),
length(ones),
replace = TRUE,
prob = c(p, 1 - p)
)
Aup
}
```

The core of the implementation above uses the `sample()`

function, which
can sample with replacement from the set \(\{0, 1\}\). The vector `ones`

contains indices of the (upper triangular part of the) adjacency matrix
containing a `1`

, and these positions are replaced by the sampled values before
the matrix is returned.

It is fairly fast to sample even a large number of random graphs this way.

```
## user system elapsed
## 1.396 0.094 1.506
```

The second function we implement checks network connectivity based on the upper triangular part of the adjacency matrix. It relies on the fact that there is a path from node 1 to node 10 consisting of \(k\) edges if and only if \((A^k)_{1,10} > 0\). We see directly that such a path needs to consist of at least \(k = 3\) edges. Also, we don’t need to check paths with more than \(k = 9\) edges as they will contain the one node multiple times and can thus be shortened.

```
discon <- function(Aup) {
A <- Aup + t(Aup)
i <- 3
Apow <- A %*% A %*% A # A%^%3
while(Apow[1, 10] == 0 & i < 9) {
Apow <- Apow %*% A
i <- i + 1
}
Apow[1, 10] == 0 # TRUE if nodes 1 and 10 not connected
}
```

We then obtain the following estimate of the probability of nodes 1 and 10 being disconnected using Monte Carlo integration.

```
seed <- 27092016
set.seed(seed)
n <- 1e5
tmp <- replicate(n, discon(sim_net(Aup, 0.05)))
mu_hat <- mean(tmp)
```

As this is a random approximation, we should report not only the Monte Carlo estimate but also the confidence interval. Since the estimate is an average of 0-1-variables, we can estimate the variance, \(\sigma^2\), of the individual terms using that \(\sigma^2 = \mu (1 - \mu)\). We could just as well have used the empirical variance, which would give almost the same numerical value as \(\hat{\mu} (1 - \hat{\mu})\), but we use the latter estimator to illustrate that any (good) estimator of \(\sigma^2\) can be used when estimating the asymptotic variance.

`## [1] 0.000226 0.000340 0.000454`

The estimated probability is low and only in about 1 of 3000 simulated graphs will node 1 and 10 be disconnected. This suggests that importance sampling can be useful if we sample from a probability distribution with a larger probability of edge failure.

To implement importance sampling we note that the point probabilities (the density w.r.t. counting measure) for sampling the 18 independent 0-1-variables \(x = (x_1, \ldots, x_{18})\) with \(P(X_i = 0) = p\) is

\[f_p(x) = p^{18 - s} (1- p)^{s}\]

where \(s = \sum_{i=1}^{18} x_i\). In the implementation, weights are computed that correspond to using probability \(p_0\) (with density \(g = f_{p_0}\)) instead of \(p\), and the weights are only computed if node 1 and 10 are disconnected.

```
weights <- function(Aup, Aup0, p0, p) {
w <- discon(Aup0)
if (w) {
s <- sum(Aup0)
w <- (p / p0)^18 * (p0 * (1 - p) / (p * (1 - p0)))^s
}
as.numeric(w)
}
```

The implementation uses the formula

\[ w(x) = \frac{f_p(x)}{f_{p_0}(x)} = \frac{p^{18 - s} (1- p)^{s}}{p_0^{18 - s} (1- p_0)^{s}} = \left(\frac{p}{p_0}\right)^{18} \left(\frac{p_0 (1- p)}{p (1- p_0)}\right)^s. \]

The importance sampling estimate of \(\mu\) is then computed.

```
set.seed(seed)
tmp <- replicate(n, weights(Aup, sim_net(Aup, 0.2), 0.2, 0.05))
mu_hat_IS <- mean(tmp)
```

And we obtain the following confidence interval using the empirical variance estimate \(\hat{\sigma}^2\).

`## [1] 0.000262 0.000296 0.000330`

The ratio of estimated variances for the plain Monte Carlo estimate and the importance sampling estimate is

`## [1] 11.22476`

Thus we need around 11 times more samples if using plain Monte Carlo integration when compared to importance sampling to obtain the same precision. A benchmark will show that the extra computing time for importance sampling is small compared to the reduction of variance. It is therefore worth the coding effort if used repeatedly, but not if it is a one-off computation.

The graph is, in fact, small enough for complete enumeration and thus the computation
of an exact solution. There are \(2^{18} = 262,144\) different networks with
any number of the edges failing. To
systematically walk through all possible combinations of edges failing,
we use the function `intToBits`

that converts an integer to its binary
representation for integers from 0 to 262,143. This is a quick and convenient
way of representing all the different fail and non-fail combinations
for the edges.

```
ones <- which(Aup == 1)
Atmp <- Aup
p <- 0.05
prob <- numeric(2^18)
for(i in 0:(2^18 - 1)) {
on <- as.numeric(intToBits(i)[1:18])
Atmp[ones] <- on
if (discon(Atmp)) {
s <- sum(on)
prob[i + 1] <- p^(18 - s) * (1 - p)^s
}
}
```

The probability that nodes 1 and 10 are disconnected can then be computed as
the sum of all the probabilities in `prob`

.

`## [1] 0.000288295`

This number should be compared to the estimates computed above. For a more complete comparison, we have used importance sampling with edge fail probability ranging from 0.1 to 0.4, see Figure 5.6. The results show that a failure probability of 0.2 is close to optimal in terms of giving an importance sampling estimate with minimal variance. For smaller values, the event that 1 and 10 become disconnected is too rare, and for larger values the importance weights become too variable. A choice of 0.2 strikes a good balance.

### 5.3.1 Object oriented implementations

Implementations of algorithms that handle graphs and simulate graphs, as in the Monte Carlo
computations above, can benefit from using an object oriented approach. To
this end we first implement a so-called *constructor*, which is a function
that takes the adjacency matrix and the edge failure probability as arguments and
returns a list with class label `network`

.

```
network <- function(A, p) {
Aup <- A
Aup[lower.tri(Aup)] <- 0
ones <- which((Aup == 1))
structure(
list(
A = A,
Aup = Aup,
ones = ones,
p = p
),
class = "network"
)
}
```

We use the constructor function `network()`

to construct and object of class
`network`

for our specific adjacency matrix.

```
## List of 4
## $ A : num [1:10, 1:10] 0 1 0 1 1 0 0 0 0 0 ...
## $ Aup : num [1:10, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
## $ ones: int [1:18] 11 22 31 41 44 52 53 63 66 73 ...
## $ p : num 0.05
## - attr(*, "class")= chr "network"
```

`## [1] "network"`

The network object contains, in addition to `A`

and `p`

, two precomputed components:
the upper triangular part of the adjacency matrix; and the indices in that matrix
containing a `1`

.

The intention is then to write two methods for the network class. A method `sim()`

that will
simulate a graph where some edges have failed, and a method `failure()`

that will
estimate the probability of node 1 and 10 being disconnected by Monte Carlo
integration. To do so we need to define the two corresponding generic functions.

The method for simulation is then implemented as a function with name
`sim.network`

.

```
sim.network <- function(x) {
Aup <- x$Aup
Aup[x$ones] <- sample(
c(0, 1),
length(x$ones),
replace = TRUE,
prob = c(x$p, 1 - x$p)
)
Aup
}
```

It is implemented using essentially the same implementation as `sim_net()`

except that `Aup`

and `p`

are extracted as components from the object `x`

instead of being arguments, and `ones`

is extracted from `x`

as well instead
of being computed. One could argue that the `sim()`

method should return an object
of class network — that would be natural. However, then we need to call the

constructor with the full adjacency matrix, this will take some time and we do
not want to do that as a default. Thus we simply return the upper triangular
part of the adjacency matrix from `sim()`

.

The `failure()`

method implements plain Monte Carlo integration as well
as importance sampling and returns a vector containing the estimate as
well as the 95% confidence interval. This implementation relies on the
already implemented functions `discon()`

and `weights()`

.

```
failure.network <- function(x, n, p0 = NULL) {
if (is.null(p0)) {
# Plain Monte Carlo simulation
tmp <- replicate(n, discon(sim(x)))
mu_hat <- mean(tmp)
se <- sqrt(mu_hat * (1 - mu_hat) / n)
} else {
# Importance sampling
p <- x$p
x$p <- p0
tmp <- replicate(n, weights(x$Aup, sim(x), p0, p))
se <- sd(tmp) / sqrt(n)
mu_hat <- mean(tmp)
}
value <- mu_hat + 1.96 * se * c(-1, 0, 1)
names(value) <- c("low", "estimate", "high")
value
}
```

We test the implementation against the previously computed results.

```
## low estimate high
## 0.000226 0.000340 0.000454
```

```
## low estimate high
## 0.000262 0.000296 0.000330
```

We find that these are the same numbers as computed above, thus the object oriented implementation concurs with the non-object oriented on this example.

We benchmark the object oriented implementation to measure if there is
any run time loss or improvement due to using objects.
One should expect a small computational overhead due to method dispatching,
that is, the procedure that R uses to look up the appropriate `sim()`

method
for an object of class `network`

. On the other hand, `sim()`

does not recompute
`ones`

every time.

```
## Unit: microseconds
## expr min lq mean median uq max neval
## sim_net(Aup, 0.05) 8.34 9.09 10.7 9.43 10.0 202 10000
## sim(my_net) 9.30 10.13 12.8 10.49 11.1 10153 10000
```

From the benchmark, the object oriented solution using `sim()`

appears to be
a bit slower than `sim_net()`

despite the latter recomputing `ones`

, and
this can be explained by method dispatching taking of the order of 1 microsecond
during these benchmark computations.

Once we have taken an object oriented approach, we can also implement methods
for some standard generic functions, e.g. the `print`

function. As this generic
function already exists, we simply need to implement a method for class `network`

.

```
print.network <- function(x) {
cat("#vertices: ", nrow(x$A), "\n")
cat("#edges:", sum(x$Aup), "\n")
cat("p = ", x$p, "\n")
}
```

```
## #vertices: 10
## #edges: 18
## p = 0.05
```

Our print method now prints out some summary information about the graph instead of just the raw list.

If you want to work more seriously with graphs, it is likely that you want to use an existing R package instead of reimplementing many graph algorithms. One of these packages is igraph, which also illustrates well an object oriented implementation of graph classes in R.

We start by constructing a new graph object from the adjacency matrix.

`## [1] "igraph"`

```
## IGRAPH 0429f5b U--- 10 18 --
## + edges from 0429f5b:
## [1] 1-- 2 1-- 4 1-- 5 2-- 3 2-- 6 3-- 6 3-- 7 3-- 8 3--10 4-- 5 4-- 8 5-- 8
## [13] 5-- 9 6-- 7 6--10 7--10 8-- 9 9--10
```

The igraph package supports a vast number of graph computation, manipulation and visualization tools. We will here illustrate how igraph can be used to plot the graph and how we can implement a simulation method for objects of class igraph.

You can use `plot(net)`

, which will call the plot method for objects of class igraph.
But before doing so, we will specify a layout of the graph.

```
# You can generate a layout ...
net_layout <- layout_(net, nicely())
# ... or you can specify one yourself
net_layout <- matrix(
c(-20, 1,
-4, 3,
-4, 1,
-4, -1,
-4, -3,
4, 3,
4, 1,
4, -1,
4, -3,
20, -1),
ncol = 2, nrow = 10, byrow = TRUE)
```

The layout we have specified makes it easy to recognize the graph.

We then use two functions from the igraph package to implement simulation of
a new graph. We need `gsize()`

, which gives the number of edges in the graph,
and we need `delete_edges()`

, which removes edges from the graph. Otherwise
the simulation is still based on `sample()`

.

```
sim.igraph <- function(x, p) {
deledges <- sample(
c(TRUE, FALSE),
gsize(x), # 'gsize()' returns the number of edges
replace = TRUE,
prob = c(p, 1 - p)
)
delete_edges(x, which(deledges))
}
```

Note that this method is also called `sim()`

, yet there is no conflict here with
the method for objects of class network because the generic `sim()`

function
will delegate the call to the correct method based on the objects class label.

If we combine our new `sim()`

method for igraphs with the plot method, we can
plot a simulated graph.

The implementation using igraph turns out to be a little slower than using
the matrix representation alone as in `sim_net()`

.

```
## user system elapsed
## 4.118 1.993 6.151
```

One could also implement the function for testing if nodes 1 and 10 are
disconnected using the `shortest_paths()`

function, but this is not faster
than the simple matrix multiplications used in `discon()`

either. However, one
should be careful not to draw any general conclusions from this. Our graph is
admittedly a small toy example, and we implemented our solutions largely
to handle this particular graph.