## Simulations of the fast probabilistic consensus protocol (FPC)

Wind extinguishes a candle and energizes fire. Likewise with randomness, uncertainty, chaos: you want to use them, not hide from them. You want to be the fire and wish for the wind.

Nassim Nicholas Taleb

#### TL;DR

The fast probabilistic consensus (FPC), introduced recently by Serguei Popov (paper) terminates in most of the times in a minimal number of steps and is much faster than the theoretical upper bounds given by Popov. It is claimed that an additional layer of randomness allows the protocol to still work in critical situations where other protocols fail. These simulations show that the FPC is indeed very robust under various attack scenarios ->.

## Introduction and set-up

We perform a first simulation-based study on the fast probabilistic consensus proposed by Serguei Popov (paper). Since the protocol is very safe and fast in non-critical situations, this study is mainly about attack strategies that try to reach disagreement of the honest nodes and considers mostly critical situations. The post is not intended to be a thorough scientific study of the FPC; it justs collects various observations in order to illustrate its functionalities.

### Question

How can a distributed network find consensus on the value of a bit?

## Donsker’s invariance principle

Donsker’s invariance principle is a classical extension of the central limit theorem. It says that the path of a discrete random walk converges, if properly rescaled, in distribution to a standard Brownian motion. More details can be found in your favorite text book or on wikipedia. Here just comes an animation:

library(dplyr)
library(gganimate)
library(transformr)
set.seed(42)
n <- 1000  # maximal steps of random walk
n_seq <- seq(from=10, 100, by = 10)  # starting phase
n_seq <- c(n_seq, seq(from=100, n, by = 100))
number_iteration <-length(n_seq)
X <- sample(c(-1,1), n, replace=TRUE, prob=c(0.5, 0.5))
S <- cumsum(X)  # random walk S_n

W <- numeric(number_iteration*n)
id <- numeric(number_iteration*n)
time <- numeric(number_iteration*n)
for (i in 1:number_iteration){
for (j in (1:n)){
index <- floor((j-1)/n * n_seq[i])
if (index==0) {
W[j+(i-1)*n] <-((((j-1)/n * n_seq[i]) -(index)) *( S[index+1]))/ sqrt(n_seq[i])
}
else {
W[j+(i-1)*n] <- (S[index] + (((j-1)/n * n_seq[i]) -(index)) *( S[index+1]-S[index]))/ sqrt(n_seq[i])
}
}
id[(1+n*(i-1)):(n*i)] <- rep(n_seq[i],n)
time[(1+n*(i-1)):(n*i)] <- (1:n)/n
}

rw <- data.frame(x = time,
y = W,
id =id)

p <- ggplot(rw, aes(x=x, y=y)) +
geom_line()+
ylab(expression(W^{(n)}~ (x)))+
transition_states(
id,
transition_length = 2,
state_length = 1) +
exit_shrink() +
ease_aes('sine-in-out')+
labs(title = "Donsker's invariance principle. n = {closest_state}")

animate(p, nframes=number_iteration*2)

## Brownian motion via Hilbert basis

Brownian motion is a central object of probability theory. One perspective to look at its construction uses Hilbert spaces.
Let be a Hilbert basis equipped with standard scalar product . Hence, every can be written in a unique way as

where and . Let be an i.i.d. sequeence of standard gaussian random variables and consider

Using the independence of we obtain

Brownian motion is then given by (a continuous modification of) , . Hence, “all” we need to do is to find a Hilbert basis. The most known is the following.

#### Haar basis

Here we refer to http://datatreker.com/haar-basis for an introduction.

haar_mother <- function(x){
(x >0 & x <= 0.5) - (x >0.5 & x <= 1)
}
haar <- function(x,j,k){
(2^(j/2) * haar_mother(2^j*x-k))
}
set.seed(42)

