By: Estadistika -- Julia

Re-posted from: https://estadistika.github.io//bayesian/probability/modeling/inference/julia/2020/03/01/Bayesian-Sequential-Learning.html

In my previous article, I discussed how we can use Bayesian approach in estimating the parameters of the model. The process revolves around solving the following conditional probability, popularly known as the Bayes’ Theorem:

$$\begin{equation} \mathbb{P}(\mathbf{w}|\mathbf{y})=\frac{\mathbb{P}(\mathbf{w})\mathbb{P}(\mathbf{y}|\mathbf{w})}{\mathbb{P}(\mathbf{y})}, \end{equation}$$

where $\mathbb{P}(\mathbf{w})$ is the a priori (prior distribution) for the objective parameters, $\mathbb{P}(\mathbf{y}|\mathbf{w})$ is the likelihood or model evidence, and $\mathbb{P}(\mathbf{y})$ is the normalizing constant with the following form:

$$\begin{equation} \mathbb{P}(\mathbf{y})=\int_{\forall\mathbf{w}\in\mathscr{P}}\mathbb{P}(\mathbf{w})\mathbb{P}(\mathbf{y}|\mathbf{w})\operatorname{d}\mathbf{w}, \end{equation}$$

where $\mathscr{P}$ is the parameter space.

Posterior Distribution

The details on the derivation of the a posteriori were also provided in the said article, but there were missing pieces, which I think is necessary for us to support our proposition, and thus we have the following result:

Simulation Experiment

The above result can be applied to any linear models (cross-sectional or time series), and I’m going to demonstrate how we can use it to model the following simulated data. I will be using Julia in Nextjournal (be sure to head over to this link for reproducibility), which already has an image available for version 1.3.1. Having said, some of the libraries are already pre-installed in the said image, for example Plots.jl. Thus, we only have to install the remaining libraries that we will need for this experiment.

Load the libraries as follows:

The theme above simply sets the theme of the plots below. Further, for reproducibility purposes, I provided a seed as well. The following function will be used to simulate a cross-sectional data with population parameters $w_0\triangleq-.3$ and $w_1\triangleq-.5$ for 20 sample size.

From the above results, the parameters of the a posteriori can be implemented as follows:

One feature of Julia is it supports unicode, making it easy to relate the codes to the math above, i.e. Σ and μ in the code are obviously $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ above, respectively. The vector operations in Julia are also cleaner compared to that in R and in Python, for example we can encode Julia’s A'A above as t(A) %*% A in R, and A.T.dot(A) in Python.

Finally, we simulate the data as follows:

We don’t really have to write down the wtrue variable above since that’s the default values of w0 and w1 arguments, but we do so just for emphasis.

While the main subject here is Bayesian Statistics, it would be better if we have an idea as to what the Frequentist solution would be. As most of you are aware of, the solution to the weights above is given by the following normal equation:
$\hat{\mathbf{w}}=(\boldsymbol{\mathfrak{A}}^{\text{T}}\boldsymbol{\mathfrak{A}})^{-1}\boldsymbol{\mathfrak{A}}^{\text{T}}\mathbf{y},$
this is a known result, and we can prove this in a separate article. The above equation is implemented as follows:

Therefore, for the current sample data, the estimate we got when assuming the weights as fixed and unknown is $\hat{\mathbf{w}}=[-0.32264, -0.59357]^{\text{T}}$. The figure below depicts the corresponding fitted line.


Not bad given that we have very small dataset. Now let’s proceed and see how we can infer the parameters in Bayesian framework. The prior distribution can be implemented as follows:

As indicated in the above proposition, the parameters are jointly modelled by a bivariate Normal distribution, as indicated by the dimension of the hyperparameter μ above. Indeed, the true parameter we are interested in is the weight vector, $\mathbf{w}$, but because we considered it to be random, then the parameters of the model we assign to it are called hyperparameters, in this case the vector $\boldsymbol{\mu}$ and the identity matrix $\mathbf{I}$ of the prior distribution.

Moreover, the likelihood of the data can be implemented as follows:

The mean of the likelihood is the specified linear model itself, which in vector form is the inner product of the transformed design matrix, $\boldsymbol{\mathfrak{A}}$, and the weights vector, $\mathbf{w}$, i.e. A[i, :]'w. This assumption is valid since the fitted line must be at the center of the data points, and that the error should be random. One of my favorite features of Julia language is the multiple dispatch. For example, the two likelihood functions defined above are not in conflict since Julia evaluates the inputs based on the type of the function arguments. The same is true for the posterior distribution implemented below. Unlike in R and in Python, I usually have to write this as a helper function, e.g. likelihood_helper.

