# Learning Julia by an example of a simple Bayesian linear regression

**URL:** <https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811>\
**Category:** Performance\
**Tags:** question, performance, distributions, bayesian-inference, linear-regression\
**Created:** [June 3, 2023, 2:13pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811 "2023-06-03T14:13:53Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 3, 2023, 2:13pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/1 "2023-06-03T14:13:53Z")

</div>

Hello everyone,

long-time lurker, very occasional poster. I decided to give Julia a try once more with 1.9. So far I have struggled with the language on and off (I have touched on this [previously](https://discourse.julialang.org/t/why-is-julia-so-great/94718/48)) but would really like to eventually get the hang of it. I feel that my biggest obstacles are that it’s hard for me to find the best way to write _julian_ code in Julia 🙂

To get better, I have written a small piece of code that is an example encompassing most of the code that I have to work with on a regular basis. The original code from which I wrote this one is a teaching code - focusing more on clarity than speed. I still think it is a good exercise. Also, the Julia code takes 30% of the time of the original code, which is in Matlab, so that is inspiring.

**My aim with this topic is to get some general guidance how can I make this code the best it can be and general feedback of what I should watch out for and how I can go about it.**

## Example summary

The structure of the code is Bayesian linear regression

Y = X \beta + \varepsilon \qquad \varepsilon \sim N(0,\Sigma)

where we have to estimate the parameters \beta and \Sigma

Thus, the code has the following steps:

1. Import data
2. Create Y and X
3. Setup a big loop where using conditional distributions:
  1. Draw \beta and save it
  2. Draw \Sigma and save it

4. Do some summaries, plots, etc.

This is pretty much it. I might have many more steps in complex codes usually but _the jist is always the same_.

## The code structure

- I have a module, named NewB, where I keep the 4 functions for my example: `mlag`, `genBeta`, `genSigma`, `gibbs` that perform the aforementioned steps 2 and 3.
- The main code, in file MainNewB.jl reads the data and defines some preliminary stuff needed for the main module.

[I have put the project on github](https://github.com/borisblagov/Julia_AR4_Bayesian_Regression) (new to this as well) and replaced the original data with random numbers. Of your fork it and you prefer the data, uncomment the lines.

I am attaching the code here.

Using `@time` macro gives me 0.07 second on my machine with 7000 draws total (vs 0.2 in Matlab so yay!).

## Specific questions

1. `gibbs` is the main function. Is saving the \beta and \Sigma optimal? Typically I have many more parameters to save and the vectors and matrices are much larger

```julia-auto
beta_dist = zeros(5,n_gibbs-burn)
if i > burn
            beta_dist[:,i-burn] = beta_d
end

```

Should I preallocate memory, as in Matlab for that? The `@time` macro shows a lot of allocations, especially if I go to 70 000 replications. I know I want to save many matrices but does this point to a problem?

```julia-auto
julia> @time include("MainNewB.jl")
  0.082845 seconds (171.56 k allocations: 68.732 MiB, 18.70% gc time)

```

1. I am not sure I understand the output of `@code_warntype` for the `gibbs` function. Highlights union of nothing and Tuple at some point. Is this something to look out for?

```julia-auto
Locals
  @_11::Union{Nothing, Tuple{Int64, Int64}}

```

1. If I have one-liners for specific tasks (such as drawing sigma), is putting it in a function still a good idea, is there overhead associated with it? If I understand the performance tip [function barriers](https://docs.julialang.org/en/v1/manual/performance-tips/#kernel-functions) this is always a good idea.
2. Is there a way for module NewB to inherit `using LinearAlgebra` or is this not a good idea? What happens if I have a module needing another module and vice versa in the same time?

## Code main file:

```julia-auto
using Revise
using LinearAlgebra
using Distributions
using DelimitedFiles
using BenchmarkTools
using NewB

# fdata = readdlm("gdp4795.txt")
# Z = 100*fdata[21:end,:]./fdata[20:end-1,:].-100
# plot(Z)
Z = rand(150,1)

n_gibbs = 7000
burn = 2000
p = 4;

(X,Y) = mlag(Z,p)

Sigma0 = I(5)*4000
BETA0 = zeros(5,1) 
sig2_d = 0.5        
nu0 =0
d0 =0

beta_d = genBeta(X,Y,BETA0,Sigma0,sig2_d)
sig2_d = genSigma(Y,X,beta_d,nu0,d0)
(beta_dist, sigma_dist) = gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)

mean(beta_dist,dims=2)
mean(sigma_dist,dims=2)

```

## Module NewB.jl:

```julia-auto
module NewB
using LinearAlgebra
using Distributions
export mlag, genBeta, genSigma, gibbs

"""
    mlag(Yfull::Matrix{Float64},p::Integer)
    Creates lags of a matrix for a VAR representation with a constant on the left
"""
function mlag(Yfull::Matrix{Float64},p::Integer)
    (Tf, n) = size(Yfull)
    X = ones(Tf-p,1)
    for i = 1:p
        X = [X Yfull[p-i+1:end-i,:]]
    end
    Y = Yfull[p+1:end,:]
    return X, Y # this changes the array passed into the function
end

function genBeta(X,Y,Beta_prior,Sigma_prior,sig2_d)
    invSig = Sigma_prior^-1
    V = (invSig + sig2_d^(-1)*(X'*X))^-1
    C = cholesky(Hermitian(V)) 
    Beta1 = V*(invSig*Beta_prior + sig2_d^(-1)*X'*Y)
    beta_d = Beta1 + C.L*randn(5,1)
    return beta_d
end

function genSigma(Y,X,beta_d,nu0,d0)
    nu1 = size(Y,1)+nu0
    d1 = d0 + only((Y-X*beta_d)'*(Y-X*beta_d)) 
    sig2_inv = rand(Gamma(nu1/2,2/d1),1)
    sig2_d = 1/only(sig2_inv)
    return sig2_d
end

function gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)
    beta_dist = zeros(5,n_gibbs-burn)
    sigma_dist = zeros(1,n_gibbs-burn)
    for i = 1:n_gibbs
        beta_d = genBeta(X,Y,BETA0,Sigma0,sig2_d)
        sig2_d = genSigma(Y,X,beta_d,nu0,d0)
        if i > burn
            beta_dist[:,i-burn] = beta_d
            sigma_dist[1,i-burn] = sig2_d
        end
    end
    return beta_dist, sigma_dist
end

end # module NewB

```

If you have gotten this far, thank you very much!

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [June 3, 2023, 3:33pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/2 "2023-06-03T15:33:07Z")

</div>

Because `genBeta` and `genSigma` are called for every iteration you probably want to benchmark them carefully and profile both execution time and memory allocations in them. In `genSigma`, for example, you evaluate the residual, `Y - X * beta_d`, twice and create a matrix product that will be 1 by 1 then extract the scalar value. That could be replaced by

```julia
    resid = Y - X * beta_d
    d1 = d0 + dot(resid, resid)

```

or

```julia
    d1 = d0 + sum(abs2, resid)

```

Similarly, `sig2_d` can be written

```julia
   sig2_d = inv(rand(Gamma(nu1 / 2, 2 / d1)))

```

However, I suspect that the time is being spent in `genBeta` and there the construction of `V` is likely the bottleneck, as that one line violates the first two of Nick Higham’s [seven sins of numerical linear algebra](https://nhigham.com/2022/10/11/seven-sins-of-numerical-linear-algebra/). The fact that a formula is written in terms of a inverse of a matrix constructed from the inverse of an `X'X`-like construction doesn’t mean it is a good idea to try to calculate it that way.

It is likely that the triangular factor `C.L` can be determined from a QR factorization of X directly.

In each call to `genBeta` you are creating the return vector, which is then copied into a column of `beta_dist` in the `gibbs` function. A common Julia idiom is to define a mutating version of a function, say `genBeta!`, which takes as an argument `beta_d` and overwrites its contents, which are then returned. That way you can reuse a working vector declared in the `gibbs` function.

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 4, 2023, 7:54am UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/3 "2023-06-04T07:54:40Z")

</div>

Oh, those or some nice tips, thank you. I will incorporate them and report here my journey.

> [@dmbates](#):
>
> A common Julia idiom is to define a mutating version of a function, say `genBeta!`, which takes as an argument `beta_d`

This tip sound invaluable, I doubt I would have found it out in a document somewhere (or payed attention). The thing is that my lord codes are just versions of that one with more things going on so such best practices are exactly what I am looking for.

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 4, 2023, 2:20pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/4 "2023-06-04T14:20:07Z")

</div>

Thank you, I had the opportunity to try out some of the suggestions.

For example, the `resid = Y - X * beta_d;` trick reduced the allocations of sigma from 7 to 44, I would watch out for such points in the future.

```julia
mine
julia> @time genSigma(Y,X,beta_d,nu0,d0)
  0.000013 seconds (7 allocations: 5.016 KiB)
0.10240049561807199

dot trick:
julia> @time genSigma(Y,X,beta_d,nu0,d0)
  0.000016 seconds (4 allocations: 2.516 KiB)
0.08709325471913029

```

and simply that reduced the total code by almost a half (and allocations by a lot)!

```julia
julia> @time include("MainNewB.jl")
  0.048585 seconds (150.56 k allocations: 51.639 MiB, 9.80% gc time)
vs
julia> @time include("MainNewB.jl")
  0.074861 seconds (171.56 k allocations: 68.732 MiB, 11.86% gc time)

```

I am at that amount of seconds that I should start using `BenchmarkTools`

I am unsure how to incorporate the `"!"` suggestion (just changing the values of `beta_d` instead of re-allocating it. I tried following the advice in [this discussion](https://discourse.julialang.org/t/how-to-define-a-function-such-that-its-arguments-can-mutate/66320/15), more specifically to use `x=f(x)` [as in this post](https://discourse.julialang.org/t/how-to-define-a-function-such-that-its-arguments-can-mutate/66320/12)

I tried pre-defining beta\_d in `gibbs` such that

```julia
function gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)
    beta_d = zeros(5,1)
    ...
    for i = 1:n_gibbs
        beta_d = genBeta(X,Y,BETA0,Sigma0,sig2_d,beta_d)
        ...
    end

```

and then in `genBeta` assigning to beta\_d’s elements but that doesn’t seem to be it.

```julia
beta_d[:,1] = Beta1 + C.L*randn(5,1)

```

alternatively I tried `beta_d .=`.

What is the correct way to go about this? Just a reference to the manual can also help, I couldn’t find easily functions with the `"!"` operator. I looked at the [`broadcast`](https://docs.julialang.org/en/v1/manual/functions/#man-vectorized) part and I thought I can follow but obviously I didn’t implement it well as I don’t see any changes to the allocations, nor time.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 4, 2023, 6:16pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/5 "2023-06-04T18:16:14Z")

</div>

Hey there @Boris,

It’s your lucky day, blabbering about Julia best practices is a favorite activity of mine. You’ll find some information in the official docs, for instance

- [Workflow Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/workflow-tips/)
- [Style Guide · The Julia Language](https://docs.julialang.org/en/v1/manual/style-guide/)
- [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/)

But it’s a bit buried and not necessarily easy to dive into. I tried to reformulate some of this advice in my (WIP) [IntroJulia](https://gdalle.github.io/IntroJulia/). In particular, check out [IntroJulia - package dev](https://gdalle.github.io/IntroJulia/package.html) and [IntroJulia - performance](https://gdalle.github.io/IntroJulia/performance.html). It sums up most of the generic advice I could give you.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 4, 2023, 6:24pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/6 "2023-06-04T18:24:23Z")

</div>

Now onto some of your specific questions:

> [@Boris](#):
>
> I am not sure I understand the output of `@code_warntype` for the `gibbs` function. Highlights union of nothing and Tuple at some point. Is this something to look out for?

Usually this denotes iteration, and you’re fine. Anything that shows up blue in `@code_warntype` is not a cause for worry. If you want an output that more closely matches your source code (instead of the low level representation that is harder to decipher), check out [GitHub - JuliaDebug/Cthulhu.jl: The slow descent into madness](https://github.com/JuliaDebug/Cthulhu.jl). A bit scary at first but remarkably easy to get a handle on. Basically a `@code_warntype` on steroids with arbitrary depth.

> [@Boris](#):
>
> If I have one-liners for specific tasks (such as drawing sigma), is putting it in a function still a good idea, is there overhead associated with it? If I understand the performance tip [function barriers](https://docs.julialang.org/en/v1/manual/performance-tips/#kernel-functions) this is always a good idea.

You also have to take into account clarity here. Separate functions are usually a good thing… up to a point. Not sure what that one-liner is, but there are many situations where `sigma = rand(Normal(0, 1))` would be clearer that `sigma = draw_sigma()` (dumb example). So if it doesn’t affect your code speed, you shouldn’t always bother.

> [@Boris](#):
>
> Is there a way for module NewB to inherit `using LinearAlgebra` or is this not a good idea? What happens if I have a module needing another module and vice versa in the same time?

What do you mean by “inherit”? You already put `using LinearAlgebra` in your `NewB` module, and if you need it in your main file then put it there again. The module is designed to be reusable by someone else, who may not import `LinearAlgebra` in their main file after all.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 4, 2023, 6:25pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/7 "2023-06-04T18:25:12Z")

</div>

As for the actual code, since you said you’re a GitHub newbie and I’m an evil person, I’m gonna open some GitHub issues to get you used to them 😈

> [@Boris](#):
>
> What is the correct way to go about this?

See the GitHub issues I opened

---

<div class="post-metadata">

**Author:** ![wc4wc4wc4](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wc4wc4wc4/32/23038_2.png) [@wc4wc4wc4](https://discourse.julialang.org/u/wc4wc4wc4)\
**Post date:** [June 4, 2023, 7:56pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/8 "2023-06-04T19:56:47Z")

</div>

Just to be sure; you do know about [Turing.jl](https://turinglang.org/dev/tutorials/05-linear-regression#model-specification), right? 🙂  
I assume that the goal of this post is to get better acquainted with the do and don’ts of Julia performance-wise, however, I just wanted to be sure that you knew it.

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 4, 2023, 8:51pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/9 "2023-06-04T20:51:43Z")

</div>

Oh boy, thank you! I will go through it, hopefully this week. I have read actually and tons of stuff about Julia over the years but because I do it sporadically (on intense intervals quite far apart) I have a lot of diminishing returns.

The docs are extremely dense - every sentence has so much information to it that it becomes tiring after a while. Especially, when every line is spending 10-15 minutes and after 30 minutes I don’t feel like I have any progress. And because I am doing all this on the side in my free time I am not getting very far ahead. Even the note taking is time consuming. I understand it doesn’t work that way, I am just annoyed that I don’t get it the first time 😃

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 4, 2023, 8:58pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/10 "2023-06-04T20:58:16Z")

</div>

I actually had no idea, thank you! I mean, I’ve seen the package name but never really checked it out.

But yes, my goal is to actually learn Julia. My research is currently about Bayesian Mixed Frequency Vectorautoregressions and using them for doing economic forecasts and nowcasts and I will be implementing those. For example I have [this working paper already implemented in Matlab.](https://joshuachan.org/papers/BVAR-MF-R1.pdf).

But I am happy to look into the codes of others, I learn a lot, so I will see what Turing has to offer and see if I can directly implement stuff from there.

---

<div class="post-metadata">

**Author:** ![Boris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/boris/32/3306_2.png) [@Boris](https://discourse.julialang.org/u/Boris)\
**Post date:** [June 4, 2023, 8:58pm UTC](https://discourse.julialang.org/t/learning-julia-by-an-example-of-a-simple-bayesian-linear-regression/99811/11 "2023-06-04T20:58:49Z")

</div>

Hahah, that’s the best way to learn, so please do, I will definitely play along!

The longer post I have to digest more and would reply later.