j_max <- 10
n_max <- 11  # maximal resolution
delta <- 2^{-n_max}
x <- (0: 2^n_max)/2^n_max
xi_0 <- rnorm(1) # xi corresponding to constant function
xi <- list()     # list of random variables xi
for (j in 0:j_max){
xi[[j+1]] <- rep(0, 2^j)
for (k  in 0:(2^j-1))
{
xi[[j+1]][k+1] <-rnorm(1)
}
}

# data.frame containing approximation of the Brownian motion
df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())
for (i in x){
y <- c(1,rep(1, i*2^n_max), rep(0, 2^n_max-i*2^n_max))
alpha <- list() # wavelet coefficients for y
bm <- sum(y)/(length(y)+1) * xi_0
for (j in 0:j_max){
alpha[[j+1]] <- rep(0, 2^j)
for (k  in 0:(2^j-1))
{
alpha[[j+1]][k+1] <-sum(haar(x,j,k)*y)*delta
bm <- bm + alpha[[j+1]][k+1]*xi[[j+1]][k+1]
}
df_new <- data.frame(x=i,
y=bm,
id=j)
df <- bind_rows(df, df_new)
}
}

#### Animation using gganimate

ggplot(df, aes(x=x, y=y)) +
geom_line()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = "Construction of Brownian motion using Haar basis. Step: {closest_state}") +
ease_aes('sine-in-out')

The result looks the same as with Levy’s construction which is not astonishing. The implementation above is not optimal and is much slower than the previous one.

It might be worth to consider a different Hilbert basis.

#### Trigonometric Hilbert basis

The trigonometric basis is defined as

together with .

set.seed(42)
N_max <-  200 # number of trigonometric basis functions used
n_max <- 11  # maximal resolution
delta <- 2^{-n_max}
x <- (0: 2^n_max)/2^n_max
xi <- rnorm(N_max+1)     # list of random variables xi

# data.frame containing approximation of the Brownian motion
df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())

for (i in x){
y <- c(1,rep(1, i*2^n_max), rep(0, 2^n_max-i*2^n_max))
alpha <- numeric() # wavelet coefficients for y
alpha[1] <- mean(y)
bm <-  alpha[1] * xi[1]
for (n in 1:N_max){
alpha[n] <-sum(sqrt(2)*cos(x*pi*n)*y)*delta
bm <- bm + alpha[n]*xi[n]
df_new <- data.frame(x=i,
y=bm,
id=n)
df <- bind_rows(df, df_new)}
}

#### Animation using gganimate

p <- ggplot(df, aes(x=x, y=y)) +
geom_line()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = "Construction of Brownian motion using trigonometric basis. Step: {closest_state}") +
ease_aes('sine-in-out')
animate(p, nframes=1000)

The above convergence is very slow. The irregularity properties of the Brownian motion are not yet very visible.

## Haar basis

Hilbert spaces play a prominent role in various fields of mathematics. An orthonormal basis of such a space is called a Hilbert basis. The purpose of this blog is to illustrate a very clasic and basic Hilbert basis – the Haar basis.

Let be a Hilbert basis of equipped with the standard scalar product . Hence, every can be written in a unique way as

where and .

A classic Hilbert basis consists of Haar functions that are supported on . They are defined using the Haar wavelet:

The Haar basis consists then of rescaled (by ) versions of shifted by ,

together with the constant function . (Note: the constant function has to be added since we consider the interval and not ). In R this looks like:

haar_mother <- function(x){
(x >0 & x <= 0.5) - (x >0.5 & x <= 1)
}
haar <- function(x,j,k){
2^(j/2) * haar_mother(2^j*x-k)
}

#### Animation of the Haar basis

j_max <- 3    # maximal depth
n_max <- 10   # resolution of grid
df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())

x <- (1: 2^n_max)/2^n_max
id<-1
for (j in 0:j_max){
for (k  in 0:(2^j-1))
{
df_new<- data.frame(x=x,
y=haar(x,j,k),
id=id)
df <- bind_rows(df, df_new)
id <- id+1
}
}