Finally, the prediction is done by sampling the weights from the posterior distribution. The center of these weights is of course the mean of the a posteriori.

For example, to sample 30 weights from the posterior distribution using all the sampled data, and return the corresponding predictions, is done as follows:

The predicted function returns both x and y values, that’s why we indexed the above result to 2, to show the predicted ys only. Further, to use only the first 10 observations of the data for calculating the $\hat{y}$, is done as follows:

So you might be wondering, what’s the motivation of only using the first 10 and not all observations? Well, we want to demonstrate how Bayesian inference learns the weights or the parameters of the model sequentially.

Visualization

At this point, we are ready to generate our main vis. The first function (dataplot) plots the generated fitted lines from the a priori or a posteriori, both of which is plotted using the extended method of contour we defined below.

Tying all the codes together, gives us this beautiful grid plot.


Since we went with style before comprehension, let me guide you then with the axes. All figures have unit square space, with contour plots having the following axes: $w_0$ (the x-axis) and $w_1$ (the y-axis). Obviously, the data space has the following axes: predictor (the x-axis) and response (the y-axis).

Discussions

We commenced the article with emphasis on the approach of Bayesian Statistics to modeling, whereby the estimation of the parameters as mentioned is based on the Bayes’ Theorem, which is a conditional probability with the following form:
\begin{equation}
\mathbb{P}(\mathbf{w}|\mathbf{y})=\frac{\mathbb{P}(\mathbf{w})\mathbb{P}(\mathbf{y}|\mathbf{w})}{\mathbb{P}(\mathbf{y})}.
\end{equation}
Now we will relate this to the above figure using some analogy, as to how the model sequentially learns the optimal estimate for the parameter of the linear model.

Consider the following: say you forgot where you left your phone, and for some reason you can’t ring it up, because it could be dead or can’t pick up some signals. Further, suppose you don’t wanna look for it, rather you let Mr. Bayes, your staff, to do the task. How would he then proceed? Well, let us consider the weight vector, $\mathbf{w}\triangleq[w_0,w_1]^{\text{T}}$, be the location of your phone. In order to find or at least best approximate the exact location, we need to first consider some prior knowledge of the event. In this case, we need to ask the following questions: where were you the last time you had the phone? Were you in the living room? Or in the kitchen? Or in the garage? And so on. In the context of modeling, this prior knowledge about the true location can be described by a probability distribution, and we refer to this as the a priori (or the prior distribution). These set of possible distributions obviously are models itself with parameters, as mentioned above, referred to as the hyperparameters, which we can tweak to describe our prior knowledge of the event. For example, you might consider the kitchen as the most probable place where you left your phone. So we adjust the location parameter of our a priori model to where the kitchen is. Hence Mr. Bayes should be in the kitchen already, assessing the coverage of his search area. Of course, you need to help Mr. Bayes on the extent of the coverage. This coverage or domain can be described by the scale parameter of your a priori. If we relate this to the main plot, we assumed the prior distribution over the weight vector to be standard bivariate Gaussian distribution, centered at zero vector with identity variance-covariance matrix. Since the prior knowledge can have broad domain or coverage on the possible values of the weight vector, the samples we get generates random fitted lines as we see in the right-most plot of the first row of the figure above.