ggplot(df, aes(x=x, y=y)) +
geom_step()+
transition_states(
id,
transition_length = 2,
state_length = 5
) +
labs(title = 'Illustration of Haar basis') +
ease_aes('sine-in-out')

#### Approximation via Haar basis

Now, every function can be written as

The coefficients are called the Haar Wavelet coefficients. Let us calculate them in the discrete setting. We give us a mesh and values and some function :

j_max <- 12
n_max <- 13  # maximal resolution
delta <- 2^{-n_max}
x <- (1: 2^n_max)/2^n_max
y <- x* sin(1/x)
alpha <- list()  # list of Haar coefficients
for (j in 0:j_max){
alpha[[j+1]] <- rep(0, 2^j)
for (k  in 0:(2^j-1))
{
alpha[[j+1]][k+1] <-sum(haar(x,j,k)*y)*delta  # approximation of scalar product
}

}

y_approx <- rep(0, length(y))  # approximated values of y

#### Calculating the approximations for different values of j

df <- data.frame(x=numeric(),
y=numeric(),
id=numeric())

y_approx <-  mean(y) # this is the contribution of the constant function
for (j in 0:j_max){
for (k  in 0:(2^j-1))
{
y_approx <- y_approx + alpha[[j+1]][k+1]*haar(x,j,k)
}
df_new <- data.frame(x=x,
y=y_approx,
id=j)
df <- bind_rows(df, df_new)
}

## Animation of convergence


ggplot(df, aes(x=x, y=y)) +
geom_step()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = "Haar basis approximation of x sin(1/x). Depth:  {closest_state}") +
ease_aes('sine-in-out')

## Lévy’s construction of Brownian motion

Brownian motion is a central object of probability theory. The idea of Lévy’s construction is to construct the Brownian motion step by step on finite sets

of dyadic points. As is dense in the Brownian motion is then obtained as the uniform limit of linear interpolation on .
It is pretty easy to illustrate this construction using R and the package gganimate. We use the same notation as in the proof of Wiener’s theorem given on page 23 in “Brownian motion” by Peter Mörters and Yuval Peres.

library(dplyr)
library(gganimate)
library(transformr)
set.seed(42)
n_max <- 14                # maximal number of steps
D <- (0: 2^n_max)/2^n_max   # this is in fact D_n_max
B <- list()
Z <- rnorm(1)
B[[1]] <- c(0, Z/2 + rnorm(1), Z)
for (n in 2:n_max){
B[[n]] <- rep(0, 2^n+1)
index_known <- seq(1,2^n+1, by=2) # indices where values are known from previous steps
B[[n]][index_known] <- B[[n-1]]
index_unknown <- seq(2, 2^n, by=2) # indices where values are not yet defined
for (d in index_unknown){
B[[n]][d] <- 0.5*(B[[n]][d-1]+B[[n]][d+1])+ rnorm(1)*2^(-(n+1)/2)
}
}

## interpolation and transformation in a data.frame
df <- data.frame(time=numeric(),
value=numeric(),
id=numeric())

for (n in 1:n_max){
D_n<-(0: 2^n)/2^n
B_interpol<- approx(D_n, B[[n]], xout = D)\$y   # interpolation
df_new <- data.frame(time=D,
value=B_interpol,
id=n)
df <- bind_rows(df, df_new)
}