Once we have the prior knowledge in place, that is, we are already in the kitchen and we know how wide the search area likely to be, we can start looking for evidence. The evidence is the realization of your true model, relating to the math above, these realizations are the $y_i$s, coming from $y_i=f(x|\mathbf{w})$, where $f(x|\mathbf{w})$ is the link function of the true model, which we attempt to approximate with our hypothesized link function, $h(x|\hat{\mathbf{w}})$, that generated the predicted $\hat{y}_i$s. For example, you may not know where exactly your phone is, but you are sure with your habits. So you inform Mr. Bayes, that the last time you were with your phone in the kitchen was drinking some coffee. Mr. Bayes will then use this as his first evidence, and assess the likelihood of each suspected location in the kitchen. That is, what is the likelihood that a particular location, formed (or realized, generated or connected to) the first evidence (taking some coffee)? For example, is it even comfortable to drink coffee in the sink? Obviously not, so very low likelihood, but likely in the dining table or close to where the coffee maker is. If we assess all possible location within our coverage using the first evidence, we get the profile likelihood, which is what we have in the first column of the grid plot above, profile likelihood for the ith evidence. Further, with the first evidence observed, the prior knowledge of Mr. Bayes needs to be updated to obtain the posterior distribution. The new distribution will have an updated location and scale parameters. If we relate to the above figure, we can see the samples of the fitted lines in the data space plot (third column, second row), starting to make guesses of possible lines given the first evidence observed. Moving on, you inform Mr. Bayes of the second evidence, that you were reading some newspaper while having some coffee. At this point, the prior belief of Mr. Bayes, for the next posterior, will be the posterior of the first evidence, and so the coverage becomes restrictive and with new location, which further help Mr. Bayes on managing the search area. The second evidence, as mentioned, will then return a new posterior. You do this again and again, informing Mr. Bayes of other evidences sequentially until the last evidence. The final evidence will end up with the final posterior distribution, which we expect to have new location parameter, closer to the exact location, and small scale parameter, covering the small circle of the exact solution. The final posterior will then be your best guess that would describe the exact location of your phone.

This may not be the best analogy, but that is how the above figure sequentially learns the optimal estimate for the weight vector in Bayesian framework.

Bayesian Deep Learning

This section deserves a separate article, but I will briefly give some motivation on how we can generalize the above discussion into complex modeling.

The intention of the article is to give the reader a low-level understanding of how the Bayes’ theorem works, and without loss of generalization, I decided to go with simple linear regression to demonstrate the above subject. However, this can be applied to any model indexed by or a function of some parameters or weights $\mathbf{w}$, with the assumption that the solution is random but govern by some probability distribution.

Complex modeling such as in Deep Learning are usually based on the assumption that the weights are fixed and unknown, which in Statistics is the Frequentist approach to inference, but without assuming some probability distribution on the error of the model. Therefore, if we are to assume some randomness on the weights, we can then use Bayesian inference to derive or at least approximate (for models with no closed-form solution) the posterior distribution. Approximate Bayesian inference are done via Markov Chain Monte Carlo (MCMC) or Variational Inference, which we can tackle in a separate post.

Libraries

There are several libraries for doing Bayesian inference, the classic and still one of the most powertful library is Stan. For Python, we have PyMC3, Pyro (based on Pytorch), and TensorFlow Porbability. For Julia, we have Turing.jl, Mamba.jl, Gen.jl, and Stan.jl. I will have a separate article for these libraries.

Next Steps

The obvious next steps for readers to try out is to model the variance as well, since in the above result, the variance of the innovation or the error is known and is equal to $\alpha$. Further, one might consider the Frequentist sequential learning as well. Or proceed with other nonlinear complex models, such as Neural Networks. We can have these in a separate article.

Software Versions

Differentiation is a central problem in many fields including deep learning, finance, scientific computing and others. The two conventional techniques used include:

• Symbolic differentiation which can result in complex and redundant expressions

• Numerical differentiation which can lead to large numerical errors

Automatic differentiation (AD) has been behind recent advances in deep learning. A Differentiable Programming System to Bridge Machine Learning and Scientific Computing, makes AD a first class feature in Julia language , helps differentiate programs, rather than simple network architectures.

Julia helps bring AD to the world of machine learning beyond neural networks. We will look at a specific use case and see how Julia allows us to compose functionalities from multiple packages seamlessly. (JuliaLang: The Ingredients for a Composable Programming Language)

Introduction

One of the advantages of Julia is its composability. Automatic differentiation can be applied in a conventional machine learning tool like xgboost. These two packages were developed independently yet they work together seamlessly. In order to illustrate this let’s look at how easily custom loss functions can be implemented in XGBoost with Zygote.

There are a lot of in-built loss functions in xgboost, but these may be suboptimal for a lot of real world problems. Consider testing for a rare disease. False positives are not good, but they are not costly mistakes in the sense that further tests can easily rule them out. On the other hand false negatives are far worse than a false positive. A person who has a condition and needs immediate care has been incorrectly identified as healthy. False negatives are far more costly than false positives. These kinds of scenarios are common in real world machine learning problems. This is where custom loss functions are useful.

The default loss functions available off-the-shelf might not be optimal for the business objectives. On the other hand, a custom loss function that closely matches the business objectives can take the asymmetry on errors into account. In the disease detection scenario discussed above, false negative errors are far more costly than false positives. Using a custom loss function that penalizes the model heavily for making false negative errors will result in a model that is averse to false negatives.