## animation
ggplot(df, aes(x=time, y=value)) +
geom_line()+
transition_states(
id,
transition_length = 2,
state_length = 1
) +
labs(title = 'Levys construction of Brownian motion.  Step: {closest_state}', x = 'time', y = 'position') +
ease_aes('sine-in-out')`

## What are the probable outcomes of the whole tournament?

There a billions of different outcomes of the 2018 FIFA World Cup™ and it is hard to get an overview on who is going to play against who in the different stages of the tournament. Most forecasts models only give the most probable outcomes. However, if this forecast is wrong at some stage (this will happen with a probability close to 1!) all other forecasts for later stages are useless. For example, let us consider a model predicting that Germany and Brazil will both win their group. If then one of the two teams do not win their group all games involving group E and F in the knock-out phase become useless. One way out of the dilemma is to forecast all different outcomes of the tournament and assign probabilities to each outcome. This is done in the following Sankey diagram. Please click on the graphic to obtain an interactive java script version. In this version you can move the different “sites” and hover over the “links” to see the corresponding probabilities. In this way you can “zoom in” the parts you’re most interested in. Continue reading “What are the probable outcomes of the whole tournament?”

## A more sophisticated forecast model

Football is a typical low-scoring game and games are frequently decided through single events in the game. These events may be extraordinary individual performances, individual errors, injuries, refereeing errors or just lucky coincidences. Moreover, during a tournament there are most of the time teams and players that are in exceptional shape and have a strong influence on the outcome of the tournament. One consequence is that every now and then alleged underdogs win tournaments and reputed favorites drop out already in the group phase.

The above effects are notoriously difficult to forecast. Despite this fact, every team has its strengths and weaknesses (e.g. defense and attack) and most of the results reflect the qualities of the teams. In order to model the random effects and the deterministic drift forecasts should be given in terms of probabilities.

A series of statistical models have been proposed in the literature for the prediction of football outcomes. They can be divided into two broad categories. The first one, the result-based model, models directly the probability of a game outcome (win/draw/loss), while the second one, the score-based model, focusses on the match score. We want to  follow the second approach since the match score is important in the group phase of the championship and it also implies a model for the first one.  The model proposed in How to impress your football fan colleagues is a first approach but it does only give the most typical outcomes.

The chances are very low (in fact almost zero) that all matches during the world cup end with these results. There are several models for this purpose and most of them involve a Poisson model. In other words, the distribution of the goals of a team is supposed to follow a Poisson distribution. This distribution is determined by one parameter called that describes the expected number of goals.

## How to impress your football fan colleagues…

This post is dedicated to Nina, she knows why.

… with zero-knowledge on football. Before and during the 2018 FIFA World Cup™ all your colleagues, friends or even family member talk about football? Who is going to win? What a surprise that team blabla lost, what do you think? Again such a close match! You do not see any way to avoid it? Impress them using the following flow chart. There are reasonable chances that your forecasts outperform the subjective and pretentious estimates of your colleagues.

The above flowchart is based on the results obtained in What are typical football results? and What are typical football results II?. Using these data you can easily cook up more sophisticated models. For more details on this series of posts see Who wins the 2018 FIFA World Cup™?.

## What are typical football results II?

We continue What are typical football results? The notion of weaker and stronger has not been made precise. Is is true that a team that has say 5 Elo points more than another team is really stronger? What might be an appropriate threshold? A glance at the current Elo ranking might give an indication that teams in within 50 points may be considered as equally strong. But is this true? At which threshold the probabilities of win, draw, lose will change?

## What are typical football results?

In this post we continue our investigation of Who wins the 2018 FIFA World Cup™? and take a first look on historical data of FIFA football matches. These are obtained from the site www.eloratings.net using the wayback machine and some copy-and-paste. Unfortunately, our data set obtained in this way is not complete and we did not obtain data on all FIFA matches in this millennium. However, we were able to retrieve all matches of the FIFA World cup 2018 participants plus the matches of Italy, the Netherlands, and Austria. Yes, Italy and the Netherlands are not qualified, but we still are convinced that these two teams are amongst the strongest teams in the world. We added Austria to pay homage to the country where we spent a lot of quality time.

We will try to answer questions like:

• What is the most probable outcome of a game? [->]
• What is the probability to have a win, a draw or a lose? [->]
• What is the probability that the stronger team wins? [->] And with what result? [->]

The answers to these questions can be found following the links after the questions. Detailed answers can be found below.