Problem formulation

We have chosen to predict the survival chances on the Titanic ocean liner using a supervised ML technique called XGBoost. It is an implementation of gradient boosting where an ensemble of weak decision tree learners are combined to produce a strong model.

In order to deal with asymmetric penalties, like the case of disease detection that was discussed above, XGBoost permits custom loss functions. First we build a simple baseline model. Let’s see how we do this in Julia.

Load Data

The datasets have been downloaded from Kaggle:

We use CSVFiles to load the training data, and DataFrames formanipulation.

julia> using CSVFiles, DataFrames
julia> df = DataFrame(CSVFiles.load("train.csv"))

The names method will give the features in this dataframe.

  julia> names(df)
  12-element Array{Symbol,1}:
  :PassengerId
  :Survived
  :Pclass
  :Name
  :Sex
  :Age
  :SibSp
  :Parch
  :Ticket
  :Fare
  :Cabin
  :Embarked

Survived is the target variable, we will have to clean and process the rest of the features before a model can be built on it. We will be building a simple model with a few features, which are Age, Embarked, Sex, Pclass, SibSp, Parch, Fare

Data preprocessing

We will have to clean up the data a bit. The following line gives the count of number of missing values in column Embarked

julia> sum(df[:,:Embarked] .== "")
  2

XGBoost does support learning on missing values, but in this case there are only 2 of them. So, it is better that we impute them with the most frequent value in the column.

Replace missing values with the most frequent port:

julia> df[df[:,:Embarked] .== "", :Embarked] = "S"

The column Age also has missing values, the following command tells us that there are 177 rows in the Age column with missing values

julia> sum(ismissing.(df[:Age]))
177

Age being a numerical feature, it is better that we use the average value of Age column to impute the missing values

Replace missing values with the average:

julia> using Statistics
julia> average_age = mean(df[.!ismissing.(df[:Age]), :Age])
julia> df[ismissing.(df[:Age]), :Age] = average_age

We can use one-hot encoding for categorical features such as Pclass and Embarked. One-hot encoding creates a binary feature for each value of the categorical feature. For instance, the feature Embarked has 3 fields, viz. S, C, Q One hot encoding will create 3 binary features, each corresponding to one of the unique values of the categorical column. The binary feature will Embarked will have 1s in rows where the value of Embarked was S, and 0s everywhere in other rows.

julia>for i in unique(df.Pclass)
	df[:,Symbol("Pclass_"*string(i))] = Int.(df.Pclass .== i)
     end

julia>for i in unique(df.Embarked)
	df[:,Symbol("Embarked_"*string(i))] = Int.(df.Embarked .== i)
	end

Gender can be encoded as a binary feature:

julia> gender_dict = Dict("male"=>1, "female"=>0);
julia> df[:Sex] = map(x->gender_dict[x], df[:Sex]);

Building our baseline model

Let’s split the dataset into training and validation set. The training
set can be created as below

julia> x_train = convert(Matrix{Float32},select(df[1:800,:],Not(:Survived)))
julia> y_train = convert(Array{Float32}, df[1:800,:Survived])

Validation set:

julia> x_val = convert(Matrix{Float32},select(df[801:end,:],Not(:Survived)))
julia> y_val = convert(Array{Float32}, df[801:end,:Survived])

Create a DMatrix:

julia> train_dmat = DMatrix(x_train, label=y_train)

Train the model:

julia> bst_base = xgboost(train_dmat,2, eta=0.3, objective="binary:logistic", eval_metric="auc")
[1] train-auc:0.893250
[2] train-auc:0.899080

Get predictions on the validation set:

julia> ŷ = predict(bst, x_val)

The following function will calculate the accuracy and weighted f score:

function evaluate(y, ŷ;threshold=0.5)
	out = zeros(Int64, 2,2)
	ŷ = Int.(ŷ.>=threshold)
	out[1,1]=sum((y.==0).&(ŷ.==0))
	out[2,2]=sum((y.==1).&(ŷ.==1))
	out[2,1]=sum((y.==1).&(ŷ.==0))
	out[1,2]=sum((y.==0).&(ŷ.==1)	
	r0 = out[1,1]/(out[1,1]+out[1,2])
	p0 = out[1,1]/(out[1,1]+out[2,1])
	f0 = 2*p0*r0/(p0+r0	
	r1 = out[2,2]/(out[2,2]+out[2,1]
	p1 = out[2,2]/(out[2,2]+out[1,2]
	f1 = 2*r1*p1/(p1+r1	
	println("Weighted f1 = ", round((sum(y .== 0.0)/length(y)) * f0 + (sum(y .== 1.0)/length(y)) * f1	digits=3))
	println("Accuracy =", (out[2,2]+out[1,1])/sum(out))
	out
end                                                                      

Let’s look at the performance of the baseline model:

julia> evaluate(y_val, ŷ)
Weighted f1 = 0.845
Accuracy = 0.8461538461538461
2×2 Array{Int64,2}:
51 6
 8 26

Let’s submit this to Kaggle to get a tangible measure of the model performance. The current model is at 12300^th^ position on the leaderboard

Custom loss function

There are 8 false negatives, in order to reduce false negatives, we can weigh false negatives higher than false positives in our loss function.

The signature of custom loss should be:

function weighted_loss(preds::Vector{Float32}, dtrain::DMatrix)
		gradients = … #calculate gradients
		hessians = … #calculate hessians
		return gradients, hessians
end

The conventional approach is to calculate the gradients and Hessians by hand and then translate them to code. Let’s start by looking at the log loss function used by the model.

Notice that the term -yln(sigma(x)) which penalizes false negatives has a coefficient of 1.

We can change this weight to make the model penalize false negatives more than false positives. With a weight w on false positives the equation will look like this:

XGBoost admits loss functions in the signature that we defined above. We will have to find the gradient and Hessian of this loss and plug them in the function.

The gradient turns out to be:

Differentiating the above again gives us the second derivative:

Now we can complete our function with a slightly higher penalty, i.e. 1.5:

function weighted_loss(preds::Vector{Float32}, dtrain::DMatrix)
	beta = 1.5
	p = 1. ./ (1 .+ exp.(-preds))
	grad = p .* ((beta .- 1) .* y .+ 1) .- beta .* y
	hess = ((beta .- 1) .* y .+ 1) .* p .* (1.0 .- p)
	return grad, hess
end

Note that this required us to do some cumbersome calculus to implement the custom loss. A smarter approach would be to invoke Julia’s automatic differentiation to calculate gradients and Hessians. Apart from not having to do the math, it gives us flexibility as well since a slight change in the function means we won’t have to redo the math again.

Automatic differentiation

XGBoost outputs scores that need to be passed through a sigmoid function. We do this inside the custom loss function that we defined above. Let’s define it here explicitly:

σ(x) = 1/(1+exp(-x))

The weighted log loss can be defined as:

weighted_logistic_loss(x,y) = -1.5 .* y*log(σ(x)) -1 .* (1-y)*log(1-σ(x))

We have added a weight of 1.5 to false negatives.

Now we are ready to take advantage of Zygote to calculate the gradient and Hessian:

gradient_logistic(x,y) = gradient(weighted_logistic_loss,x,y)[1]

The gradient method differentiates the weighted_logistic_loss function. The parameters of the weighted_logistic_loss function are also passed alongside. We then take the first element with [1] because we are interested in the derivative with respect to the first parameter x.

Similarly, the Hessian can be calculated by differentiating the above gradient function:

hess_logistic(x,y) = gradient(grad_logistic,x,y)[1]

We can define the custom loss function using this gradient and Hessian:

function custom_objective(preds::Vector{Float32}, dtrain::DMatrix)
  y = get_info(dtrain, "label")
  grad = grad_logistic.(preds,y)
  hess = hess_logistic.(preds,y)
  return grad,hess
end

This can be passed to the XGBoost trainer:

julia> bst = xgboost(train_dmat, 2,eta=0.3, eval_metric="auc", obj=custom_objective)
[1] train-auc:0.892897
[2] train-auc:0.899565

At this point we can evaluate the model, and compare it against the baseline:

julia> ŷ = predict(bst, x_val)
julia> evaluate(y_val, ŷ)
Weighted f1 = 0.878
Accuracy =0.8791208791208791
2×2 Array{Int64,2}:
53 4
7 27

That’s quite an improvement for a small change in the loss function. We didn’t have to do any of the taxing math to try out this objective function. All we had to do was define the function and let Zygote take care of calculating the gradient and Hessian.

At this point, we can submit the results from the new model to Kaggle. We have been able to move 6400 places up on the leaderboard with very little effort.

This exercise illustrates the composability of Julia packages and the flexibility that it gives data scientists.