Julia Community 🟣: Steven Siew

How to use Matrices to perform Least Square Fitting of experimental data using non-linear model

Steven Siew — Wed, 03 Apr 2024 01:39:53 +0000

Explaining how the algorithm works and the concepts involved is quite difficult and so I shall do it slowly and section by section. I shall refer to the reader as Alice as in the famous character in the Novel "Alice through the Looking Glass". The reason for this is because I need to take the reader into the Looking Glass Universe where things are view differently

Assumptions

In this article I assume that You know

How to add a package to Julia
What an XY graph is
What a Matrix is
The basic idea of what an inverse of a Matrix is
The basic idea of what a transpose of a Matrix is
Basic maths. Addition Subtration. Multiplication. Division. Power.

In this article, I will show you how to perform Least Square Fitting to obtain the parameters of a non-linear model on experimental data.

Step 0 Talk abut the model

In this article we shall model the experimental data using the non-linear model

Y = A exp(B * X) + C

where X,Y are the measurement values from the experiment and

A,B,C are the unknown constants.

The unknown constants are also refer to as the parameters

We generalise this by saying

Y = f(x) and

f(x) = A exp(B * x) + C

Now this function f is very important and it represent the relationship between X and Y. That is to say: Given X, what is Y?

Later on in the article we will use the function f to get a "Looking Glass" function g by stepping through the Looking Glass. Naturally the function g is highly related to the function f . We shall talk about this later.

Step 1 Obtaining experimental data

We could perform the experiment and obtain lots of values of X and Y but for this article we shall generate the values using a deterministic random number generator.

Here is the source code for generating the values of X and Y

using Distributions, StableRNGs

#=
We want to generate Float64 data deterministically
=#

function generatedata(formula; input=[Float64(_) for _ in 0:(8-1)], seed=1234, digits=0,
    preamble="rawdata = ", noisedist=Normal(0.0, 0.0), verbose=true)
    vprint(verbose, str) = print(verbose == true ? str : "")
    local myrng = StableRNG(seed)
    local array = formula.(input) + rand(myrng, noisedist, length(input))
    if digits > 0
        array = round.(array, digits=digits)
    end
    vprint(verbose, preamble)
    # Now print the begining array symbol '['
    vprint(verbose, "[")
    local counter = 0
    local firstitem = true
    for element in array
        if firstitem == true
            firstitem = false
        else
            vprint(verbose, ",")
        end
        counter += 1
        if counter % 5 == 1
            vprint(verbose, "\n  ")
        end
        vprint(verbose, element)
    end
    vprint(verbose, "]\n")
    return array
end

A_actual = 1.5
B_actual = -0.25
C_actual = 3.5
f(x) = A_actual * exp(B_actual * x) + C_actual
numofsamples = 400+1
X = [Float64(_) for _ in 0:0.01:(numofsamples-1)/100]
Y = generatedata(f, input=X, digits=5, noisedist=Normal(0.0, 0.0005), preamble="Y = ",verbose=false)
println("\n\nX = ", X, "\nY = ", round.(Y,digits=5), "\n\n")

with the output

X = [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4, 1.41, 1.42, 1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5, 1.51, 1.52, 1.53, 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6, 1.61, 1.62, 1.63, 1.64, 1.65, 1.66, 1.67, 1.68, 1.69, 1.7, 1.71, 1.72, 1.73, 1.74, 1.75, 1.76, 1.77, 1.78, 1.79, 1.8, 1.81, 1.82, 1.83, 1.84, 1.85, 1.86, 1.87, 1.88, 1.89, 1.9, 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97, 1.98, 1.99, 2.0, 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09, 2.1, 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19, 2.2, 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3, 2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4, 2.41, 2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5, 2.51, 2.52, 2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6, 2.61, 2.62, 2.63, 2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7, 2.71, 2.72, 2.73, 2.74, 2.75, 2.76, 2.77, 2.78, 2.79, 2.8, 2.81, 2.82, 2.83, 2.84, 2.85, 2.86, 2.87, 2.88, 2.89, 2.9, 2.91, 2.92, 2.93, 2.94, 2.95, 2.96, 2.97, 2.98, 2.99, 3.0, 3.01, 3.02, 3.03, 3.04, 3.05, 3.06, 3.07, 3.08, 3.09, 3.1, 3.11, 3.12, 3.13, 3.14, 3.15, 3.16, 3.17, 3.18, 3.19, 3.2, 3.21, 3.22, 3.23, 3.24, 3.25, 3.26, 3.27, 3.28, 3.29, 3.3, 3.31, 3.32, 3.33, 3.34, 3.35, 3.36, 3.37, 3.38, 3.39, 3.4, 3.41, 3.42, 3.43, 3.44, 3.45, 3.46, 3.47, 3.48, 3.49, 3.5, 3.51, 3.52, 3.53, 3.54, 3.55, 3.56, 3.57, 3.58, 3.59, 3.6, 3.61, 3.62, 3.63, 3.64, 3.65, 3.66, 3.67, 3.68, 3.69, 3.7, 3.71, 3.72, 3.73, 3.74, 3.75, 3.76, 3.77, 3.78, 3.79, 3.8, 3.81, 3.82, 3.83, 3.84, 3.85, 3.86, 3.87, 3.88, 3.89, 3.9, 3.91, 3.92, 3.93, 3.94, 3.95, 3.96, 3.97, 3.98, 3.99, 4.0]
Y = [5.00026, 4.99671, 4.99168, 4.98814, 4.98473, 4.98102, 4.97771, 4.97393, 4.96957, 4.96567, 4.96252, 4.96015, 4.9556, 4.95183, 4.94868, 4.94439, 4.94139, 4.93792, 4.93397, 4.93044, 4.92662, 4.92324, 4.91925, 4.91613, 4.91287, 4.90931, 4.90609, 4.90147, 4.89842, 4.89508, 4.89157, 4.8876, 4.88476, 4.88175, 4.87781, 4.87462, 4.87074, 4.86746, 4.86411, 4.86027, 4.8575, 4.85359, 4.84989, 4.84682, 4.84404, 4.83992, 4.83673, 4.83395, 4.83077, 4.82657, 4.82314, 4.82066, 4.81735, 4.81311, 4.80956, 4.80774, 4.80396, 4.80025, 4.79692, 4.79453, 4.79056, 4.78834, 4.78521, 4.78146, 4.77902, 4.77504, 4.77291, 4.76972, 4.7647, 4.76251, 4.75913, 4.75643, 4.75262, 4.7496, 4.74612, 4.74393, 4.74016, 4.73748, 4.73441, 4.73105, 4.72772, 4.72505, 4.72198, 4.71788, 4.71669, 4.71308, 4.71046, 4.70721, 4.70349, 4.70032, 4.69721, 4.69474, 4.69145, 4.68943, 4.68538, 4.68237, 4.68024, 4.67703, 4.67413, 4.67229, 4.66888, 4.66589, 4.66264, 4.65933, 4.65577, 4.65449, 4.65059, 4.64823, 4.6461, 4.6431, 4.6392, 4.63667, 4.63245, 4.63044, 4.62859, 4.62477, 4.62197, 4.62009, 4.61724, 4.61425, 4.6117, 4.6092, 4.60616, 4.60314, 4.59975, 4.5978, 4.59554, 4.59203, 4.58831, 4.58639, 4.58263, 4.58192, 4.57885, 4.5761, 4.57306, 4.57052, 4.56772, 4.56495, 4.56203, 4.55984, 4.55683, 4.55441, 4.55088, 4.54908, 4.54606, 4.54433, 4.54214, 4.53836, 4.5362, 4.53241, 4.53141, 4.52867, 4.52551, 4.52338, 4.52062, 4.51883, 4.51532, 4.51249, 4.51068, 4.50698, 4.50629, 4.50336, 4.50098, 4.49811, 4.49508, 4.49317, 4.48983, 4.48791, 4.48504, 4.48337, 4.48136, 4.4778, 4.47559, 4.47316, 4.4703, 4.46867, 4.46528, 4.46452, 4.45991, 4.45953, 4.45685, 4.45328, 4.45157, 4.44927, 4.44621, 4.44493, 4.4421, 4.43922, 4.43794, 4.43555, 4.43322, 4.43068, 4.42809, 4.4264, 4.42316, 4.42128, 4.41914, 4.41739, 4.41406, 4.41157, 4.40969, 4.4076, 4.40453, 4.40268, 4.40096, 4.39855, 4.39654, 4.39399, 4.39193, 4.3892, 4.38635, 4.38445, 4.38312, 4.38146, 4.37867, 4.37581, 4.37449, 4.37264, 4.3692, 4.36779, 4.36556, 4.36366, 4.36077, 4.3589, 4.3578, 4.35527, 4.35167, 4.34997, 4.3491, 4.3457, 4.34382, 4.34193, 4.3396, 4.33865, 4.33591, 4.33407, 4.33229, 4.32937, 4.32821, 4.3259, 4.32237, 4.32102, 4.31951, 4.31657, 4.31483, 4.31291, 4.31125, 4.30921, 4.30641, 4.30484, 4.3028, 4.301, 4.29897, 4.29732, 4.29448, 4.29281, 4.28975, 4.28915, 4.2868, 4.2858, 4.28295, 4.28058, 4.28009, 4.27766, 4.27502, 4.27358, 4.27174, 4.2687, 4.26786, 4.26582, 4.26392, 4.26183, 4.25931, 4.25732, 4.25707, 4.25473, 4.25133, 4.24999, 4.24787, 4.24655, 4.24507, 4.24296, 4.24149, 4.23917, 4.23702, 4.2357, 4.23442, 4.23145, 4.23031, 4.22895, 4.22615, 4.22453, 4.22319, 4.22086, 4.21921, 4.21783, 4.21609, 4.21384, 4.21163, 4.21076, 4.20879, 4.2065, 4.20493, 4.20258, 4.20095, 4.19956, 4.19858, 4.19528, 4.19366, 4.19288, 4.1901, 4.18912, 4.1877, 4.18586, 4.18422, 4.18253, 4.18112, 4.17955, 4.17818, 4.17609, 4.17396, 4.17281, 4.17105, 4.16854, 4.16699, 4.16667, 4.16367, 4.16227, 4.16073, 4.15925, 4.15801, 4.15625, 4.15331, 4.15193, 4.15034, 4.14909, 4.14783, 4.14581, 4.14392, 4.14243, 4.1407, 4.13978, 4.13744, 4.13692, 4.13441, 4.13327, 4.13116, 4.13048, 4.12872, 4.12701, 4.1243, 4.12402, 4.12258, 4.11982, 4.11875, 4.11768, 4.11669, 4.1146, 4.11309, 4.11081, 4.10998, 4.10841, 4.10661, 4.10511, 4.10405, 4.10142, 4.10056, 4.09889, 4.09851, 4.09635, 4.09528, 4.09188, 4.09195, 4.0895, 4.08858, 4.08764, 4.0859, 4.08428, 4.08316, 4.08134, 4.07975, 4.07809, 4.07754, 4.07452, 4.0756, 4.07306, 4.07144, 4.06975, 4.06911, 4.06743, 4.06596, 4.06462, 4.06311, 4.06288, 4.05962, 4.05898, 4.05684, 4.05688, 4.05526, 4.05294, 4.05135]

Step 2 Doing the Least Square Fitting on the experimental data

We show you the source code for the algorithmn and the output first and later in the article we shall explain how the algorithmn works.

"""
Provide the solution to the nonlinear Least Square
"""
function nonlinearLSQ(xarray, yarray; maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ], stepsize=1.0, verbose=true)
    len_xarray = length(xarray)
    len_yarray = length(yarray)
    if len_xarray != len_yarray
        println("Error in function nonlinearLSQ: length of xarray != length of yarray")
        return nothing
    end
    # Parameter array P = [ A , B , C ]
    P = startpos  # Initial Starting Position
    for iteration = 1:maxiter
        R = zeros(len_yarray) # a len_yarray size array
        # Let's create the J Matrix (Jacobian Matrix)
        J = zeros(len_yarray, 3) # a len_yarray*3 matrix
        for i = 1:len_yarray
            J[i, 1] = exp(P[2] * X[i])
            J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]
            J[i, 3] = 1.0
            R[i] = yarray[i] - (  P[1] * exp(P[2] * xarray[i]) + P[3]  )
        end
        # Dont forget to reshape values from 
        # a 1D vector R into a len*1 column matrix RR
        local RR = reshape(R, len_yarray, 1) # a len_yarray * 1 Matrix
        DeltaP = inv(J' * J) * J' * RR 
        # At this point the varible DeltaP is a len_yarray * 1 Matrix
        # We turn it into an 1 dimensional vector using vec() function
        P = P + stepsize .* vec( DeltaP )
        if verbose == true
            println("Iteration is ",iteration)
            println(" DeltaP is ",round.(DeltaP,digits=5))
            println(" P is ",round.(P,digits=5))
        end
    end
    println()
    return P
end

solution = nonlinearLSQ(X, Y, maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ])
println("Trying with startpos [ 1.0 , -0.1 , 1.0 ] stepsize=1.0 maxiter=8")
println("The estimate for A is ", round.(solution[1],digits=5))
println("The estimate for B is ", round.(solution[2],digits=5))
println("The estimate for C is ", round.(solution[3],digits=5))

With the output

Iteration is 1
 DeltaP is [-1.58274; -0.42322; 4.57972;;]
 P is [-0.58274, -0.52322, 5.57972]
Iteration is 2
 DeltaP is [1.82102; -0.53788; -1.82729;;]
 P is [1.23828, -1.0611, 3.75242]
Iteration is 3
 DeltaP is [-0.34721; 0.78483; 0.28928;;]
 P is [0.89107, -0.27627, 4.0417]
Iteration is 4
 DeltaP is [0.601; 0.04232; -0.53397;;]
 P is [1.49208, -0.23395, 3.50773]
Iteration is 5
 DeltaP is [0.00431; -0.01652; -0.00425;;]
 P is [1.49638, -0.25047, 3.50349]
Iteration is 6
 DeltaP is [0.00429; 0.00069; -0.00425;;]
 P is [1.50067, -0.24979, 3.49923]
Iteration is 7
 DeltaP is [1.0e-5; -0.0; -1.0e-5;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 8
 DeltaP is [0.0; 0.0; -0.0;;]
 P is [1.50068, -0.24979, 3.49923]

Trying with startpos [ 1.0 , -0.1 , 1.0 ] stepsize=1.0 maxiter=8
The estimate for A is 1.50068
The estimate for B is -0.24979
The estimate for C is 3.49923

Step 3 Discussing how Jacobian Matrix is constructed in the code.

In the source code, it was mentioned that you need to create the Jacobian Matrix. Now you may want to read the wikipedia article below on how to create the Jacobian Matrix but I found the wikipedia article to be confusing for newbies.

https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant

The reason the wikipedia article is confusing is that it talks about a function f but that f is referred to a particular generic function related to the issue at hand.

I have modified the content of the wikipedia article to better reflect our problem and it looks like this.

Now time for some explaination. The Matrix p is a 3x1 parameter matrix. It consists of

Now the function g on the otherhand is a family of functions.

Because we have 401 data points from the experiment. Namely

Now Alice, we come back to the function f which we mentioned at the start of the article. As you recall, we had the non-linear model

Y = A exp(B * X) + C

where X,Y are the measurement values from the experiment and

A,B,C are the unknown constants.

Which we generalise by saying

Y = f(x) and

f(x) = A exp(B * x) + C

Now Alice, understand that this is TRUE before we perform the experiment. That is to say X and Y are variables while A, B and C are the unknown constants.

That is to say
X can be any value(s)
Y can be any value(s)
but
A is a fixed value (but currently unknown)
B is a fixed value (but currently unknown)
C is a fixed value (but currently unknown)

And the WHOLE POINT of the experiment is to determine the values of A, B and C

Now it is time to enter the Looking Glass by performing the experiment. Once we had perform the experiment, we are in the MAGICAL REALM of the LOOKING GLASS and things change:

X and Y are now a list of FIXED value pairs

x1 , y1
x2 , y2
.
.
.
x401 , y401

We represent this by saying x and y is a 401x1 column matrix

And we represent it by

y = g(p)

Where g is a family of functions consists of

g is a 1x401 matrix of

g1(A,B,C) = A exp(B * x1) + C
g2(A,B,C) = A exp(B * x2) + C
g3(A,B,C) = A exp(B * x3) + C
.
.
.
g401(A,B,C) = A exp(B * x401) + C

Which in this case, we can say g(A,B,C) is a function related to the function of f(x)

Now when we put in the values of x1 , x2 , x3 , ... , x401 from the experimental results, we get

g1(A,B,C) = A exp(B * 0.0) + C
g2(A,B,C) = A exp(B * 0.01) + C
g3(A,B,C) = A exp(B * 0.02) + C
.
.
.
g401(A,B,C) = A exp(B * 4.0) + C

It is important to know WHAT HAS CHANGED

That is to say X and Y are constants while A, B and C are now variables

That is to say
X is FIXED (by the values obtained from the experiment)
Y is FIXED (by the values obtained from the experiment)
but
A can be any value
B can be any value
C can be any value

We need to do this because we wanted to differentiate the variables A , B and C and we cannot (technically) differentiate a constant.

So back to generating the Jacobian Matrix

We need to find the partial derivative of g1 in respect to A

in particular

g1(A,B,C) = A exp(B * x1) + C

There are 3 ways

Using Mathematica
Using Wolfram Alpha
Doing it by hand

Using Mathematica

Using Wolfram Alpha

using url https://www.wolframalpha.com

Doing it by hand

Too lazy to do it by hand, Here are the basic steps

The equivalent part of the source code is

J[i, 1] = exp(P[2] * X[i])

Remember P[2] is the variable B

Partial derivative with respect to B

Next We need to find the partial derivative of g1 in respect to B

in particular

g1(A,B,C) = A exp(B * x1) + C

The answer is A * exp(B * x1) * x1

Which is reflected in the source code

J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]

Partial derivative with respect to C

Finally We need to find the partial derivative of g1 in respect to C

in particular

g1(A,B,C) = A exp(B * x1) + C

The answer is 1

Which is reflected in the source code

J[i, 3] = 1.0

Step 4 Discussing how the algorithm works. Multiple explainations. Pick one that you like

How does the algorithmn work?

There are many explanation of the algorithm works. Here are some of the explainations

First we calculate

r_i = y_i - g_i(p)

We want g_i(p) to be exactly y_i so

r_i = 0 (since g_i(p) is exactly y_i)

Now this is a function r_i(p) which we can diffentiate using the Jacobian matrix.

Now, what is a Jacobian? A Jacobian is the equivalence of a multidimensional derivative analog of a function expressed in a matrix format. With the importance on the fact that it is multidimensional.

Jacobian(r) = Jacobian(y - g(p))
Jacobian(r) = Jacobian(-g(p)) since y is a constant
Jacobian(r) = -1 * Jacobian(g(p))

therefore

p_1 = p_0 - Jacobian(r(p_0))
p_1 = p_0 + Jacobian(g(p_0))

Please note that given p_i the Jacobian(r(p_i)) points to a direction that is the maximum increases of r(p_i) (ie maximum gradient).

in general, if you start with a parameter estimate p_i , you can get a better estimate with.

p_(i+1) = p_i - Jacobian(r(p_i))
thus
p_(i+1) = p_i + Jacobian(g(p_i))

So we can continously update the value of p

Step 5 Discussing the effect of the step size

Our starting pos is

solution = nonlinearLSQ(X, Y, maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ])

And what happens next is

pos[0] is [ 1.0 , -0.1 , 1.0 ]

Our first step is
DeltaP is [-1.58274; -0.42322; 4.57972;;]
pos[1] is [-0.58274, -0.52322, 5.57972]

Our second step is
DeltaP is [1.82102; -0.53788; -1.82729;;]
pos[2] is [1.23828, -1.0611, 3.75242]

Our third step is
DeltaP is [-0.34721; 0.78483; 0.28928;;]
pos[3] is [0.89107, -0.27627, 4.0417]

Our fourth step is
DeltaP is [0.601; 0.04232; -0.53397;;]
pos[4] is [1.49208, -0.23395, 3.50773]

Our fifth step is
DeltaP is [0.00431; -0.01652; -0.00425;;]
pos[5] is [1.49638, -0.25047, 3.50349]

Our sixth step is
DeltaP is [0.00429; 0.00069; -0.00425;;]
pos[6] is [1.50067, -0.24979, 3.49923]

Our seventh step is
DeltaP is [1.0e-5; -0.0; -1.0e-5;;]
pos[7] is [1.50068, -0.24979, 3.49923]

Our eigth step is
DeltaP is [0.0; 0.0; -0.0;;]
pos[8] is [1.50068, -0.24979, 3.49923]

Noticed that our step are gradually getting smaller and smaller.

So why do we start with startpos = [ 1.0 , -0.1 , 1.0 ]? Because I knew the answer is 1.5 , -0.25 , 3.5 . I know the answer is that because I builted up the raw data in the first place.

But what if we do not know what the answer is?

Then I woud have started the program with the startpos of [ 1.0 , -1.0 , 1.0 ] being a generic starting point. I tend to avoid values of 0.0

So what happens when we start with that startpos instead? We get

Trying with startpos [ 1.0 , -1.0 , 1.0 ] stepsize=1.0 maxiter=8

pos[0] is [ 1.0 , -1.0 , 1.0 ]
Iteration is 1
DeltaP is [-0.08045; 0.88542; 3.02142;;]
pos[1] is [0.91955, -0.11458, 4.02142]
Iteration is 2
DeltaP is [-0.71245; -0.37299; 0.76899;;]
pos[2] is [0.20711, -0.48757, 4.79041]
Iteration is 3
DeltaP is [1.06606; 1.33297; -1.0707;;]
pos[3] is [1.27317, 0.8454, 3.71972]
Iteration is 4
DeltaP is [-1.46498; 0.03173; 1.33635;;]
pos[4] is [-0.19181, 0.87712, 5.05607]
Iteration is 5
DeltaP is [0.01914; -0.19153; -0.0255;;]
pos[5] is [-0.17267, 0.68559, 5.03056]
Iteration is 6
DeltaP is [-0.15717; -0.37635; 0.19506;;]
pos[6] is [-0.32984, 0.30924, 5.22562]
Iteration is 7
DeltaP is [-1.23965; -0.65218; 1.30354;;]
pos[7] is [-1.5695, -0.34294, 6.52916]
Iteration is 8
DeltaP is [3.00069; -0.07667; -2.96156;;]
pos[8] is [1.43119, -0.41961, 3.5676]

The estimate for A is 1.43119
The estimate for B is -0.41961
The estimate for C is 3.5676

The important thing is that we get the WRONG answer. But why? The ANSWER is that we keep overstepping the correct answer and we didn't take enough steps. If we tried with MORE steps:

solution = nonlinearLSQ(X, Y, startpos=[ 1.0 , -1.0, 1.0 ],
maxiter=18,
verbose=true)

Trying with startpos [ 1.0 , -1.0 , 1.0 ] stepsize=1.0 maxiter=18
Iteration is 1
 DeltaP is [-0.08045; 0.88542; 3.02142;;]
 P is [0.91955, -0.11458, 4.02142]
Iteration is 2
 DeltaP is [-0.71245; -0.37299; 0.76899;;]
 P is [0.20711, -0.48757, 4.79041]
Iteration is 3
 DeltaP is [1.06606; 1.33297; -1.0707;;]
 P is [1.27317, 0.8454, 3.71972]
Iteration is 4
 DeltaP is [-1.46498; 0.03173; 1.33635;;]
 P is [-0.19181, 0.87712, 5.05607]
Iteration is 5
 DeltaP is [0.01914; -0.19153; -0.0255;;]
 P is [-0.17267, 0.68559, 5.03056]
Iteration is 6
 DeltaP is [-0.15717; -0.37635; 0.19506;;]
 P is [-0.32984, 0.30924, 5.22562]
Iteration is 7
 DeltaP is [-1.23965; -0.65218; 1.30354;;]
 P is [-1.5695, -0.34294, 6.52916]
Iteration is 8
 DeltaP is [3.00069; -0.07667; -2.96156;;]
 P is [1.43119, -0.41961, 3.5676]
Iteration is 9
 DeltaP is [-0.08597; 0.14295; 0.08345;;]
 P is [1.34522, -0.27666, 3.65105]
Iteration is 10
 DeltaP is [0.14663; 0.02842; -0.14309;;]
 P is [1.49185, -0.24823, 3.50796]
Iteration is 11
 DeltaP is [0.00879; -0.00157; -0.00869;;]
 P is [1.50064, -0.2498, 3.49926]
Iteration is 12
 DeltaP is [4.0e-5; 1.0e-5; -4.0e-5;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 13
 DeltaP is [-0.0; -0.0; 0.0;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 14
 DeltaP is [0.0; 0.0; -0.0;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 15
 DeltaP is [0.0; -0.0; 0.0;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 16
 DeltaP is [-0.0; -0.0; 0.0;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 17
 DeltaP is [0.0; 0.0; -0.0;;]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 18
 DeltaP is [0.0; 0.0; -0.0;;]
 P is [1.50068, -0.24979, 3.49923]

The estimate for A is 1.50068
The estimate for B is -0.24979
The estimate for C is 3.49923

Then we get the right ANSWER

Or we could try with a half step size

solution = nonlinearLSQ(X, Y, startpos=[ 1.0 , -1.0, 1.0 ],
maxiter=(8*2), stepsize=(1.0/2),
verbose=true )

Trying with startpos [ 1.0 , -1.0 , 1.0 ] stepsize=0.5 maxiter=(8*2)
Iteration is 1
 DeltaP is [-0.08045; 0.88542; 3.02142;;]
 P is [0.95978, -0.55729, 2.51071]
Iteration is 2
 DeltaP is [0.24708; 0.36408; 1.2708;;]
 P is [1.08332, -0.37525, 3.14611]
Iteration is 3
 DeltaP is [0.31175; 0.14463; 0.45672;;]
 P is [1.23919, -0.30294, 3.37447]
Iteration is 4
 DeltaP is [0.2326; 0.05843; 0.15327;;]
 P is [1.35549, -0.27372, 3.45111]
Iteration is 5
 DeltaP is [0.13803; 0.02526; 0.05519;;]
 P is [1.42451, -0.26109, 3.47871]
Iteration is 6
 DeltaP is [0.07442; 0.01163; 0.02225;;]
 P is [1.46172, -0.25528, 3.48983]
Iteration is 7
 DeltaP is [0.03853; 0.00557; 0.00982;;]
 P is [1.48098, -0.25249, 3.49474]
Iteration is 8
 DeltaP is [0.01959; 0.00272; 0.00459;;]
 P is [1.49078, -0.25113, 3.49704]
Iteration is 9
 DeltaP is [0.00987; 0.00135; 0.00222;;]
 P is [1.49571, -0.25046, 3.49814]
Iteration is 10
 DeltaP is [0.00496; 0.00067; 0.00109;;]
 P is [1.49819, -0.25012, 3.49869]
Iteration is 11
 DeltaP is [0.00248; 0.00033; 0.00054;;]
 P is [1.49943, -0.24995, 3.49896]
Iteration is 12
 DeltaP is [0.00124; 0.00017; 0.00027;;]
 P is [1.50005, -0.24987, 3.49909]
Iteration is 13
 DeltaP is [0.00062; 8.0e-5; 0.00013;;]
 P is [1.50037, -0.24983, 3.49916]
Iteration is 14
 DeltaP is [0.00031; 4.0e-5; 7.0e-5;;]
 P is [1.50052, -0.24981, 3.49919]
Iteration is 15
 DeltaP is [0.00016; 2.0e-5; 3.0e-5;;]
 P is [1.5006, -0.2498, 3.49921]
Iteration is 16
 DeltaP is [8.0e-5; 1.0e-5; 2.0e-5;;]
 P is [1.50064, -0.24979, 3.49922]

The estimate for A is 1.50064
The estimate for B is -0.24979
The estimate for C is 3.49922

We can even try with one-tenth the step size

#=
  Try it with [ 1.0 , -1.0 , 1.0 ] stepsize=0.1
=#

solution = nonlinearLSQ(X, Y, startpos=[ 1.0 , -1.0, 1.0 ], 
    maxiter=(8*10), stepsize=(1.0/10),
    verbose=false)
println("Trying with startpos [ 1.0 , -1.0 , 1.0 ] stepsize=0.1 maxiter=(8*10)")
println("The estimate for A is ", round.(solution[1],digits=5))
println("The estimate for B is ", round.(solution[2],digits=5))
println("The estimate for C is ", round.(solution[3],digits=5))

With the result

Trying with startpos [ 1.0 , -1.0 , 1.0 ] stepsize=0.1 maxiter=(8*10)
The estimate for A is 1.50016
The estimate for B is -0.24992
The estimate for C is 3.49907

Step 6 Further Improvements

One obvious improvement is to compute the two step sizes

full step
half step

and then "decide" which step is the appropriate step.

We can do this by comparing the fitness value of each step by applying a fitness function to each candidate step size. We then choose the step with the greatest fitness value (smallest error).

That is to say "fitness value" = 1.0/abs(error)

"""
Provide the solution to the linear Least Square Version 2
"""
function nonlinearLSQ_ver2(xarray, yarray; maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ], verbose=true)
    len_xarray = length(xarray)
    len_yarray = length(yarray)
    if len_xarray != len_yarray
        println("Error in function nonlinearLSQ: length of xarray != length of yarray")
        return nothing
    end
    # Now define our error function
    function abserrorfunc(xarray,yarray,P)
        local i
        local abserrorvalue = 0.0
        for i = 1:length(yarray)
            abserrorvalue += (yarray[i] - (  P[1] * exp(P[2] * xarray[i]) + P[3]  ))^2.0
        end
        return abserrorvalue
    end
    # Now define our fitness function
    function fitness(xarray,yarray,P)
        return 1.0/abserrorfunc(xarray,yarray,P)
    end
    # Parameter array P = [ A , B , C ]
    P = startpos  # Initial Starting Position
    for iteration = 1:maxiter
        R = zeros(len_yarray) # a len_yarrayx1 column matrix
        # Let's create the J Matrix (Jacobian Matrix)
        J = zeros(len_yarray, 3) # a len_yarray*3 matrix
        for i = 1:len_yarray
            J[i, 1] = exp(P[2] * X[i])
            J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]
            J[i, 3] = 1.0
            R[i] = yarray[i] - (  P[1] * exp(P[2] * xarray[i]) + P[3]  )
        end
        # Dont forget to reshape values from 
        # a 1D vector R into a len*1 column matrix RR
        local RR = reshape(R, len_yarray, 1) # a len_yarray * 1 Matrix
        DeltaP = inv(J' * J) * J' * RR 
        # Define our two stepsizes
        local fullstep = 1.0 .* vec( DeltaP )
        local halfstep = 0.5 .* vec( DeltaP )
        # Calculate fitness value for both stepsize
        local fitnessvalue_fullstep = fitness(xarray,yarray,P + fullstep)
        local fitnessvalue_halfstep = fitness(xarray,yarray,P + halfstep)
        # Choose the step with the greatest fitness value
        local appropriatestep
        if fitnessvalue_fullstep >= fitnessvalue_halfstep
            appropriatestep = "fullstep"
            P = P + fullstep
        else
            appropriatestep = "halfstep"
            P = P + halfstep
        end
        if verbose == true
            println("Iteration is ",iteration)
            println(" DeltaP is ",round.(DeltaP,digits=5))
            println("The appropriate step size is ",appropriatestep)
            println(" P is ",round.(P,digits=5))
        end
    end
    println()
    return P
end

#=
  Try it with [ 1.0 , -1.0 , 1.0 ] 
=#
println("Trying with Version 2 with startpos [ 1.0 , -1.0 , 1.0 ] maxiter=8")
solution = nonlinearLSQ_ver2(X, Y, startpos=[ 1.0 , -1.0, 1.0 ], 
    maxiter=18, 
    verbose=true)
println("")
println("The estimate for A is ", round.(solution[1],digits=5))
println("The estimate for B is ", round.(solution[2],digits=5))
println("The estimate for C is ", round.(solution[3],digits=5))
print("\n\n")

with output

Trying with Version 2 with startpos [ 1.0 , -1.0 , 1.0 ] maxiter=8
Iteration is 1
 DeltaP is [-0.08045; 0.88542; 3.02142;;]
The appropriate step size is fullstep
 P is [0.91955, -0.11458, 4.02142]
Iteration is 2
 DeltaP is [-0.71245; -0.37299; 0.76899;;]
The appropriate step size is halfstep
 P is [0.56333, -0.30108, 4.40592]
Iteration is 3
 DeltaP is [0.91013; 0.12439; -0.87981;;]
The appropriate step size is fullstep
 P is [1.47345, -0.17669, 3.52611]
Iteration is 4
 DeltaP is [-0.13223; -0.0919; 0.13188;;]
The appropriate step size is fullstep
 P is [1.34122, -0.26859, 3.65799]
Iteration is 5
 DeltaP is [0.15487; 0.02025; -0.15423;;]
The appropriate step size is fullstep
 P is [1.4961, -0.24834, 3.50376]
Iteration is 6
 DeltaP is [0.00455; -0.00146; -0.0045;;]
The appropriate step size is fullstep
 P is [1.50064, -0.2498, 3.49926]
Iteration is 7
 DeltaP is [3.0e-5; 1.0e-5; -3.0e-5;;]
The appropriate step size is fullstep
 P is [1.50068, -0.24979, 3.49923]
Iteration is 8
 DeltaP is [-0.0; -0.0; 0.0;;]
The appropriate step size is fullstep
 P is [1.50068, -0.24979, 3.49923]


The estimate for A is 1.50068
The estimate for B is -0.24979
The estimate for C is 3.49923

If you have 3 or more step sizes then you need version 3 of the code below

"""
Provide the solution to the linear Least Square Version 3
"""
function nonlinearLSQ_ver3(xarray, yarray; maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ], verbose=true)
    len_xarray = length(xarray)
    len_yarray = length(yarray)
    if len_xarray != len_yarray
        println("Error in function nonlinearLSQ: length of xarray != length of yarray")
        return nothing
    end
    # Now define our error function
    function abserrorfunc(xarray,yarray,P)
        local i
        local abserrorvalue = 0.0
        for i = 1:length(yarray)
            abserrorvalue += (yarray[i] - (  P[1] * exp(P[2] * xarray[i]) + P[3]  ))^2.0
        end
        return abserrorvalue
    end
    # Parameter array P = [ A , B , C ]
    P = startpos  # Initial Starting Position
    for iteration = 1:maxiter
        R = zeros(len_yarray) # a len_yarrayx1 column matrix
        # Let's create the J Matrix (Jacobian Matrix)
        J = zeros(len_yarray, 3) # a len_yarray*3 matrix
        for i = 1:len_yarray
            J[i, 1] = exp(P[2] * X[i])
            J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]
            J[i, 3] = 1.0
            R[i] = yarray[i] - (  P[1] * exp(P[2] * xarray[i]) + P[3]  )
        end
        # Dont forget to reshape values from 
        # a 1D vector R into a len*1 column matrix RR
        local RR = reshape(R, len_yarray, 1) # a len_yarray * 1 Matrix
        DeltaP = inv(J' * J) * J' * RR 
        # Define our various stepsizes
        local fullstep         = (1.0/1.0)   .* vec( DeltaP )
        local halfstep         = (1.0/2.0)   .* vec( DeltaP )
        local quarterstep      = (1.0/4.0)   .* vec( DeltaP )
        local eighthstep       = (1.0/8.0)   .* vec( DeltaP )
        local sixteenthstep    = (1.0/16.0)  .* vec( DeltaP )
        # Calculate error value for various stepsize
        local errorvalue_fullstep        = abserrorfunc(xarray,yarray,P + fullstep)
        local errorvalue_halfstep        = abserrorfunc(xarray,yarray,P + halfstep)
        local errorvalue_quarterstep     = abserrorfunc(xarray,yarray,P + quarterstep)
        local errorvalue_eighthstep      = abserrorfunc(xarray,yarray,P + eighthstep)
        local errorvalue_sixteenthstep   = abserrorfunc(xarray,yarray,P + sixteenthstep)
        # Create an array of tuples
        local resultarray = Tuple{Float64, Vector{Float64}, String}[]
        push!(resultarray,  (errorvalue_fullstep,fullstep,"fullstep")  )
        push!(resultarray,  (errorvalue_halfstep,halfstep,"halfstep")  )
        push!(resultarray,  (errorvalue_quarterstep,quarterstep,"quarterstep")  )
        push!(resultarray,  (errorvalue_eighthstep,eighthstep,"eighthstep")  )
        push!(resultarray,  (errorvalue_sixteenthstep,sixteenthstep,"sixteenthstep")  )
        # Choose the step with the smallest error value

        local sortedresultarray = sort(resultarray,lt=(A,B)->A[1]<B[1])
        local appropriatestep = sortedresultarray[1][3]
        local correctstep = sortedresultarray[1][2]

        P = P + correctstep

        if verbose == true
            println("Iteration is ",iteration)
            println(" DeltaP is ",round.(DeltaP,digits=5))
            println("The appropriate step size is ",appropriatestep)
            println("The correct step is ",round.(correctstep,digits=5))
            println(" P is ",round.(P,digits=5))
        end
    end
    println()
    return P
end

#=
  Try it with [ 1.0 , -1.0 , 1.0 ] 
=#
println("Trying with Version 3 with startpos [ 1.0 , -1.0 , 1.0 ] maxiter=8")
solution = nonlinearLSQ_ver3(X, Y, startpos=[ 1.0 , -1.0, 1.0 ], 
    maxiter=8, 
    verbose=true)
println("")
println("The estimate for A is ", round.(solution[1],digits=5))
println("The estimate for B is ", round.(solution[2],digits=5))
println("The estimate for C is ", round.(solution[3],digits=5))
print("\n\n")

with the result

Trying with Version 3 with startpos [ 1.0 , -1.0 , 1.0 ] maxiter=8
Iteration is 1
 DeltaP is [-0.08045; 0.88542; 3.02142;;]
The appropriate step size is fullstep
The correct step is [-0.08045, 0.88542, 3.02142]
 P is [0.91955, -0.11458, 4.02142]
Iteration is 2
 DeltaP is [-0.71245; -0.37299; 0.76899;;]
The appropriate step size is quarterstep
The correct step is [-0.17811, -0.09325, 0.19225]
 P is [0.74144, -0.20783, 4.21367]
Iteration is 3
 DeltaP is [0.72118; -0.09448; -0.67662;;]
The appropriate step size is fullstep
The correct step is [0.72118, -0.09448, -0.67662]
 P is [1.46263, -0.30231, 3.53705]
Iteration is 4
 DeltaP is [0.00973; 0.04897; -0.00986;;]
The appropriate step size is fullstep
The correct step is [0.00973, 0.04897, -0.00986]
 P is [1.47236, -0.25334, 3.52719]
Iteration is 5
 DeltaP is [0.02814; 0.0036; -0.02778;;]
The appropriate step size is fullstep
The correct step is [0.02814, 0.0036, -0.02778]
 P is [1.50049, -0.24975, 3.49941]
Iteration is 6
 DeltaP is [0.00018; -4.0e-5; -0.00018;;]
The appropriate step size is fullstep
The correct step is [0.00018, -4.0e-5, -0.00018]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 7
 DeltaP is [0.0; 0.0; -0.0;;]
The appropriate step size is fullstep
The correct step is [0.0, 0.0, -0.0]
 P is [1.50068, -0.24979, 3.49923]
Iteration is 8
 DeltaP is [-0.0; -0.0; 0.0;;]
The appropriate step size is fullstep
The correct step is [-0.0, -0.0, 0.0]
 P is [1.50068, -0.24979, 3.49923]


The estimate for A is 1.50068
The estimate for B is -0.24979
The estimate for C is 3.49923

How to use Matrices to perform Simple Least Square fitting of experimental data

Steven Siew — Tue, 12 Mar 2024 00:40:28 +0000

How to use Matrices to perform Simple Least Square fitting of experimental data

Assumptions
In this article I assume that You know

How to add a package to Julia
What the equation of a straight line is
What a quadratic equation is
What an XY graph is
What a Matrix is
The basic idea of what an inverse of a Matrix is
The basic idea of what a transpose of a Matrix is
Basic maths. Addition Subtration. Multiplication. Division. Power.

Let's make things simple. Let's just say that you are a Medical Researcher and you are performing an experiment. You want to measure a person's

Body Mass Index (BMI)
What percentage of the population in that BMI value has Type 2 Diabetes (You get this number from Medical Census)

You would have an XY graph where

X axis is the BMI
Y axis is the percentage of population in that BMI value that has Type 2 Diabetes

You want to get a linear equation of

Y = M * X + C where both M and C are numeric constants

and both variable X and Y are of various values

In this article, the actual values of X and Y does not matter. What I wanted to show is how to get the value of M and C using matrices in Julia.

First thing we need to do is to get some raw data. Rather then getting REAL DATA, we can just generate them using code from another article of mine earlier.

Let's write some code to generate experimental data X and Y

using Distributions,StableRNGs

#=
We want to generate Float64 data deterministically
=#

function generatedata(formula;input=[Float64(_) for _ in 0:(8-1)],seed=1234,digits=0,
    preamble="rawdata = ",noisedist=Normal(0.0,0.0),verbose=true)
    vprint(verbose,str) = print( verbose == true ? str : "")
    local myrng = StableRNG(seed)
    local array = formula.(input) + rand(myrng,noisedist,length(input))
    if digits > 0
        array = round.(array,digits=digits)
    end
        vprint(verbose,preamble)
    # Now print the begining array symbol '['
    vprint(verbose,"[")
    local counter = 0
    local firstitem = true
    for element in array
        if firstitem == true
            firstitem = false
        else
            vprint(verbose,",")
        end
        counter +=1
        if counter % 5 == 1
            vprint(verbose,"\n  ")
        end
        vprint(verbose,element)
    end
    vprint(verbose,"]\n")
    return array
end

M_actual = 0.4
C_actual = 0.3
f(x) = M_actual * x + C_actual 
numofsamples = 3
X = [ Float64(_) for _ in 0:(numofsamples-1)]
Y = generatedata(f,input=X,digits=2,noisedist=Normal(0.0,0.2),preamble="Y = ")
print("\n\nX = ",X,"\nY = ",Y,"\n\n")

and the output is

Y = [
  0.4,0.88,0.76]


X = [0.0, 1.0, 2.0]
Y = [0.4, 0.88, 0.76]

So we now have three equations (using Y = M * X + C )

0.40 = M * 0.0 + C
0.88 = M * 1.0 + C
0.76 = M * 2.0 + C

we can now rewrite the three equations as

0.40 = M * 0.0 + C * 1.0
0.88 = M * 1.0 + C * 1.0
0.76 = M * 2.0 + C * 1.0

Because multiplying anything by one changes nothing.

Next we change the order of the terms on the right hand side

0.40 = C * 1.0 + M * 0.0 
0.88 = C * 1.0 + M * 1.0 
0.76 = C * 1.0 + M * 2.0

Next we change the order of the multiplication

0.40 = 1.0 * C + 0.0 * M
0.88 = 1.0 * C + 1.0 * M
0.76 = 1.0 * C + 2.0 * M

Now we got it in a form that can be easily converted to Matrix format.

We shall use the following

Y to represent [ 0.4, 0.88, 0.76 ]
K to represent [ 1.0 0.0; 1.0 1.0; 1.0 2.0 ]
P to represent [ C , M ]

julia> K = [ 1.0 0.0; 1.0 1.0; 1.0 2.0 ]
3×2 Matrix{Float64}:
 1.0  0.0
 1.0  1.0
 1.0  2.0

So in Matrix form we have Y = K P
Where

K is the 3x2 matrix form of X
P is the 2x1 matrix form of the parameters C and M

Here is my rough sketch on paper

Now that we got it in Matrix format, all we need to do is to solve for P

This is the WRONG WAY to do it

Show wrong way to solve for P

And here is why it is wrong!

julia> K = [ 1.0 0.0; 1.0 1.0; 1.0 2.0 ]
3×2 Matrix{Float64}:
 1.0  0.0
 1.0  1.0
 1.0  2.0

julia> inv(K)
ERROR: DimensionMismatch: matrix is not square: dimensions are (3, 2)
Stacktrace:
 [1] checksquare
   @ /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:241 [inlined]
 [2] inv(A::Matrix{Float64})
   @ LinearAlgebra /Applications/Julia-1.10.app/Contents/Resources/julia/share/julia/stdlib/v1.10/LinearAlgebra/src/dense.jl:911
 [3] top-level scope
   @ REPL[3]:1

And this is the correct way to do it.

Let's write some code to implement our solution

"""
Provide the solution to the linear Least Square
"""
function linearLSQ(xarray,yarray)
    len_xarray = length(xarray)
    len_yarray = length(yarray)
    if len_xarray != len_yarray
        println("Error in function linearLSQ: length of xarray != length of yarray")
        return nothing
    end
    # Let's create the K Matrix
    K = zeros(len_yarray,2) # a len*2 matrix
    for i = 1:len_yarray
        K[i,1] = 1.0
        K[i,2] = xarray[i]
    end
    # Dont forget to reshape cases from a 1D vector yarray into a len*1 column matrix Y
    local Y = reshape(yarray,len_yarray,1)
    # In Julia K' is the Transpose of K
    P = inv(K'*K) * K' * Y
    return P
end

solution = linearLSQ(X,Y)
println("The estimate for C is ",solution[1])
println("The estimate for M is ",solution[2])

when we run it, we have the solution

The estimate for C is 0.5
The estimate for M is 0.17999999999999994

The reason we do not get the correct values is because we only have 3 experimental samples from our experimental data. We can get a better estimate if we do the experiment again with 100 samples.

Also those of you who is observant will noticed that I use the Adjoint operator and not the transpose operator. This is because when used on REAL NUMBER matrices the adjoint operator is the same as a Transpose operator.

Let's do the experiment again with 100 samples this time

numofsamples = 100
X2 = [ Float64(_) for _ in 0:(numofsamples-1)]
Y2 = generatedata(f,input=X2,digits=2,noisedist=Normal(0.0,0.2),preamble="Y = ")

solution2 = linearLSQ(X2,Y2)
println("The estimate for C is ",solution2[1])
println("The estimate for M is ",solution2[2])

And the results are

Y = [
  0.4,0.88,0.76,1.24,1.76,
  2.16,2.72,3.08,3.21,3.52,
  4.12,5.03,5.07,5.42,6.01,
  6.14,6.78,7.23,7.49,7.91,
  8.21,8.68,8.91,9.48,9.99,
  10.37,10.9,10.85,11.43,11.89,
  12.28,12.48,13.14,13.71,13.92,
  14.42,14.64,15.09,15.52,15.75,
  16.4,16.59,16.86,17.38,18.02,
  18.11,18.57,19.2,19.66,19.71,
  20.06,20.79,21.18,21.2,21.5,
  22.48,22.67,22.89,23.25,24.0,
  24.1,24.9,25.33,25.52,26.22,
  26.31,27.13,27.52,27.18,27.97,
  28.28,28.86,28.99,29.43,29.69,
  30.45,30.59,31.16,31.56,31.85,
  32.15,32.71,33.1,33.08,34.23,
  34.39,34.96,35.27,35.39,35.72,
  36.07,36.68,36.96,37.74,37.71,
  38.09,38.82,39.11,39.53,40.37]
The estimate for C is 0.23681584158414637
The estimate for M is 0.4009188718871888

Now the estimated value for M is almost spot on but the estimated value for C is a bit off.

Using linear least square to predict/interpolate the next value in the sequence.

If the values are subject to only a little bit of noise, we can use linear least square to predict the next value in the sequence. For example:

numofsamples = 3
X3 = [ Float64(_) for _ in 0:(numofsamples-1)]
Y3 = generatedata(f,input=X3,digits=2,noisedist=Normal(0.0,0.02),preamble="Y = ",verbose=false)

println("\n\nX3 = ",X3,"\nY3 = ",Y3,"\n\n")

"""
Interpolate the next value in the sequence
using linear interpolation
"""
function LSQ_linear_nextseq(data; verbose=false)
    len = length(data)
    day = [ _ for _ in 0:(len-1)]

    #=
    using the Ordinary Least Squares Algorithm
    https://en.wikipedia.org/wiki/Linear_least_squares
    https://en.wikipedia.org/wiki/Ordinary_least_squares
    =#
    K = zeros(len,2) # a len*2 matrix
    for i = 1:len
        K[i,1] = 1.0
        K[i,2] = day[i]
    end
    # Dont forget to reshape cases from a 1D vector into a len*1 column matrix
    P = inv(K'*K) * K' * reshape(data,len,1)
    linear_func(t) = P[1] + P[2] * t
    next_value_in_sequence = linear_func(len)

    linear_param_rounded = round.(P,digits=4)
    if verbose
        println("linear model is $(linear_param_rounded[1]) + $(linear_param_rounded[2]) * t")
    end
    return next_value_in_sequence
end

println("The next value in the sequence is interpolated as ",
LSQ_linear_nextseq(Y3,verbose=true))

with the result

X3 = [0.0, 1.0, 2.0]
Y3 = [0.31, 0.72, 1.07]


linear model is 0.32 + 0.38 * t
The next value in the sequence is interpolated as 1.4600000000000002

So that is how to predict/interpolate the next value in the sequence but what about the previous value in the sequence? The easiest way is just to reverse the array and call the next value in the sequence again.

"""
Interpolate the previous value in the sequence
using linear interpolation
"""
function LSQ_linear_prevseq(data; verbose=false)
    newdata = reverse(data)
    value = LSQ_linear_nextseq(newdata)
    return value
end

Quadratic Least Squares

Linear Least Squares is basically modelling the data with a STRAIGHT LINE. But straight line is not the only way you can model experimental data. In high school we learn about quadratic equaion like this

Y = A * X^2 + B * X + C
or
Y = C + B * X + A * X^2

Using the data from https://en.wikipedia.org/wiki/Ordinary_least_squares (see picture below)

Here is the code

X4 = [ 1.47,1.50,1.52,1.55,1.57,1.60,1.63,1.65,
       1.68,1.70,1.73,1.75,1.78,1.80,1.83 ]

Y4 = [ 52.21,53.12,54.48,55.84,57.20,58.57,
       59.93,61.29,63.11,64.47,66.28,68.10,69.92,
       72.19,74.46 ]

"""
Provide the solution to the quadratic Least Square
"""
function quadLSQ(xarray,yarray)
    len_xarray = length(xarray)
    len_yarray = length(yarray)
    if len_xarray != len_yarray
        println("Error in function linearLSQ: length of xarray != length of yarray")
        return nothing
    end
    # Let's create the K Matrix
    K = zeros(len_yarray,3) # a len*3 matrix
    for i = 1:len_yarray
        K[i,1] = 1.0
        K[i,2] = xarray[i]
        K[i,3] = xarray[i]^2.0
    end
    # Dont forget to reshape cases from a 1D vector yarray into a len*1 column matrix Y
    local Y = reshape(yarray,len_yarray,1)
    P = inv(K'*K) * K' * Y
    return P
end

solution4 = quadLSQ(X4,Y4)
print("\n\n")
println("The estimate for C is ",solution4[1])
println("The estimate for B is ",solution4[2])
println("The estimate for A is ",solution4[3])

With the output

The estimate for C is 128.81280358089134
The estimate for B is -143.16202286910266
The estimate for A is 61.96032544388436

Compare our result with the result on the wikipedia page below

Generating Test Data

Steven Siew — Mon, 11 Mar 2024 00:08:24 +0000

Generating Test Data for Julia Source Code

Sometimes you need to constuct some data for your julia program and you want to do it easily. You want the following

Deterministic. You need the same data regardless of which version of Julia you are using. You also want it to be reproducible so that anyone who runs the source code get the exactly same data.
Looks "nice" on your source code. You do not want 15 decimals digits of precision for each data point.
Easy to generate. You do not want to waste 2 hours of your time to generate 200 data points. You want it in 2 minutes.

In this article, we will show some code to do this quickly. This code requires the following packages

StableRNGs
Distributions

Here is the Source code

using Distributions,StableRNGs

#=
We want to generate Float64 data deterministically
=#

function generatedata(formula;input=[Float64(_) for _ in 0:(8-1)],seed=1234,digits=0,
preamble="rawdata = ",noisedist=Normal(0.0,0.0),verbose=true)
    vprint(verbose,str) = print( verbose == true ? str : "")
    local myrng = StableRNG(seed)
    local array = formula.(input) + rand(myrng,noisedist,length(input))
    if digits > 0
        array = round.(array,digits=digits)
    end
        vprint(verbose,preamble)
    # Now print the begining array symbol '['
    vprint(verbose,"[")
    local counter = 0
    local firstitem = true
    for element in array
        if firstitem == true
            firstitem = false
        else
            vprint(verbose,",")
        end
        counter +=1
        if counter % 5 == 1
            vprint(verbose,"\n  ")
        end
        vprint(verbose,element)
    end
    vprint(verbose,"]\n")
    return array
end

generatedata(x->0.4*x+0.3,input=[ Float64(_) for _ in 0:(10-1)],digits=2,noisedist=Normal(0.0,0.2));
println("The end.")

We generate the rawdata using the formula f(x)= 0.4*x + 0.3 using the input data of [0.0, 1.0, 2.0, ... , 9.0] with the additional noise distribution of NormalDistribution(mu=0.0 , sigma=0.2). We also make sure there is only 5 data points per line and each data points have only 2 digits after the decimal point.

Here is the output

rawdata = [
  0.4,0.88,0.76,1.24,1.76,
  2.16,2.72,3.08,3.21,3.52]
The end.

Now you can cut and paste the generated data straight into your julia source code. If you are even lazier, you can also do this:

rawdata = generatedata(x->0.4*x+0.3,input=[ Float64(_) for _ in 0:(10-1)],digits=2,noisedist=Normal(0.0,0.2))

Sets under Julia

Steven Siew — Sun, 26 Nov 2023 04:02:25 +0000

Sets are useful concepts in Julia but few people uses them so in this article we shall describe the basic concept of Sets and how to use them

Construction

julia>  evenunderten = Set([2,4,6,8])
Set{Int64} with 4 elements:
  4
  6
  2
  8

julia> oddunderten = Set([1,3,5,7,9])
Set{Int64} with 5 elements:
  5
  7
  9
  3
  1

You can use basic Set operators

julia> numbersunderten = union(evenunderten,oddunderten)
Set{Int64} with 9 elements:
  5
  4
  6
  7
  2
  9
  8
  3
  1

julia> issubset(evenunderten,numbersunderten)
true

julia> bothevenandoddunderten = intersect(evenunderten,oddunderten)
Set{Int64}()

julia> primesunderten = Set([2,3,5,7])
Set{Int64} with 4 elements:
  5
  7
  2
  3

julia> evenandprime = intersect(evenunderten,primesunderten)
Set{Int64} with 1 element:
  2

julia> nonprimeevens = setdiff(evenunderten,primesunderten)
Set{Int64} with 3 elements:
  4
  6
  8

Get the size of an existing Set

julia> length(numbersunderten)
9

Checking if a set contains a specific element

julia> 16 in numbersunderten
false

julia> 6 in numbersunderten
true

Adding an element to an existing Set

julia> push!(evenunderten,0)
Set{Int64} with 5 elements:
  0
  4
  6
  2
  8

julia> evenunderten
Set{Int64} with 5 elements:
  0
  4
  6
  2
  8

Deleting a specific element from an existing Set

julia> delete!(evenunderten,0)
Set{Int64} with 4 elements:
  4
  6
  2
  8

julia> evenunderten
Set{Int64} with 4 elements:
  4
  6
  2
  8

Get all elements of a Set

@inline function elementsofset(S::Set)
    return [ element for element in S ]
end

julia> elementsofset(evenunderten)
4-element Vector{Int64}:
 4
 6
 2
 8

Pop out the smallest element of a set

function popsmallest!(S::Set)
    if length(S) == 0
        return nothing
    end
    local smallest = minimum(elementsofset(S))
    delete!(S,smallest)
    return smallest
end

julia> popsmallest!(evenunderten)
2

julia> evenunderten
Set{Int64} with 3 elements:
  4
  6
  8

Find all even elements of a set

julia> filter(iseven,elementsofset(numbersunderten))
4-element Vector{Int64}:
 4
 6
 2
 8

Find all odd elements of a set

julia> filter(!iseven,elementsofset(numbersunderten))
5-element Vector{Int64}:
 5
 7
 9
 3
 1

In Statistics you can do sampling without replacement, here we can do it with Sets

julia> function samplingwithoutreplacement(S::Set)
           local sample = rand(S)
           delete!(S,sample)
           return sample
       end
samplingwithoutreplacement (generic function with 1 method)

julia> a = samplingwithoutreplacement(evenunderten)
4

julia> b = samplingwithoutreplacement(evenunderten)
6

julia> evenunderten
Set{Int64} with 1 element:
  8

We can use this to deal out 4 cards out of a pack of 52 playing cards


julia> playingcardset = Set(  reduce(vcat, [ "$(k) of $(suite)" for k = 1:13 , suite = ["Diamond","Heart","Spade","Club"] ] )  )
Set{String} with 52 elements:
  "2 of Club"
  "1 of Diamond"
  "5 of Club"
  "6 of Spade"
  "4 of Club"
  "7 of Heart"
  "2 of Diamond"
  "8 of Diamond"
  "10 of Diamond"
  "2 of Spade"
  "11 of Club"
  "13 of Diamond"
  "5 of Diamond"
  ⋮ 

julia> firstcard = samplingwithoutreplacement(playingcardset)
"1 of Diamond"

julia> secondcard = samplingwithoutreplacement(playingcardset)
"6 of Heart"

julia> thirdcard = samplingwithoutreplacement(playingcardset)
"6 of Spade"

julia> fourthcard = samplingwithoutreplacement(playingcardset)
"4 of Spade"

Now we can simulate Adam and Eve playing poker.

julia> adam = Set([])
Set{Any}()

julia> eve = Set([])
Set{Any}()

julia> playingcardset = Set(  reduce(vcat, [ "$(k) of $(suite)" for k = 1:13 , suite = ["Diamond","Heart","Spade","Club"] ] )  )
Set{String} with 52 elements:
  "2 of Club"
  "1 of Diamond"
  "5 of Club"
  "6 of Spade"
  "4 of Club"
  "7 of Heart"
  "2 of Diamond"
  "8 of Diamond"
  "10 of Diamond"
  "2 of Spade"
  "11 of Club"
  "13 of Diamond"
  "5 of Diamond"
  ⋮ 

julia> for card = 1:5
         push!(adam,samplingwithoutreplacement(playingcardset))
         push!(eve,samplingwithoutreplacement(playingcardset))
       end

julia> adam
Set{Any} with 5 elements:
  "6 of Club"
  "9 of Heart"
  "5 of Heart"
  "12 of Club"
  "9 of Club"

julia> eve
Set{Any} with 5 elements:
  "2 of Spade"
  "10 of Spade"
  "6 of Heart"
  "1 of Club"
  "7 of Heart"

Understanding Markov Chain Monte Carlo with Metropolis Algorithm

Steven Siew — Mon, 10 Apr 2023 01:03:33 +0000

What is Markov Chain Monte Carlo

Markov Chain Monte Carlo is a method by which (additional) samples can be generated (from the last sample) such that the probability density of samples (in total) is proportional to a known function.

What Markov Chain Monte Carlo is used for is parameter estimation (such as means, variances, expected values) and exploration of the posterior distribution of Bayesian models.

Another way to think of Markov Chain Monte Carlo is this:

Imagine a process which produces data, for example, a bias coin producing a series of coin toss results.

H,T,T,H,T,T,T,T,H,T

What Markov Chain Monte Carlo does is to take the result of the coin toss and produce a series of parameter (or parameters) which could be the ACTUAL parameter of the bias coin. The more likely the parameter to be the ACTUAL PARAMETER, the more often it is produced by the Markov Chain Monte Carlo simulation.

So what you need to know is this:

ACTUAL PARAMETER in a process produces the data (a sequence of datum) in an experiment

We collect the actual data and feed it to Markov Chain Monte Carlo Simulation so that it produces a list of Parameter candidates where the most common parameter candidate is the one to be most likely to be the ACTUAL PARAMETER.

What is Metropolis Algorithm

In 1952, Arianna Rosenbluth together with her husband Marshall Rosenbluth worked on what is now known as the Metropolis Hastings algorithm.

They worked on (aka used) the supercomputer MANIAC on Los Alamos National Laboratory in 1952.

Metropolis Hastings algorithm is the very first MCMC algorithm.

https://youtu.be/evyKPVS_M7E
She's a MANIAC, MANIAC on the computing floor
And she's computing like she's never computed before

If you need to read a paper on Metropolis Hastings algorithm, you can get the PDF paper here

https://link.springer.com/article/10.3758/s13423-016-1015-8

I am sorry that I ask you to read the paper above but they can explain MCMC (Markov Chain Monte Carlo) better than I can.

MCMC is a form of Bayesian inference.

Bayesian inference uses the information provided by observed data about a (set of) parameter(s), formally the likelihood, to update a prior state of beliefs about a (set of) parameter(s) to become a posterior state of beliefs about a (set of) parameter(s).

If the above sounds like gobbledygook to you then as a newbie you are NOT ALONE.

What it means is that

You have an idea about the process that generates results in an experiment.
The process have a set of ACTUAL parameters with ACTUAL (but unknown) values.
You have a vague idea about the values of the parameters BEFORE you start the experiment. This idea is called "prior parameter values" or simply as "prior".
You perform the experiment (or ask someone else to do so)
After the experiment, you feed the data from the experiment to MCMC
Then you have to tell MCMC how to calculate the likelihood of obtaining a specific datum from a specific set of parameter values. This is done by providing MCMC with the likelihood function.
The MCMC tool returns the result which is the most likely values of the parameters. This result is called "posterior parameter values" or simply as "posterior"

So, newbie, what this means is that

For MCMC to tell you what the parameter values are, you must first tell MCMC what you THINK the parameter values are BEFORE you start the experiment. Yes, I know it sounds CRAZY!!!

Because the only reason you perform the experiment is to find out the value of the parameters in the first place. What if you know NOTHING about the value of the parameter? What do you tell MCMC? "Hello MCMC, I know nothing about the possible value of the parameter. Let me perform the experiment first, then I can tell you what the values of the parameter might be." But NO. MCMC wants you to tell it what you know about the value of the parameter BEFORE you start the experiment.

Welcome to the crazy world of Bayesian Inference.

Watch a video about Markov Chain Monte Carlo

Please watch this video about Markov Chain Monte Carlo

https://youtu.be/Qqz5AJjyugM

Metropolis Hastings Algorithm

Step 1 : Gather the data from an experiment.

An experiment can be anything. It can be as simple as asking all the students in a high school class to step on a weight scale and record down the weight of each student in kilograms.

data[1] = 55.8 kg
.
.
.
data[10] = 63.5 kg

Step 2 : Come up with a model process which produces the data

The Markov Chain Monte Carlo is NOT a magic tool. It cannot give you a model process. Instead you have to come up with the model process by yourself. You have to use your knowledge, the current scientific knowledge or even pure dumb guesses to come up with a process model.

Say your model process is that

The weight of each student is generated by a Gaussian/Normal Distribution with the mean μ (aka mu) and the standard deviation σ (aka sigma)

In Turing.jl it would look something like this

data[k] ~ Normal(μ,σ)

Step 3 : Pick a random starting value for the parameter

In this step, we need to inform MCMC of our belief about the vague values of our parameters μ (mu) and σ (sigma) before the EXPERIMENT starts.

Here for μ (mu), we need to express our real world knowledge into a mathematical distribution.

We know the lightest person on earth (via www.google.com) is about 6.0 kg and the heaviest person on earth is about 600.0 kg (,remember I said VAGUE values). So we can say

Our belief about the average weight μ (mu) is represented by

prior_μ_distr = Distributions.Uniform(6.0,600.0)

But our daily common sense says most people weights about 60.0 kg so we can also express it as:

prior_μ_distr = Distributions.Normal(60, 10)

And for σ (sigma), we know that it cannot be NEGATIVE. So we can use the Exponential distribution because it is always positive.

prior_μ_distr = Distributions.Normal(60, 10)
prior_σ_distr = Distributions.Exponential(20)
p_μ[1] = rand(prior_μ_distr)
p_σ[1] = rand(prior_σ_distr)

or you can pick any value out of the thin air

prior_μ_distr = Distributions.Uniform(30,90)
prior_σ_distr = Distributions.Uniform(0,50)
p_μ[1] = 40.0
p_σ[1] = 5.0

Step 4 : Pick a range for the parameter(s)

One of the reason, we need Upper Bound and Lower Bound (and JumpingWidth) is to generate the standard deviation for jumping distance.

standard_deviation_for_jumping_distance is
JumpingWidth*(μ_ub-μ_lb)

Later in the Module SimpleMCMC, we don't need Upper Bound and Lower Bound because we can automatically generate a standard deviation to determine the jumping distance.

μ_lb = 30       # Lower bound
μ_ub = 90       # Upper bound

σ_lb = 0        # Lower bound
σ_ub = 50       # Upper bound

Step 5 : get a candidate for new parameter

p_μ_current = p_μ[1]
p_σ_current = p_σ[1]
JumpingWidth = 0.01
p_μ_new = rand(  Normal( p_μ_current , JumpingWidth*(μ_ub-μ_lb) )  )
p_σ_new = rand(  Normal( p_σ_current , JumpingWidth*(σ_ub-σ_lb) )  )

Step 6 : Create a function to determine the probability that the data collected is a result of a particular parameter

# The likelihood function which returns the probability,
# given a fixed mu and a fixed sigma, of generating 
# the sequence of datapoints
function likelihood(data::Vector,mu::Real,sigma::Real)
    numofdatapoints = length(data)
    result=1.0
    for k = 1:numofdatapoints
        result *= pdf(Normal(mu, sigma),data[k])
    end
    return result
end

where

? pdf

pdf(d::Distribution{ArrayLikeVariate{N}}, x) where {N}

Evaluate the probability density function of d at every element in a collection x.

Step 7 : Starting at parameter[1] , generate a fix number of new parameters using the following algorithm:

Generate a candidate.
Find the prior for the current_parameter. (priorμ0 and priorσ0)
Find the prior for the candidate_for_new_parameter. (priorμ1 and priorσ1)
Find the likelihood of current_parameter.
Find the likelihood of candidate_for_new_parameter.
Calculate q0 and q1.
If q1 > q0 then select the candidate_for_new_parameter. Otherwise generate a random number between 0 and 1.
If the random number is less than q1/q0 then select the candidate_for_new_parameter otherwise select current_parameter

n_samples = 1_000_000

for i in 2:n_samples
    # Generate a candidate. 
    p_μ_new = rand(  Normal( p_μ[i-1] , JumpingWidth*(μ_ub-μ_lb) )  )
    p_σ_new = rand(  Normal( p_σ[i-1] , JumpingWidth*(σ_ub-σ_lb) )  )

    priorμ0 = Distributions.pdf(prior_μ_distr,p_μ[i-1])
    priorσ0 = Distributions.pdf(prior_σ_distr,p_μ[i-1])    

    priorμ1 = Distributions.pdf(prior_μ_distr,p_μ_new)
    priorσ1 = Distributions.pdf(prior_σ_distr,p_σ_new)

    # Find the likelihood of current_parameter
    L0 = likelihood(data, p_μ[i-1], p_σ[i-1])

    # Calculate q0
    q0 = L0 * priorμ0 * priorσ0

    # Find the likelihood of candidate_for_new_parameter.
    L1 = likelihood(data, p_μ_new, p_σ_new)

    # Calculate q1
    q1 =  L1 * priorμ1 * priorσ1

    if rand() < q1/q0  # Then jump to the NEW POSITION
        p_μ[i] = p_μ_new
        p_σ[i] = p_σ_new
    else  # Stick with the CURRENT POSITION
        p_μ[i] = p_μ[i-1]
        p_σ[i] = p_σ[i-1]
    end
end

Step 8 : Now we have the array of parameters. We can plot the density graphs.

array of p_μ
array of p_σ

# Next we need to display the histogram of the array p
# this is the shape of the posterior produced by MCMC
histogram(p_μ,legend=false,title="Posterior of mu") |> display
histogram(p_σ,legend=false,title="Posterior of sigma") |> display

Here is a diagram from a youtube video that illustrate what we are trying to do.

Simple Julia implementation of Metropolis Algorithm

using Distributions
using StatsPlots

# Here are the weights of the students (in kilograms) in a classrom
# Assuming the weight of the students are generated by a biological
# process with a normal distribution with mean μ and standard dev σ
# find the posterior distribution of those two parameters
data = [ 55.8, 64.1, 61.2, 55.2, 58.0, 69.8, 53.2, 70.0, 63.3, 63.5 ]

# define our prior for μ as a normal centered on 60 ( required for MCMC )
prior_μ_distr = Distributions.Normal(60, 10)
prior_μ(p) = Distributions.pdf(prior_μ_distr,p)
μ_lb = 30       # Lower bound
μ_ub = 90       # Upper bound

# define our prior for σ as a exponential with mean of 20 ( required for MCMC )
prior_σ_distr = Distributions.Exponential(20)
prior_σ(p) = Distributions.pdf(prior_σ_distr,p)
σ_lb = 0        # Lower bound
σ_ub = 50       # Upper bound

# define the number of samples used in MCMC
# And create array parameter with that many elements (filled with 0.0)
n_samples = 1_000_000
p_μ = zeros(n_samples)
p_σ = zeros(n_samples)

# Set the starting value to a random number between lb and ub
p_μ[1] = rand(prior_μ_distr)
while p_μ[1] < μ_lb || p_μ[1] > μ_ub
    p_μ[1] = rand(prior_μ_distr)
end

p_σ[1] = rand(prior_σ_distr)
while p_σ[1] < σ_lb || p_σ[1] > σ_ub
    p_σ[1] = rand(prior_σ_distr)
end

# The likelihood function which returns the probability,
# given a fixed mu and a fixed sigma, of the sequence of datapoints
function likelihood(data::Vector,mu::Real,sigma::Real)
    numofdatapoints = length(data)
    result=1.0
    for k = 1:numofdatapoints
        result *= pdf(Normal(mu, sigma),data[k])
    end
    return result
end

JumpingWidth = 0.01
for i in 2:n_samples
    # p_new has a value around the vicinity of p[i-1]
    p_μ_new = rand(  Normal( p_μ[i-1] , JumpingWidth*(μ_ub-μ_lb) )  )
    # make sure p_μ_new is between μ_lb and μ_ub
    if p_μ_new < μ_lb
        p_μ_new = μ_lb + abs( p_μ_new - μ_lb )
    elseif p_μ_new > μ_ub
        p_μ_new = μ_ub - abs( p_μ_new - μ_ub )
    end

    p_σ_new = rand(  Normal( p_σ[i-1] , JumpingWidth*(σ_ub-σ_lb) )  )
    # make sure p_σ_new is between σ_lb  and σ_ub
    if p_σ_new < σ_lb
        p_σ_new = σ_lb + abs( p_σ_new - σ_lb )
    elseif p_σ_new > σ_ub
        p_σ_new = σ_ub - abs( p_σ_new - σ_ub )
    end
    " q0 is posterior0 = likelihood0 * prior0        "
    " q1 is posterior1 = likelihood1 * prior1        "
#=
q0  is the probability for the CURRENT POSITION
q1  is the probability for the CANDIDATE for the NEW POSITION
=#
    priorμ0 = prior_μ(p_μ[i-1])
    priorσ0 = prior_σ(p_σ[i-1])
    priorμ1 = prior_μ(p_μ_new)
    priorσ1 = prior_σ(p_σ_new)

    q0 = likelihood(data, p_μ[i-1], p_σ[i-1]) * priorμ0 * priorσ0
    q1 = likelihood(data, p_μ_new, p_σ_new) * priorμ1 * priorσ1
    # The value of parameter[i] depends on whether the
    # random number is less than q1/q0
#=
If  the probability for the CANDIDATE for the NEW POSITION (q1)
    is GREATER than
the probability for the CURRENT POSITION (q0)
then jump to the NEW POSITION

Otherwise

Only jump to the CANDIDATE for the NEW POSITION if
a random number (between 0.0 and 1.0) is less than   q1/q0
=#
    if rand() < q1/q0  # Then jump to the NEW POSITION
        p_μ[i] = p_μ_new
        p_σ[i] = p_σ_new
    else  # Stick with the CURRENT POSITION
        p_μ[i] = p_μ[i-1]
        p_σ[i] = p_σ[i-1]
    end


end

# Next we need to display the histogram of the array p
# this is the shape of the posterior produced by MCMC
histogram(p_μ,legend=false,title="Posterior of mu") |> display
histogram(p_σ,legend=false,title="Posterior of sigma") |> display

This shows that the ACTUAL PARAMETER value for mu is around 61.4

This shows that the ACTUAL PARAMETER value for mu is around 6.66

A optimized version of Metropolis Algorithm using Log function and burned in samples

Here we enhanced the algorithm by using two further optimization

Use logarithm addition instead of standard arithmetic multiplication

# The likelihood function which returns the probability,
# given a fixed mu and a fixed sigma, of the sequence of datapoints
function likelihood(data::Vector,mu::Real,sigma::Real)
    numofdatapoints = length(data)
    result=1.0
    for k = 1:numofdatapoints
        result *= pdf(Normal(mu, sigma),data[k])
    end
    return result
end

Because we are doing lots and lots and lots of multiplication, even Float64 can run out of room. Therefore most MCMC programs uses Logarithm and this ensures that Float64 will not hit its limitation.

So we do

# The likelihood function which returns the probability,
# given a fixed mu and a fixed sigma, of the sequence of datapoints
function loglikelihood(data::Vector,mu::Real,sigma::Real)
    numofdatapoints = length(data)
    result=0.0
    for k = 1:numofdatapoints
        result += log(  pdf(Normal(mu, sigma),data[k])  )
    end
    return result
end

Dispose of the burned_in_samples

We should get rid of the initial burned in samples because they are quite far away from the expected proper probability distribution of the parameter

# define the number of samples used in MCMC
# And create array parameter with that many elements (filled with 0.0)
burnedinsamples = 1000
n_samples = 1_000_000 + burnedinsamples

#= Lots of code =#

# Finally we must not forget to remove the burned in samples
deleteat!(p_μ,1:burnedinsamples)
deleteat!(p_σ,1:burnedinsamples)

So here is the complete source code

using Distributions
using StatsPlots

# Here are the weights of the students (in kilograms) in a classrom
# Assuming the weight of the students are generated by a biological
# process with a normal distribution with mean μ and standard dev σ
# find the posterior distribution of those two parameters
data = [ 55.8, 64.1, 61.2, 55.2, 58.0, 69.8, 53.2, 70.0, 63.3, 63.5 ]

# define our prior for μ as a normal centered on 60 ( required for MCMC )
prior_μ_distr = Distributions.Normal(60, 10)
prior_μ(p) = Distributions.pdf(prior_μ_distr,p)
μ_lb = 30       # Lower bound
μ_ub = 90       # Upper bound

# define our prior for σ as a exponential with mean of 20 ( required for MCMC )
prior_σ_distr = Distributions.Exponential(20)
prior_σ(p) = Distributions.pdf(prior_σ_distr,p)
σ_lb = 0        # Lower bound
σ_ub = 50       # Upper bound

# define the number of samples used in MCMC
# And create array parameter with that many elements (filled with 0.0)
burnedinsamples = 1000
n_samples = 1_000_000 + burnedinsamples
p_μ = zeros(n_samples)
p_σ = zeros(n_samples)

# Set the starting value to a random number between lb and ub
p_μ[1] = rand(prior_μ_distr)
while p_μ[1] < μ_lb || p_μ[1] > μ_ub
    p_μ[1] = rand(prior_μ_distr)
end

p_σ[1] = rand(prior_σ_distr)
while p_σ[1] < σ_lb || p_σ[1] > σ_ub
    p_σ[1] = rand(prior_σ_distr)
end


# The likelihood function which returns the probability,
# given a fixed mu and a fixed sigma, of the sequence of datapoints
function loglikelihood(data::Vector,mu::Real,sigma::Real)
    numofdatapoints = length(data)
    result=0.0
    for k = 1:numofdatapoints
        result += log(  pdf(Normal(mu, sigma),data[k])  )
    end
    return result
end

JumpingWidth = 0.01
for i in 2:n_samples
    # p_new has a value around the vicinity of p[i-1]
    p_μ_new = rand(  Normal( p_μ[i-1] , JumpingWidth*(μ_ub-μ_lb) )  )
    # make sure p_μ_new is between μ_lb and μ_ub
    if p_μ_new < μ_lb
        p_μ_new = μ_lb + abs( p_μ_new - μ_lb )
    elseif p_μ_new > μ_ub
        p_μ_new = μ_ub - abs( p_μ_new - μ_ub )
    end

    p_σ_new = rand(  Normal( p_σ[i-1] , JumpingWidth*(σ_ub-σ_lb) )  )
    # make sure p_σ_new is between σ_lb  and σ_ub
    if p_σ_new < σ_lb
        p_σ_new = σ_lb + abs( p_σ_new - σ_lb )
    elseif p_σ_new > σ_ub
        p_σ_new = σ_ub - abs( p_σ_new - σ_ub )
    end
    " q0 is posterior0 = likelihood0 * prior0        "
    " q1 is posterior1 = likelihood1 * prior1        "
#=
q0  is the probability for the CURRENT POSITION
q1  is the probability for the CANDIDATE for the NEW POSITION
=#
    logpriorμ0 = log(  prior_μ(p_μ[i-1])  )
    logpriorσ0 = log(  prior_σ(p_σ[i-1])  )
    logpriorμ1 = log(  prior_μ(p_μ_new)   )
    logpriorσ1 = log(  prior_σ(p_σ_new)   )

    logq0 = loglikelihood(data, p_μ[i-1], p_σ[i-1]) + logpriorμ0 + logpriorσ0
    logq1 = loglikelihood(data, p_μ_new, p_σ_new) + logpriorμ1 + logpriorσ1
    # The value of parameter[i] depends on whether the
    # log(random number) is less than logq1 - logq0
#=
If  the probability for the CANDIDATE for the NEW POSITION (q1)
    is GREATER than
the probability for the CURRENT POSITION (q0)
then jump to the NEW POSITION

Otherwise

Only jump to the CANDIDATE for the NEW POSITION if
a random number (between 0.0 and 1.0) is less than   q1/q0
=#
    if log(rand()) < logq1 - logq0    # Then jump to the NEW POSITION
        p_μ[i] = p_μ_new
        p_σ[i] = p_σ_new
    else        # Stick with the CURRENT POSITION
        p_μ[i] = p_μ[i-1]
        p_σ[i] = p_σ[i-1]
    end


end

# Finally we must not forget to remove the burned in samples
deleteat!(p_μ,1:burnedinsamples)
deleteat!(p_σ,1:burnedinsamples)

# Next we need to display the histogram of the array p
# this is the shape of the posterior produced by MCMC
histogram(p_μ,legend=false,title="Posterior of mu") |> display
histogram(p_σ,legend=false,title="Posterior of sigma") |> display

How to make this runs faster and easier to use

According to performance tricks for Julia

Any code that is performance critical should be inside a function. Code inside functions tends to run much faster than top level code, due to how Julia's compiler works.

The second problem is that each time we want to perform MCMC on a set of rawdata, we need to write the entire program again in Julia Programming Language.

A much better solution is to write it as a module and call the module to use MCMC.

So this version of MCMC is written as a module and it requires the following packages

Distributions
StatsPlots
Statistics
Printf
PrettyTables
DoubleFloats
CSV
DataFrames
Random

The best thing to do is to split up the code in two separate parts.

The SimpleMCMC module in a file called "SimpleMCMC.jl"
The main program in a file called "MCMC_example01.jl"
Put these two files in their own directory. Which in my case is the directory /Users/ssiew/juliascript/MCMC_dir on my Macbook
Add code to SimpleMCMC modules to save the graphics to the current directory. This would make it simple to have a copy of the graphs if we need to look at them at a later time.
Use the directory for our own enviroment when running julia.

Notes: vi is a program on unix/MacOS to edit a textfile

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ vi SimpleMCMC.jl

# Filename is "SimpleMCMC.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"
using Pkg
Pkg.add(["Distributions","StatsPlots","Statistics",
"Printf","PrettyTables","DoubleFloats",
"CSV","DataFrames","Random"])

module SimpleMCMC
    export MetropolisHastings_MCMC
    export Gibbs_MCMC
    export describe_paramvec,describe_multiparamvec
    export find_best_param
    export vec2pdfcdf, find_invcdf, find_cdf, find_pdf
    export BayesFactorStrengthToString
    export findHPDI, ComparePost2Prior
    export CreateTrainTest

    using Distributions, StatsPlots, Statistics, Printf, PrettyTables
    using DoubleFloats, CSV, DataFrames, Random

    # Version 7
    #  v7: Add find_best_param, vec2pdfcdf, find_invcdf, find_cdf, find_pdf
    #  v6: Add Gibbs_MCMC and change module name from MHMCMC to SimpleMCMC
    #  v5: turn it into a module
    #  v4: remove prior from Parameter struct and the struct itself
    #  v3: create JumpingWidthVec, lbVec, ubVec  to improve efficiency
    #  v2: Remove q0 from the main loop to improve efficiency
    #  v1: Turn it into a version that uses pure function(s)

    struct Bin
        lb::Double64
        ub::Double64
        center::Double64
    end

    function MetropolisHastings_MCMC(data,lop::Vector{Distribution},loglikelihood;samples,burnedinsamples,JumpingWidth=0.01)
        # lop is "List Of Parameters"
        # Find the number of parameters
        numofparam = length(lop)
        # calc the n_samples
        n_samples = samples + burnedinsamples
        # Create a parameter Matrix p with n_samples rows and numofparam cols
        p = zeros(n_samples,numofparam)
        # JumpingWidthVector, lowerboundVector and upperboundbVector
        JumpingWidthVec = zeros(numofparam)
        lbVec = zeros(numofparam)
        ubVec = zeros(numofparam)
        # Create a Vector of our logprior function
        logpriorVec = Array{Function}(undef,numofparam)

        # Set the starting value to a random number between lb and ub
        for k = 1:numofparam  # for each parameter k
            # Set the starting value to a random number between lb and ub
            p[1,k] = rand(lop[k])

            # get 100 random values and find the practical max and practical min
            HundredRandValues = rand(lop[k],100)
            prac_min = minimum(HundredRandValues)
            prac_max = maximum(HundredRandValues)
            JumpingWidthVec[k] = JumpingWidth * (prac_max - prac_min)

            lbVec[k] = minimum(lop[k])
            ubVec[k] = maximum(lop[k])

            # Array of function prior(x) = pdf(distribution,x)
            logpriorVec[k] = x->logpdf(lop[k],x)
        end

        p_old = vec(p[1,:])
        q0 = loglikelihood(data,p_old...) + sum([logpriorVec[k](p_old[k]) for k = 1:numofparam ])

        # prepare the p_new array
        p_new = zeros(numofparam)
        # Main loop for MetropolisHastings MCMC
        @inbounds for i = 2:n_samples
            # p_new has a value around the vicinity of p[i-1]
            for k = 1:numofparam  # for each parameter k
                # p_new has a value around the vicinity of p[i-1]
                p_new[k] = rand(  Normal( p[i-1,k] , JumpingWidthVec[k] )  )
                # make sure p_new is between lb and ub
                if p_new[k] < lbVec[k]
                    p_new[k] = lbVec[k] + abs( p_new[k] - lbVec[k] )
                elseif p_new[k] > ubVec[k]
                    p_new[k] = ubVec[k] - abs( p_new[k] - ubVec[k] )
                end
            end

            # Calc the two posterior
            " q0 is posterior0 = likelihood0 * prior0        "
            " q1 is posterior1 = likelihood1 * prior1        "
            q1 = loglikelihood(data,p_new...) + sum([logpriorVec[k](p_new[k]) for k = 1:numofparam ])
            # The value of p[i] depends on whether the
            # random number is less than q1/q0
            p[i,:] .= log(rand()) < q1-q0 ? (q0=q1; p_old=p_new[:]) : p_old
        end
        # Finally we must not forget to remove the burned in samples
        return p[(burnedinsamples+1):end,:]
    end

    function Gibbs_MCMC(data,lop::Vector{Distribution},loglikelihood;samples,burnedinsamples)
        # lop is "List Of Parameters"
        # Find the number of parameters
        numofparam = length(lop)
        # calc the n_samples
        n_samples = samples + burnedinsamples
        # Create a parameter Matrix p with n_samples rows and numofparam cols
        p = zeros(n_samples,numofparam)
        # Number of holes to drill
        numofdrillholes = 1001
        drillhole = Array{Float64}(undef,numofdrillholes)
        mycdf = Array{Double64}(undef,numofdrillholes+1)
        # prepare the vectors: startposvec stopposvec chunksizevec
        startposvec = zeros(numofparam)
        stopposvec = zeros(numofparam)
        chunksizevec = zeros(numofparam)

        # Set the starting value to a random number between lb and ub
        for k = 1:numofparam  # for each parameter k
            # Set the starting value to a random number between lb and ub
            p[1,k] = rand(lop[k])

            # get 400 random values and find the practical max and practical min
            FourHundredRandValues = rand(lop[k],400)
            startposvec[k] = minimum(FourHundredRandValues)
            stopposvec[k]  = maximum(FourHundredRandValues)
            chunksizevec[k] = (stopposvec[k]-startposvec[k])/(numofdrillholes-1)
        end

        # Main loop for Gibbs MCMC
        @inbounds for i = 2:n_samples
            # prepare the p_new vector
            p_new = vec(p[i-1,:])
            # Select a new value for each parameter of p_new
            for k = 1:numofparam  # for each parameter k
                # Calculate chunk size
                startpos,stoppos,chunksize = startposvec[k],stopposvec[k],chunksizevec[k]
                mycdf[1] = Double64(0.0)
                # Start drilling the drill holes for parameter k
                for n = 1:numofdrillholes
                    p_new[k] = startpos + (n-1) * chunksize
                    drillhole[n] = exp(  loglikelihood(data,p_new...)  )
                    mycdf[n+1] = mycdf[n] + drillhole[n]
                end
                mycdf /= mycdf[end]
                # Now we perform an inverse CDF sampling
                x = 0.0
                while x == 0.0
                    x = rand()  # Make sure x is not zero
                end
                counter = 2
                while !(mycdf[counter-1] < x <= mycdf[counter])
                    counter += 1
                end
                counter -= 1
                # now find the kth parameter value
                pos = startpos + (counter-1) * chunksize
                # This is the value we are going to jump to
                p_new[k] = pos
                p[i,k] = pos
            end
        end
        # Finally we must not forget to remove the burned in samples
        return p[(burnedinsamples+1):end,:]
    end

    function describe_paramvec(name,namestr,v::Vector)
        m,s,v_sorted = mean(v),std(v),sort(v)
        pc5 = v_sorted[Int64(round(0.055*length(v_sorted)))]
        pc94 = v_sorted[Int64(round(0.945*length(v_sorted)))]
        println(
            "Parameter $(name)   mean ",@sprintf("%6.3f",round(m,sigdigits=3)),
            "   sd ",@sprintf("%6.3f",round(s,sigdigits=3)),
            "   5.5% ",@sprintf("%6.3f",round(pc5,sigdigits=3)),
            "   94.5% ",@sprintf("%6.3f",round(pc94,sigdigits=3))
        )
        histogram(v,legend=false,title="Posterior of $(namestr)") |> display
        return (m,s,pc5,pc94)
    end

    function describe_multiparamvec(multi::Vector;savegraph=false)
        # name,namestr,v::Vector
        numofparams = length(multi)
        table = Array{Any}(undef,numofparams,5)
        for i = 1:numofparams
            v = multi[i][3]
            m,s,v_sorted = mean(v),std(v),sort(v)
            pc5 = v_sorted[Int64(round(0.055*length(v_sorted)))]
            pc94 = v_sorted[Int64(round(0.945*length(v_sorted)))]
            table[i,1]=multi[i][1]
            table[i,2],table[i,3],table[i,4],table[i,5]=m,s,pc5,pc94
            histogram(v,legend=false,title="Posterior of $(multi[i][2])") |> display
            # Save the graphics to the current directory
            if savegraph
                savefig("Posterior of $(multi[i][2]).png")
            end
        end
        pretty_table(table;
            header = tuple(["Parameter","Mean","Stddev","5.5%","94.5%"]),
            alignment=:l,
            #formatter=Dict(0 => (v,i) -> typeof(v) <: Number ? round(v,sigdigits=3) : v )
            formatters=((v,i,j) -> typeof(v) <: Number ? round(v,sigdigits=3) : v,)
        )
        return table
    end

    function find_best_param(data,param_mat,loglikelihood)
        bestpos = 0
        bestscore = -Inf
        for k = 1:size(param_mat,1)
            currscore = loglikelihood(data,param_mat[k,:]...)
            if currscore > bestscore
                bestscore = currscore
                bestpos = k
            end
        end
        return (bestpos,bestscore)
    end

    function vec2pdfcdf(vec::Vector,nbins::Int64=20)
        vec_sorted = sort(vec)
        len = length(vec)
        pc0001_pos = Int64(round(len * 0.0001,RoundDown))
        pc0001_pos = pc0001_pos > 0 ? pc0001_pos : 1
        pc9999_pos = Int64(round(len * 0.9999,RoundDown))
        pc9999_pos = pc9999_pos < len ? pc9999_pos : len

        chunksize = (Double64(vec_sorted[pc9999_pos]) -
                     Double64(vec_sorted[pc0001_pos])) / Double64(nbins-2)
        bins = Vector{Bin}(undef,nbins)
        binscount = Int64[ 0 for k = 1:nbins ]
        bins[1]   = Bin(Double64(vec_sorted[pc0001_pos]) - chunksize,
                        Double64(vec_sorted[pc0001_pos]),
                        Double64(vec_sorted[pc0001_pos]) - (chunksize/df64"2.0"))
        bins[end] = Bin(Double64(vec_sorted[pc9999_pos]),
                        Double64(vec_sorted[pc9999_pos]) + chunksize,
                        Double64(vec_sorted[pc9999_pos]) + (chunksize/df64"2.0"))
        # Create and specify each bin
        for k = 2:nbins-1
            bins[k] = Bin(bins[1].ub+(k-2)*chunksize,
                          bins[1].ub+(k-1)*chunksize,
                          bins[1].ub+(k-1.5)*chunksize)
        end
        # Now perform the histogram function
        for n = 1:len
            ((vec[n] < bins[1].lb) || (vec[n] >= bins[end].ub)) && continue
            k=1
            while !(bins[k].lb <= vec[n] < bins[k].ub)
                k +=1
            end
            binscount[k] += 1
        end
        binscount_total = sum(binscount)
        # Now perform normalization
        pdfvalue = Float64[ Float64(binscount[k]/binscount_total) for k = 1:nbins ]
        # println("Sum of binsvalue is ",sum(binsvalue))
        centervec = Float64[ Float64(bins[k].center) for k = 1:nbins ]
        # Now create cdf
        cdfcount = Array{Int64}(undef,nbins)
        cdfcount[1] = binscount[1]
        for k = 2:nbins
            cdfcount[k] = cdfcount[k-1] + binscount[k]
        end
        cdfvalue = Float64[ Float64(cdfcount[k]/binscount_total) for k = 1:nbins ]
        return (centervec,pdfvalue,cdfvalue)
    end

    function find_invcdf(gridvec,cdfvec,x)
        if x < cdfvec[1]
            return gridvec[1]
        end
        counter = 2
        while !(cdfvec[counter-1] < x <= cdfvec[counter])
            counter += 1
        end
        return gridvec[counter]
    end

    function find_cdf(gridvec,cdfvec,x)
        len = length(gridvec)
        width = gridvec[2] - gridvec[1]
        # estimate counter
        startpos = gridvec[1] - 1.5 * width
        k = Int64(round((x - startpos)/width,RoundDown))
        if k < 1
            return 0.0
        elseif  k > len
            return 1.0
        else
            return cdfvec[k]
        end
    end

    function find_pdf(gridvec,pdfvec,x)
        len = length(gridvec)
        width = gridvec[2] - gridvec[1]
        # estimate counter
        startpos = gridvec[1] - 1.5 * width
        k = Int64(round((x - startpos)/width,RoundDown))
        if k < 1 || k > len
            return 0.0
        else
            return pdfvec[k]
        end
    end

    function BayesFactorStrengthToString(x)
        if 10.0^0.0 <= x < 10.0^0.5
            return "Barely worth mentioning"
        elseif 10.0^0.5 <= x < 10.0^1.0
            return "Substantial"
        elseif 10.0^1.0 <= x < 10.0^1.5
            return "Strong"
        elseif 10.0^1.5 <= x < 10.0^2.0
            return "Very Strong"
        elseif 10.0^2.0 <= x
            return "Decisive"
        else
            return "AGAINST at " * BayesFactorStrength(1.0/x)
        end
    end

    # Now we create the function findHPDI()
    function findHPDI(grid::Vector,pdf::Vector,mass,depth)
        len_grid,len_pdf = length(grid),length(pdf)
        @assert len_grid == len_pdf "Grid vector and pdf vector does not have the same length"
        result = findHPDI_coarsegrain(depth,mass,1,len_pdf,pdf)
        # We have done the course grain, so now we need to
        # do the fine grain adjustment
        while true
            # Create a list of candidates
            candidates = []
            for startadj = -10:10
                for stopadj = -10:10
                    cand = (  max(result[1]+startadj,1),
                              min(result[2]+stopadj,len_pdf)
                           )
                    push!(candidates,cand)
                end
            end
            # We need to evaluate the candidates by one criterion
            # The criterion is mininum mass
            eligible = []
            for candidate in candidates
                candidate_mass = sum(pdf[candidate[1]:candidate[2]])
                candidate_width = candidate[2]-candidate[1]
                if  candidate_mass >= mass && candidate_width > 0
                    candidate_density = candidate_mass / candidate_width
                    push!(eligible,[candidate_density,candidate])
                    #= println("FG: candidate ",candidate," has mass ",
                    candidate_mass,
                    " and density ",candidate_density)
                    =#
                end
            end
            sort!(eligible)
            # println(eligible)
            # Now extract the best candidate
            best = eligible[end]
            best_density = best[1]
            best_candidate = best[2]
            best_mass = sum(pdf[best_candidate[1]:best_candidate[2]])
            #= println("FG: Best candidate ",best_candidate,
            "has mass ",best_mass,
            " has density ",best_density)
            =#
            # Now call itself recursively
            if best_candidate[1] == result[1] && best_candidate[2] == result[2]
                # There is no point recursing down any further as
                # we reached the limit of recursion
                return result
            else
                # Subtract the depth level and recurse down further
                result = best_candidate
            end
        end
    end

    function findHPDI_coarsegrain(depth::Integer,mass,start,stop,pdf)
        @assert start <= stop "start is greater than stop"
        if depth < 1 || stop == start || stop == start + 1
            return (start,stop)
        end
        len = stop - start
        chunk = len ÷ 3
        # Handle special case where chunk is zero
        if chunk == 0
            # This is the limit of Coarse Grain algorithmn
            return (start,stop)
        end
        # There are four points
        point1 = start
        point2 = start + chunk
        point3 = start + 2 * chunk
        point4 = stop
        # From these four points,
        # we can produce 5 candidates
        candidate1 = (point1,point2)
        candidate2 = (point1,point3)
        candidate3 = (point2,point3)
        candidate4 = (point2,point4)
        candidate5 = (point3,point4)
        candidate6 = (point1,point4)
        # Create a list of candidates
        candidates = [ candidate1, candidate2, candidate3, candidate4, candidate5, candidate6 ]
        # We need to evaluate the candidates by one criterion
        # The criterion is mininum mass
        eligible = []
        for candidate in candidates
            candidate_mass = sum(pdf[candidate[1]:candidate[2]])
            candidate_width = candidate[2]-candidate[1]
            candidate_density = candidate_mass / candidate_width
            if  candidate_mass >= mass
                push!(eligible,[candidate_density,candidate])
                #= println("candidate ",candidate," has mass ",
                candidate_mass,
                " and density ",candidate_density)
                =#
            end
        end
        sort!(eligible)
        # println(eligible)
        # Now extract the best candidate
        best = eligible[end]
        best_density = best[1]
        best_candidate = best[2]
        #= println("Best candidate ",best_candidate,
        " has density ",best_density)
        =#
        # Now call itself recursively
        if best_candidate[1] == start && best_candidate[2] == stop
            # There is no point recursing down any further as
            # we reached the limit of recursion
            return (start,stop)
        else
            # Subtract the depth level and recurse down further
            return findHPDI_coarsegrain(depth-1,mass,best_candidate[1],best_candidate[2],pdf)
        end
    end

    function ComparePost2Prior(prior,vec,nbins=1001;kwargs...)
        lb = minimum(prior)
        ub = maximum(prior)
        if lb == -Inf || ub == Inf
            FourHundredRandValues = rand(prior,400)
            lb = lb == -Inf ? minimum(FourHundredRandValues) : lb
            ub = ub ==  Inf ? maximum(FourHundredRandValues) : ub
        end
        prior_range = range(lb,length=1001,stop=ub)
        (centervec,pdfvalue,cdfvalue) = vec2pdfcdf(vec,nbins)
        prior_vec = [ pdf(prior,x) for x in prior_range ]
        pdf_vec = [ find_pdf(centervec,pdfvalue,x) for x in prior_range ]
        scaling_factor = sum(prior_vec)/sum(pdf_vec)
        pdf_vec = pdf_vec * scaling_factor
        plot(prior_range,[prior_vec,pdf_vec],
            xlimit=(lb,ub),legend=false;kwargs...
        ) |> display
    end

    function CreateTrainTest(df::DataFrame,prop=0.5,randomseed=1234)
        df_training = similar(df,0)
        df_testing  = similar(df,0)

        # Now split the df into df_training and df_testing
        df_size = size(df,1)
        training_proportion = prop
        trainingsize = round(df_size*training_proportion)

        # Create a random permutaion vector
        randvec = randperm!(MersenneTwister(randomseed),
                            Vector{Int64}(undef,df_size))

        for k in axes(df)[1]
            push!(  k ≤ trainingsize ?
                    df_training : df_testing ,
                    df[randvec[k],:]
            )
        end
        return (df_training,df_testing)
    end

end
#= End of module SimpleMCMC =#

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ vi MCMC_example01.jl

# Filename is "MCMC_example01.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"

# include the source code for the SimpleMCMC module
include("/Users/ssiew/juliascript/MCMC_dir/SimpleMCMC.jl")

using Distributions, StatsPlots, Statistics
using .SimpleMCMC

function MCMC_example()

    # Here are the weights of the students (in kilograms) in a classrom
    # Assuming the weight of the students are generated by a biological
    # process with a normal distribution with mean μ and standard dev σ
    # find the posterior distribution of those two parameters
    data = [ 55.8, 64.1, 61.2, 55.2, 58.0, 69.8, 53.2, 70.0, 63.3, 63.5 ]

    # define our prior for μ as a normal centered on 60 ( required for MCMC )
    μ = Normal(60, 10)
    # define our prior for σ as a exponential with mean of 20 ( required for MCMC )
    σ = Exponential(20)

    # Put all our parameters into a list
    ListOfParams = Distribution[μ,σ]

    # The loglikelihood function which returns the log(probability),
    # given a fixed mu and a fixed sigma, of the sequence of datapoints
    function loglikelihood(data::Vector,mu,sigma)
        numofdatapoints = length(data)
        result=0.0
        for k = 1:numofdatapoints
            result += logpdf(Normal(mu,sigma),data[k])
        end
        return result
    end

    p = MetropolisHastings_MCMC(
            data,
            ListOfParams,
            loglikelihood,
            samples=1_000_000,
            burnedinsamples=1000
    )

    # Next we need to display the histogram of the column p[:,1]
    # this is the shape of the posterior produced by MCMC
    vec_μ=vec(p[:,1])  # Convert a Matrix column into a 1-D Vector
    vec_σ=vec(p[:,2])  # Convert a Matrix column into a 1-D Vector
    #describe_paramvec("μ","mu",vec_μ)
    #describe_paramvec("σ","sigma",vec_σ)
    describe_multiparamvec([("μ","mu",vec_μ),("σ","sigma",vec_σ)])
    histogram2d(vec_μ,vec_σ,nbins=100,xlabel="mu",ylabel="sigma") |> display
    # save the graphics to the current directory
    savefig("histogram2D_mu_sigma.png")
end

# We call the main function to kick off the example
MCMC_example()

Now we can call Julia on the terminal to run the script

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ julia --project=. MCMC_example01.jl

Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
[ Info: Precompiling SimpleMCMC [top-level]
┌───────────┬──────┬────────┬──────┬───────┐
│ Parameter │ Mean │ Stddev │ 5.5% │ 94.5% │
├───────────┼──────┼────────┼──────┼───────┤
│ μ         │ 61.4 │ 2.09   │ 58.1 │ 64.6  │
│ σ         │ 6.66 │ 1.85   │ 4.44 │ 9.89  │
└───────────┴──────┴────────┴──────┴───────┘

What can MCMC do that ordinary high school statistics can't do?

So far, you have been underwhelm by the MCMC tool. We got a bunch of student weights and found the mean and standard deviation of the students. Big deal! Ordinary high school statistics can do that!

Why go through so much complex programming just to find the mean and standard deviation? It really really felt much ado about nothing. Please show us something that MCMC can do and high school statistics CANNOT DO.

Okay. Consider this problem. There are 100 students in a class. 55% of the students are male and 45% of the students are female. The teacher measured the heights of the students in metres. The teacher DID NOT record down the name or the gender of each student.

Your mission, if you can even do it, is to find the mean and standard deviation of

The male students
The female students

Now this really is mission impossible as the male students and female students are mixed together. You cannot simply tell from a height if the student is male or if the student is female and yet you have to "somehow" seperate them into two gender and then calculate their gender mean and gender standard deviation.

This is truely something that high school statistics cannot do.

Generate 100 datapoints

First thing first, we need to generate 100 datapoints to represent the height of the students. We also need to make sure that we get the same datapoints each time. To do so, we use StableRNGs which guarranty that we get the same random numbers everytime (even if the Julia version changes).

In the code below, we use the following ACTUAL PARAMETERS

μm_actual = 1.76 (male mean)
σm_actual = 0.0635 (male std dev)
μf_actual = 1.62 (female mean)
σf_actual = 0.0559 (female std dev)

These will be the numbers that we can trying to estimate using MCMC

using StableRNGs,Distributions

function generatedata(;numofdatum = 100)
    rng = StableRNG(12345)
    data = zeros(numofdatum)
    # 55% male with 1.76 m and std dev of 0.0635
    for k = 1:Int64(55//100*numofdatum)
        data[k] = rand(rng,Normal(1.76,0.0635))
    end
    # 45% female with 1.62 m and std dev of 0.0559
    for k = (1+Int64(55//100*numofdatum)):numofdatum
        data[k] = rand(rng,Normal(1.62,0.0559))
    end
    return round.(data,digits=2)
end

data = generatedata()

Error Score

As we know what the Actual Parameter values are, we can compute how wrong our estimates are.

We defined a function called errorscore to compute the error score value where the large the error score is, the more wrong the estimated values are when compare to the actual values.

First attempt

# Filename is "bothgender01.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"
include("/Users/ssiew/juliascript/MCMC_dir/SimpleMCMC.jl")

using Distributions, StatsPlots, Statistics
using .SimpleMCMC
using StableRNGs,Distributions

const μm_actual = 1.76
const σm_actual = 0.0635
const μf_actual = 1.62
const σf_actual = 0.0559

ActualParameters = [μm_actual,σm_actual,μf_actual,σf_actual]

function generatedata(;numofdatum = 100)
    rng = StableRNG(12345)
    data = zeros(numofdatum)
    # 55% male with 1.76 m and std dev of 0.0635
    for k = 1:Int64(55//100*numofdatum)
        data[k] = rand(rng,Normal(μm_actual,σm_actual))
    end
    # 45% female with 1.62 m and std dev of 0.0559
    for k = (1+Int64(55//100*numofdatum)):numofdatum
        data[k] = rand(rng,Normal(μf_actual,σf_actual))
    end
    return round.(data,digits=2)
end

# Generate a numeric score so that we can compare how MCMC is performing
function errorscore(actual_parameter, estimated_parameter, scalingfactor, weighting)
    local temp
    temp = actual_parameter .- estimated_parameter
    temp = temp ./ scalingfactor
    temp = temp .^ 2
    temp = temp .* weighting
    return round(sum(temp),sigdigits=3)
end

function MCMC_example()

    data = generatedata()

    μm = Normal(1.6, 0.5)
    σm = Uniform(0.0,1.0)
    μf = Normal(1.6, 0.5)
    σf = Uniform(0.0,1.0)

    # Put all our parameters into a list
    ListOfParams = Distribution[μm,σm,μf,σf]

    # The loglikelihood function which returns the log(probability),
    # given a fixed mu and a fixed sigma, of the sequence of datapoints
    function loglikelihood(data::Vector,μm,σm,μf,σf)
        numofdatapoints = length(data)
        result=0.0
        @inbounds for k = 1:numofdatapoints
            mmm = 0.55*pdf(Normal(μm,σm),data[k])
            fff = 0.45*pdf(Normal(μf,σf),data[k])
            result += log(   mmm + fff   )
        end
        return result
    end

    p = MetropolisHastings_MCMC(
            data,
            ListOfParams,
            loglikelihood,
            samples=1_000_000,
            burnedinsamples=1000
    )

    # Next we need to display the histogram of the column p[:,1]
    # this is the shape of the posterior produced by MCMC
    vec_μm=vec(p[:,1])  # Convert a Matrix column into a 1-D Vector
    vec_σm=vec(p[:,2])  # Convert a Matrix column into a 1-D Vector
    vec_μf=vec(p[:,3])  # Convert a Matrix column into a 1-D Vector
    vec_σf=vec(p[:,4])  # Convert a Matrix column into a 1-D Vector
    #describe_paramvec("μ","mu",vec_μ)
    #describe_paramvec("σ","sigma",vec_σ)
    t = describe_multiparamvec([("μm","mu male",vec_μm),
                            ("σm","sigma male",vec_σm),
                            ("μf","mu female",vec_μf),
                            ("σf","sigma female",vec_σf)
                            ],savegraph=true)
    histogram2d(vec_μm,vec_σm,nbins=100,xlabel="mu male",ylabel="sigma male") |> display
    savefig("histogram2D_mu_male_sigma_male.png")
    histogram2d(vec_μf,vec_σf,nbins=100,xlabel="mu female",ylabel="sigma female") |> display
    savefig("histogram2D_mu_female_sigma_female.png")
    histogram2d(vec_μm,vec_μf,nbins=100,xlabel="mu male",ylabel="mu female") |> display
    savefig("histogram2D_mu_male_mu_female.png")
    histogram2d(vec_σm,vec_σf,nbins=100,xlabel="sigma male",ylabel="sigma female") |> display
    savefig("histogram2D_sigma_male_sigma_female.png")

    estimated = round.(t[:,2],sigdigits=3)

    println("error score is ",
        errorscore(
        ActualParameters,
        estimated,
        [0.01,0.01,0.01,0.01],
        [1/4,1/4,1/4,1/4]
        )
    )
end

# We call the main function to kick off the example
MCMC_example()

Prior

Now we need to explain our prior.

We knows most humans are around 1.6 metres tall and their standard deviation is a lot less than 0.5 metres. So we set the prior for mean for height to be:

Normal(1.6,0.5) for both males and females.

We set the prior for standard deviation for height to be

Uniform(0.0,1.0)

Because it has to be greater than zero but less than 1 metre

loglikelihood function

The pdf (probability density function) of the distribution of height for the students is the weighted combination of the pdf distribution of both the genders of the students.

0.55 * pdf(male height distribution) + 0.45 * pdf(female height distribution)

function loglikelihood(data::Vector,μm,σm,μf,σf)
    numofdatapoints = length(data)
    result=0.0
    @inbounds for k = 1:numofdatapoints
        mmm = 0.55*pdf(Normal(μm,σm),data[k])
        fff = 0.45*pdf(Normal(μf,σf),data[k])
        result += log(   mmm + fff   )
    end
    return result
end

results

So here are the results

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ julia --project=. bothgender01.jl

Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
[ Info: Precompiling SimpleMCMC [top-level]
┌───────────┬────────┬────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev │ 5.5%   │ 94.5%  │
├───────────┼────────┼────────┼────────┼────────┤
│ μm        │ 1.67   │ 0.0649 │ 1.6    │ 1.77   │
│ σm        │ 0.071  │ 0.0136 │ 0.0525 │ 0.0957 │
│ μf        │ 1.71   │ 0.0781 │ 1.59   │ 1.79   │
│ σf        │ 0.0614 │ 0.016  │ 0.0419 │ 0.0921 │
└───────────┴────────┴────────┴────────┴────────┘
error score is 45.5

Analysis of First attempt

The first thing we noticed about results is that MCMC thinks that the mean of the male height has two peaks. They are 1.77 m and 1.62 m.

It also thinks that the mean of the female height has two peaks. They are also 1.77 m and 1.62 m.

This woukld be clearer when we look at the 2D histogram of “mu male” vs “mu female”. It shows two distinct peaks.

one peak is (mu male=1.62,mu female=1.77)
the other peak is (mu male=1.77,mu female=1.66)

After thinking about it, it became clear that these are the two possibilities that will generate the distribution of the overall student heights. The MCMC does NOT know which of the two possibilities is responsible for the datapoints.

Where as we humans know (through our common sense) that the average height of males is higher than the average height of females.

So what we need to do is to inform MCMC of the fact that mean of male height is ALWAYS higher than the mean of the female height.

Second attempt

Let's alter the source code to let the MCMC know that mean of males are higher than mean of females.

We do this by changing the loglikelihood function to reduce each component of the result (per datapoint) by -1.0 if μf > μm.

we use (the reverse logic)

(μm > μf ? 0.0 : -1.0)

# Filename is "bothgender02.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"
include("/Users/ssiew/juliascript/MCMC_dir/SimpleMCMC.jl")

using Distributions, StatsPlots, Statistics
using .SimpleMCMC
using StableRNGs,Distributions

const μm_actual = 1.76
const σm_actual = 0.0635
const μf_actual = 1.62
const σf_actual = 0.0559

ActualParameters = [μm_actual,σm_actual,μf_actual,σf_actual]

function generatedata(;numofdatum = 100)
    rng = StableRNG(12345)
    data = zeros(numofdatum)
    # 55% male with 1.76 m and std dev of 0.0635
    for k = 1:Int64(55//100*numofdatum)
        data[k] = rand(rng,Normal(μm_actual,σm_actual))
    end
    # 45% female with 1.62 m and std dev of 0.0559
    for k = (1+Int64(55//100*numofdatum)):numofdatum
        data[k] = rand(rng,Normal(μf_actual,σf_actual))
    end
    return round.(data,digits=2)
end

# Generate a numeric score so that we can compare how MCMC is performing
function errorscore(actual_parameter, estimated_parameter, scalingfactor, weighting)
    local temp
    temp = actual_parameter .- estimated_parameter
    temp = temp ./ scalingfactor
    temp = temp .^ 2
    temp = temp .* weighting
    return round(sum(temp),sigdigits=3)
end

function MCMC_example()

    data = generatedata()

    μm = Normal(1.6, 0.5)
    σm = Uniform(0.0,1.0)
    μf = Normal(1.6, 0.5)
    σf = Uniform(0.0,1.0)

    # Put all our parameters into a list
    ListOfParams = Distribution[μm,σm,μf,σf]

    # The loglikelihood function which returns the log(probability),
    # given a fixed mu and a fixed sigma, of the sequence of datapoints
    function loglikelihood(data::Vector,μm,σm,μf,σf)
        numofdatapoints = length(data)
        result=0.0
        @inbounds for k = 1:numofdatapoints
            mmm = 0.55*pdf(Normal(μm,σm),data[k])
            fff = 0.45*pdf(Normal(μf,σf),data[k])
            result += (  log(   mmm + fff   )     +
                         (μm > μf ? 0.0 : -1.0)
                      )
        end
        return result
    end

    p = MetropolisHastings_MCMC(
            data,
            ListOfParams,
            loglikelihood,
            samples=1_000_000,
            burnedinsamples=1000
    )

    # Next we need to display the histogram of the column p[:,1]
    # this is the shape of the posterior produced by MCMC
    vec_μm=vec(p[:,1])  # Convert a Matrix column into a 1-D Vector
    vec_σm=vec(p[:,2])  # Convert a Matrix column into a 1-D Vector
    vec_μf=vec(p[:,3])  # Convert a Matrix column into a 1-D Vector
    vec_σf=vec(p[:,4])  # Convert a Matrix column into a 1-D Vector
    #describe_paramvec("μ","mu",vec_μ)
    #describe_paramvec("σ","sigma",vec_σ)
    t = describe_multiparamvec([("μm","mu male",vec_μm),
                            ("σm","sigma male",vec_σm),
                            ("μf","mu female",vec_μf),
                            ("σf","sigma female",vec_σf)
                            ],savegraph=true)
    histogram2d(vec_μm,vec_σm,nbins=100,xlabel="mu male",ylabel="sigma male") |> display
    savefig("histogram2D_mu_male_sigma_male.png")
    histogram2d(vec_μf,vec_σf,nbins=100,xlabel="mu female",ylabel="sigma female") |> display
    savefig("histogram2D_mu_female_sigma_female.png")
    histogram2d(vec_μm,vec_μf,nbins=100,xlabel="mu male",ylabel="mu female") |> display
    savefig("histogram2D_mu_male_mu_female.png")
    histogram2d(vec_σm,vec_σf,nbins=100,xlabel="sigma male",ylabel="sigma female") |> display
    savefig("histogram2D_sigma_male_sigma_female.png")

    estimated = round.(t[:,2],sigdigits=3)

    println("error score is ",
        errorscore(
        ActualParameters,
        estimated,
        [0.01,0.01,0.01,0.01],
        [1/4,1/4,1/4,1/4]
        )
    )
end

# We call the main function to kick off the example
MCMC_example()

results Second Attempt

So here are the results

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ julia --project=. bothgender02.jl

Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
┌───────────┬────────┬────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev │ 5.5%   │ 94.5%  │
├───────────┼────────┼────────┼────────┼────────┤
│ μm        │ 1.75   │ 0.0195 │ 1.71   │ 1.78   │
│ σm        │ 0.0664 │ 0.0142 │ 0.0487 │ 0.0938 │
│ μf        │ 1.61   │ 0.0209 │ 1.59   │ 1.66   │
│ σf        │ 0.0702 │ 0.0151 │ 0.0524 │ 0.1    │
└───────────┴────────┴────────┴────────┴────────┘
error score is 1.03

Analysis of Second attempt

Now the results are so much better than the first attempt. For one, we got a single peak for the mean height for males and a single peak for the mean height for females.

And the 2D histogram of “mu male” vs “mu female” also shows a single peak. This is very good news indeed!

Fit booster using Turbo Cycle

Is there any other method where the error score can be further reduced? Well, we could use the turbo mode.

How Turbo Cycle works

In a car, the basic concept of a Turbo Charger is to use the ultra high pressure of the exhaust gas to power an air compressor in order to compress the intake air to an higher air pressure before it enters an engine.

But why do turbo chargers work? The secret is that the amount of energy liberated by the internal combustion engine is dependant on two things.

amount of fuel molecules
amount of oxygen molecules

While we have no problems with putting in as much fuel as we want. The only way to increase the amount of oxygen molecules is to increase the air pressure of the intake air using an air compressor. Otherwise we are LIMITED by

18% oxygen gas in the atmosphere (the rest is mostly Nitrogen gas)
1.0 standard air pressure in the atmosphere

Back to explaination, we use use the POWER OF THE EXHAUST to increase the QUALITY OF THE AIR INTAKE. or you can say "We use the HIGH PRESSURE of the exhaust gas to further increase the air pressure of the air intake."

With MCMC, we have the exactly same issue.

The PRIOR of the input to MCMC is basically very LOW QUALITY.
The output of the MCMC is HIGHER QUALITY POSTERIOR

The idea is to use the "HIGHER QUALITY POSTERIOR" as the PRIOR for the next cycle of the MCMC input.

Limitation of Turbo Cycle

The limiting factor of course is the Fuel or Data from the experiment.

No matter how much oxygen molecules you have in the HIGH PRESSURE AIR INTAKE, if there is only so much fuel then you can only get so much power from the engine.

In MCMC, the data from the experiment is the same limiting factor on your POSTERIOR.

This is why you should only use a MaxTurboCycle of 2 (or a low number).

Implementing Turbo Cycle

This is the code to implement turbo cycle.

# Filename is "bothgender02_turbo.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"
include("/Users/ssiew/juliascript/MCMC_dir/SimpleMCMC.jl")

using Distributions, StatsPlots, Statistics
using .SimpleMCMC
using StableRNGs,Distributions

const μm_actual = 1.76
const σm_actual = 0.0635
const μf_actual = 1.62
const σf_actual = 0.0559

ActualParameters = [μm_actual,σm_actual,μf_actual,σf_actual]

function generatedata(;numofdatum = 100)
    rng = StableRNG(12345)
    data = zeros(numofdatum)
    # 55% male with 1.76 m and std dev of 0.0635
    for k = 1:Int64(55//100*numofdatum)
        data[k] = rand(rng,Normal(μm_actual,σm_actual))
    end
    # 45% female with 1.62 m and std dev of 0.0559
    for k = (1+Int64(55//100*numofdatum)):numofdatum
        data[k] = rand(rng,Normal(μf_actual,σf_actual))
    end
    return round.(data,digits=2)
end

# Generate a numeric score so that we can compare how MCMC is performing
function errorscore(actual_parameter, estimated_parameter, scalingfactor, weighting)
    local temp
    temp = actual_parameter .- estimated_parameter
    temp = temp ./ scalingfactor
    temp = temp .^ 2
    temp = temp .* weighting
    return round(sum(temp),sigdigits=3)
end

function MCMC_example()

    data = generatedata()

    μm = Normal(1.6, 0.5)
    σm = Uniform(0.0,1.0)
    μf = Normal(1.6, 0.5)
    σf = Uniform(0.0,1.0)

    # Put all our parameters into a list
    ListOfParams = Distribution[μm,σm,μf,σf]

    # The loglikelihood function which returns the log(probability),
    # given a fixed mu and a fixed sigma, of the sequence of datapoints
    function loglikelihood(data::Vector,μm,σm,μf,σf)
        numofdatapoints = length(data)
        result=0.0
        @inbounds for k = 1:numofdatapoints
            mmm = 0.55*pdf(Normal(μm,σm),data[k])
            fff = 0.45*pdf(Normal(μf,σf),data[k])
            result += (  log(   mmm + fff   )     +
                         (μm > μf ? 0.0 : -1.0)
                      )
        end
        return result
    end

    local vec_μm,vec_σm,vec_μf,vec_σf
    local MaxTurboCycle = 2
    for TurboCycle = 0:MaxTurboCycle
        println("==================")
        println("TurboCycle is ",TurboCycle)
        p = MetropolisHastings_MCMC(
                data,
                ListOfParams,
                loglikelihood,
                samples=1_000_000,
                burnedinsamples=1000
        )

        # Next we need to display the histogram of the column p[:,1]
        # this is the shape of the posterior produced by MCMC
        vec_μm=vec(p[:,1])  # Convert a Matrix column into a 1-D Vector
        vec_σm=vec(p[:,2])  # Convert a Matrix column into a 1-D Vector
        vec_μf=vec(p[:,3])  # Convert a Matrix column into a 1-D Vector
        vec_σf=vec(p[:,4])  # Convert a Matrix column into a 1-D Vector
        #describe_paramvec("μ","mu",vec_μ)
        #describe_paramvec("σ","sigma",vec_σ)
        t = describe_multiparamvec([("μm","mu male",vec_μm),
                                ("σm","sigma male",vec_σm),
                                ("μf","mu female",vec_μf),
                                ("σf","sigma female",vec_σf)
                                ],savegraph=true)

        estimated = round.(t[:,2],sigdigits=3)

        println("error score is ",
            errorscore(
            ActualParameters,
            estimated,
            [0.01,0.01,0.01,0.01],
            [1/4,1/4,1/4,1/4]
            )
        )

        ListOfParams = Distribution[ Normal(round(t[k,2],sigdigits=3),round(t[k,3],sigdigits=3))  for k=1:4 ]
    end

    histogram2d(vec_μm,vec_σm,nbins=100,xlabel="mu male",ylabel="sigma male") |> display
    savefig("histogram2D_mu_male_sigma_male.png")
    histogram2d(vec_μf,vec_σf,nbins=100,xlabel="mu female",ylabel="sigma female") |> display
    savefig("histogram2D_mu_female_sigma_female.png")
    histogram2d(vec_μm,vec_μf,nbins=100,xlabel="mu male",ylabel="mu female") |> display
    savefig("histogram2D_mu_male_mu_female.png")
    histogram2d(vec_σm,vec_σf,nbins=100,xlabel="sigma male",ylabel="sigma female") |> display
    savefig("histogram2D_sigma_male_sigma_female.png")

end

# We call the main function to kick off the example
MCMC_example()

Turbo Cycle result

So here are the results

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ julia --project=. bothgender02_turbo.jl

Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
==================
TurboCycle is 0
┌───────────┬────────┬────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev │ 5.5%   │ 94.5%  │
├───────────┼────────┼────────┼────────┼────────┤
│ μm        │ 1.75   │ 0.0196 │ 1.71   │ 1.78   │
│ σm        │ 0.0667 │ 0.0193 │ 0.0486 │ 0.094  │
│ μf        │ 1.61   │ 0.0211 │ 1.59   │ 1.66   │
│ σf        │ 0.0701 │ 0.0151 │ 0.0523 │ 0.0997 │
└───────────┴────────┴────────┴────────┴────────┘
error score is 1.03
==================
TurboCycle is 1
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.76   │ 0.0107  │ 1.74   │ 1.77   │
│ σm        │ 0.0638 │ 0.00883 │ 0.0507 │ 0.0788 │
│ μf        │ 1.61   │ 0.0114  │ 1.59   │ 1.63   │
│ σf        │ 0.0669 │ 0.00871 │ 0.054  │ 0.0819 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.553
==================
TurboCycle is 2
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.76   │ 0.00807 │ 1.75   │ 1.77   │
│ σm        │ 0.0622 │ 0.00628 │ 0.0526 │ 0.0726 │
│ μf        │ 1.61   │ 0.00827 │ 1.59   │ 1.62   │
│ σf        │ 0.0657 │ 0.00639 │ 0.0559 │ 0.0763 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.494

This is the limit on what we can know about the parameters without further specialized knowledge about the genders.

The best way to increase the accuracy of MCMC

The best way to increase the accuracy of MCMC is to get more datapoints.

For example: If we got 1000 datapoints instead of 100 datapoints, we can get much better results as seen below

$ julia --project=. bothgender02_turbo_d1000.jl
==================
turbonum is 0
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.76   │ 0.0497  │ 1.75   │ 1.77   │
│ σm        │ 0.0643 │ 0.00578 │ 0.0598 │ 0.0689 │
│ μf        │ 1.62   │ 0.00548 │ 1.62   │ 1.63   │
│ σf        │ 0.0533 │ 0.00327 │ 0.0492 │ 0.0573 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.0185
==================
turbonum is 1
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.76   │ 0.00372 │ 1.75   │ 1.77   │
│ σm        │ 0.0641 │ 0.00252 │ 0.0602 │ 0.0682 │
│ μf        │ 1.62   │ 0.0028  │ 1.62   │ 1.63   │
│ σf        │ 0.0528 │ 0.00189 │ 0.0498 │ 0.0559 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.0249
==================
turbonum is 2
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.76   │ 0.00257 │ 1.76   │ 1.76   │
│ σm        │ 0.064  │ 0.00181 │ 0.0611 │ 0.0669 │
│ μf        │ 1.62   │ 0.00205 │ 1.62   │ 1.62   │
│ σf        │ 0.0526 │ 0.00148 │ 0.0502 │ 0.0549 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.0278

Further increasing in fitness using specialized knowledge about genders

In a nutshell, the more we know about the characteristic of the gender of the students, the more we can further increase th fitness of the estimated parameters.

Let us suppose that we know the following about the genders of the students BEFORE we start the experiment.

The sigma_of_male_student is greater than sigma_of_female_student. In other words, the variation of the height of the male students is greater than the variation of the height of the female students.
mu_of_male_student - mu_of_female_student = 0.14 In other words the average height of the male students is 14 cm higher than the average height of the female students

We can represent it in the loglikelihood function as

    function loglikelihood(data::Vector,μm,σm,μf,σf)
        numofdatapoints = length(data)
        result=0.0
        @inbounds for k = 1:numofdatapoints
            mmm = 0.55*pdf(Normal(μm,σm),data[k])
            fff = 0.45*pdf(Normal(μf,σf),data[k])
            result += (  log(   mmm + fff   )     +
                         (μm > μf ? 0.0 : -1.0)   +
                         (σm > σf ? 0.0 : -1.0)   +
                         (-1.0 * min( abs((μm - (μf+0.14))/0.05),1.0 )  )
                      )
        end
        return result
    end

So our sourcecode is

# Filename is "bothgender04_turbo.jl"
# Directory is "/Users/ssiew/juliascript/MCMC_dir"
include("/Users/ssiew/juliascript/MCMC_dir/SimpleMCMC.jl")

using Distributions, StatsPlots, Statistics
using .SimpleMCMC
using StableRNGs,Distributions

const μm_actual = 1.76
const σm_actual = 0.0635
const μf_actual = 1.62
const σf_actual = 0.0559

ActualParameters = [μm_actual,σm_actual,μf_actual,σf_actual]

function generatedata(;numofdatum = 100)
    rng = StableRNG(12345)
    data = zeros(numofdatum)
    # 55% male with 1.76 m and std dev of 0.0635
    for k = 1:Int64(55//100*numofdatum)
        data[k] = rand(rng,Normal(μm_actual,σm_actual))
    end
    # 45% female with 1.62 m and std dev of 0.0559
    for k = (1+Int64(55//100*numofdatum)):numofdatum
        data[k] = rand(rng,Normal(μf_actual,σf_actual))
    end
    return round.(data,digits=2)
end

# Generate a numeric score so that we can compare how MCMC is performing
function errorscore(actual_parameter, estimated_parameter, scalingfactor, weighting)
    local temp
    temp = actual_parameter .- estimated_parameter
    temp = temp ./ scalingfactor
    temp = temp .^ 2
    temp = temp .* weighting
    return round(sum(temp),sigdigits=3)
end

function MCMC_example()

    data = generatedata()

    μm = Normal(1.6, 0.5)
    σm = Uniform(0.0,1.0)
    μf = Normal(1.6, 0.5)
    σf = Uniform(0.0,1.0)

    # Put all our parameters into a list
    ListOfParams = Distribution[μm,σm,μf,σf]

    # The loglikelihood function which returns the log(probability),
    # given a fixed mu and a fixed sigma, of the sequence of datapoints
    function loglikelihood(data::Vector,μm,σm,μf,σf)
        numofdatapoints = length(data)
        result=0.0
        @inbounds for k = 1:numofdatapoints
            mmm = 0.55*pdf(Normal(μm,σm),data[k])
            fff = 0.45*pdf(Normal(μf,σf),data[k])
            result += (  log(   mmm + fff   )     +
                         (μm > μf ? 0.0 : -1.0)   +
                         (σm > σf ? 0.0 : -1.0)   +
                         (-1.0 * min( abs((μm - (μf+0.14))/0.01),1.0 )  )
                      )
        end
        return result
    end

    local vec_μm,vec_σm,vec_μf,vec_σf
    local MaxTurboCycle = 2
    for TurboCycle = 0:MaxTurboCycle
        println("==================")
        println("TurboCycle is ",TurboCycle)
        p = MetropolisHastings_MCMC(
                data,
                ListOfParams,
                loglikelihood,
                samples=1_000_000,
                burnedinsamples=1000
        )

        # Next we need to display the histogram of the column p[:,1]
        # this is the shape of the posterior produced by MCMC
        vec_μm=vec(p[:,1])  # Convert a Matrix column into a 1-D Vector
        vec_σm=vec(p[:,2])  # Convert a Matrix column into a 1-D Vector
        vec_μf=vec(p[:,3])  # Convert a Matrix column into a 1-D Vector
        vec_σf=vec(p[:,4])  # Convert a Matrix column into a 1-D Vector
        #describe_paramvec("μ","mu",vec_μ)
        #describe_paramvec("σ","sigma",vec_σ)
        t = describe_multiparamvec([("μm","mu male",vec_μm),
                                ("σm","sigma male",vec_σm),
                                ("μf","mu female",vec_μf),
                                ("σf","sigma female",vec_σf)
                                ],savegraph=true)

        estimated = round.(t[:,2],sigdigits=3)

        println("error score is ",
            errorscore(
            ActualParameters,
            estimated,
            [0.01,0.01,0.01,0.01],
            [1/4,1/4,1/4,1/4]
            )
        )

        ListOfParams = Distribution[ Normal(round(t[k,2],sigdigits=3),round(t[k,3],sigdigits=3))  for k=1:4 ]
    end

    histogram2d(vec_μm,vec_σm,nbins=100,xlabel="mu male",ylabel="sigma male") |> display
    savefig("histogram2D_mu_male_sigma_male.png")
    histogram2d(vec_μf,vec_σf,nbins=100,xlabel="mu female",ylabel="sigma female") |> display
    savefig("histogram2D_mu_female_sigma_female.png")
    histogram2d(vec_μm,vec_μf,nbins=100,xlabel="mu male",ylabel="mu female") |> display
    savefig("histogram2D_mu_male_mu_female.png")
    histogram2d(vec_σm,vec_σf,nbins=100,xlabel="sigma male",ylabel="sigma female") |> display
    savefig("histogram2D_sigma_male_sigma_female.png")

end

# We call the main function to kick off the example
MCMC_example()

So here are the results

$ echo "This is the prompt on Terminal on MacOS"
$ cd /Users/ssiew/juliascript/MCMC_dir
$ julia --project=. bothgender04_turbo.jl

Starting Julia...
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 1.8.5 (2023-01-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |
==================
TurboCycle is 0
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.75   │ 0.0094  │ 1.73   │ 1.76   │
│ σm        │ 0.0718 │ 0.00942 │ 0.0592 │ 0.0889 │
│ μf        │ 1.61   │ 0.0094  │ 1.59   │ 1.62   │
│ σf        │ 0.0624 │ 0.00804 │ 0.0514 │ 0.0771 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.778
==================
TurboCycle is 1
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.75   │ 0.00542 │ 1.74   │ 1.76   │
│ σm        │ 0.0703 │ 0.00581 │ 0.0615 │ 0.0801 │
│ μf        │ 1.61   │ 0.00542 │ 1.6    │ 1.62   │
│ σf        │ 0.0621 │ 0.00496 │ 0.0544 │ 0.0703 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.712
==================
TurboCycle is 2
┌───────────┬────────┬─────────┬────────┬────────┐
│ Parameter │ Mean   │ Stddev  │ 5.5%   │ 94.5%  │
├───────────┼────────┼─────────┼────────┼────────┤
│ μm        │ 1.75   │ 0.00334 │ 1.74   │ 1.76   │
│ σm        │ 0.0698 │ 0.00432 │ 0.0631 │ 0.0769 │
│ μf        │ 1.61   │ 0.00334 │ 1.6    │ 1.62   │
│ σf        │ 0.0623 │ 0.00388 │ 0.0561 │ 0.0686 │
└───────────┴────────┴─────────┴────────┴────────┘
error score is 0.702

Using Riemann Zeta Function and DoubleFloats to calculate Pi

Steven Siew — Sat, 11 Mar 2023 03:16:05 +0000

Riemann Zeta Function

It has the following particular values

DoubleFloats

DoubleFloats uses two Float64 instances to represent a single floating point number with 85 significant bits.

In this article we use it to calculate Pi because it can calculate pi to about 30 decimal digits where as normal Float64 only have 15 decimal digits.

Program

using DoubleFloats

function calc_zeta_02(n)
    local accumulator = Double64(0)
    for k = n:-1:1
        accumulator += Double64(1)/(Double64(k)^2)
    end
    result = Double64(6) * accumulator
    result = sqrt(result)
    return result
end

function calc_zeta_04(n)
    local accumulator = Double64(0)
    for k = n:-1:1
        accumulator += Double64(1)/(Double64(k)^4)
    end
    result = Double64(90) * accumulator
    result = sqrt(sqrt(result))
    return result
end

function calc_zeta_08(n)
    local accumulator = Double64(0)
    for k = n:-1:1
        accumulator += Double64(1)/(Double64(k)^8)
    end
    result = Double64(9450) * accumulator
    result = sqrt(sqrt(sqrt(result)))
    return result
end

function calc_zeta_16(n)
    local accumulator = Double64(0)
    for k = n:-1:1
        accumulator += Double64(1)/(Double64(k)^16)
    end
    result = Double64("325641566250") / Double64(3617) * accumulator
    result = sqrt(sqrt(sqrt(sqrt(result))))
    return result
end

println("Start of program")
println("calc_zeta_02(40) = ",string(calc_zeta_02(40))  )
println("calc_zeta_04(40) = ",string(calc_zeta_04(40))  )
println("calc_zeta_08(40) = ",string(calc_zeta_08(40))  )
println("calc_zeta_16(40) = ",string(calc_zeta_16(40))  )
println("              pi = ",string(Double64(pi))     )

It gives the following output

Start of program
calc_zeta_02(40) = 3.11792619829937803643280481368424029
calc_zeta_04(40) = 3.14158901347571579803199355470915640
calc_zeta_08(40) = 3.14159265358948106767954153673911126
calc_zeta_16(40) = 3.14159265358979323846264337322256472
              pi = 3.14159265358979323846264338327950588

When you are adding up lots of floating point numbers, always add up the smallest number first before adding the bigger numbers. That is why we start at the highest value of k first and systematically reduce k until it reaches the value one.

Why stop at Riemann Zeta 16 when you can go even higher

When something is worth doing, it's worth overdoing.

We will just do one more to show you how powerful Riemann Zeta 32 is

function calc_zeta_32(n)
    local accumulator = Double64(0)
    for k = n:-1:1
        accumulator += Double64(1)/(Double64(k)^32)
    end
    result = Double64("62490220571022341207266406250") / Double64("7709321041217") * accumulator
    result = sqrt(sqrt(sqrt(sqrt(sqrt(result)))))
    return result
end

println("")
println("calc_zeta_16(40) = ",string(calc_zeta_16(40))  )
println("calc_zeta_32(40) = ",string(calc_zeta_32(40))  )
println("              pi = ",string(Double64(pi))     )

with the output of


calc_zeta_16(40) = 3.14159265358979323846264337322256472
calc_zeta_32(40) = 3.14159265358979323846264338327950588
              pi = 3.14159265358979323846264338327950588

Amazingly it only took 9 iterations to get to the same value of the constant pi in doublefloats format

julia> calc_zeta_32(9) |> string
"3.14159265358979323846264338327950588"

Are there any alternative to writing sqrt(sqrt(sqrt(sqrt(sqrt(result)))))?

Of course there are! First we find the result of sqrt(sqrt(sqrt(sqrt(sqrt(123_456_789.0)))))

julia> sqrt(sqrt(sqrt(sqrt(sqrt(123_456_789.0)))))
1.7900280769789536

Please notice that there are 5 sqrt in the statement above.

The solution is to create a function f such that

    #=
       func_compose(func,n)

    Compose a function multiple times
    =#
    function func_compose(func,n)
        newfunc = func
        local loop = n
        while loop > 1
            #  help?> ∘
            #  "∘" can be typed by \circ<tab>
            newfunc = func ∘ newfunc
            loop -= 1
        end
        return newfunc
    end

f = func_compose(sqrt,5)
println("f(123_456_789.0) = ",f(123_456_789.0))

with the output

f(123_456_789.0) = 1.7900280769789536

Note that '∘' is the Function composition operator

You don't even need to save the result of func_compose(sqrt,5) to the variable f

You can use it directly like this

julia> func_compose(sqrt,5)(123_456_789.0)
1.7900280769789536

You can also turn func_compose() into a binary operator like this

help?> ⋏
"⋏" can be typed by \curlywedge<tab>

julia> ⋏(a,b)=func_compose(a,b)
⋏ (generic function with 1 method)

julia> (sqrt ⋏ 5)(123_456_789.0)
1.7900280769789536

julia> g = sqrt ⋏ 5
sqrt ∘ sqrt ∘ sqrt ∘ sqrt ∘ sqrt

julia> g(123_456_789.0)
1.7900280769789536

How to customize existing base Julia functions

Steven Siew — Sun, 05 Mar 2023 00:23:36 +0000

How to customize existing base Julia function

Have you ever wish you could create a customize version of an existing Base function.

Here is a simple trick to do it.

julia> myprint(a,b...;kwargs...) = begin print("> "); print(a,b...;kwargs...); end
myprint (generic function with 1 method)

julia> myprint("hello world\n")
> hello world

julia> myprintln(a,b...;kwargs...) = begin print("> "); println(a,b...;kwargs...); end
myprintln (generic function with 1 method)

julia> myprintln("hello ",56.7," world")
> hello 56.7 world

Now you have the functions myprint() and myprintln() which behaves how you wanted it with all the functionality of print() and println()

You can also create it this way

function myprint(a,b...;kwargs...)
    print("> ")
    print(a,b...;kwargs...)
end

Using Measurements.jl for Triangulation in Julia

Steven Siew — Fri, 24 Feb 2023 07:56:33 +0000

Finding the location of Cerberus

HMVS Cerberus is a warship that was "sunk" in blackrock. You can see the ship from shore. The objective is to find the location of HMS Cerberus without a boat.

You can measure the bearing to the ship from two locations

Shore : bearing to ship is 263.3 degrees
Pier : bearing to ship is 305.8 degrees

Using Google maps, you can get the location of the two shore locations. To get the location, right click on the spot you are interested in.

Shore : location is -37.96705678343727 (latitude), 145.01158473574677 (longtitude)
Pier : location is -37.968305948144405 (latitude), 145.00959897290315 (longtitude)

To get the location into an XY coordinate system, we use this website and set the earth geographic datum as WGS84.

https://awsm-tools.com/lat-long-to-utm?form%5Bellipsoid%5D=WGS+84

This gives us the WGS84 GPS coordinates system where the numbers are in units of meters/metres. The three numbers are

UTM ZONE (eg. 55H which is zone number 55 with Latitude Band H)
Easting (number of meters east from a virtual base location)
Northing (number of meters north from a virtual base location)

You can read more about the UTM coordinate system below
https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system

Our details are

Shore : location is 325335.77856296976 (E) 5795975.086601857 (N) 55H (UTM Zone)
Pier : location is 325164.2938595414 (E) 5795832.741503723 (N) 55H (UTM Zone)

Now that we got all the information that we needed, we can calculate the location of the Cerebus/target

Julia source code to calculate location

First we give you the entire source code in Julia so that you can copy and paste it into a textfile.

This source code requires the following packages

Printf
Measurements

using Printf

#=  Crosspoint_print(x)
    Print coordinates to two decimal places
=#
function Crosspoint_print(x)
    if length(x) == 1
        print(rstrip(@sprintf("%-10.2f", x)))
        return
    end
    local flag = false
    print("[")
    for ele in x
        print(flag==true ? "," : "")
        print(rstrip(@sprintf("%-10.2f", ele)))
        flag = true
    end
    print("]")
end

Crosspoint_println(x) = begin Crosspoint_print(x); println(); end

#=  intersect2lines(A,a,B,b)

    Return the intersection of two lines from
    two known points A,B using sin_degree() and
    cos_degree() only. a and b are the bearings
    from the unknown position to the two known points A,B
=#
function intersect2lines(A,a,B,b)
    den = sind(a - b)    # sind() function is Sine in degrees, cosd() is similar
    tempX = B[1]*cosd(b)*sind(a) - (A[1]*cosd(a) + (B[2] - A[2])*sind(a))*sind(b)
    tempY = A[2]*cosd(b)*sind(a) - cosd(a)*((A[1] - B[1])*cosd(b) + B[2]*sind(b))
    newX = tempX/den
    newY = tempY/den
    return [newX,newY]
end

#=  triangulation(A,a,B,b)

    Return the position of an unknown location
    by triagulation from two known points A,B and
    the bearing from those two known points A,B to
    the unknown location
=#
triangulation(A,a,B,b) = intersect2lines(A,a+180.0,B,b+180.0)


Pos_Shore = (325335.78,5795975.09)    # Easting 325335.78 Northing 5795975.09
Pos_Shore_bearing_to_target = 263.3   # 263.3 degrees

Pos_Pier  = (325164.29,5795832.74)    # Easting 325164.29 Northing 5795832.74
Pos_Pier_bearing_to_target  = 305.8   # 305.8 degrees


println("program start")

println("The coordinates for :")
print("Shore is ")
Crosspoint_println(Pos_Shore)
print("Pier is ")
Crosspoint_println(Pos_Pier)
println("")
println("Bearing from shore to cerberus/target is ",round(Pos_Shore_bearing_to_target,digits=1))
println("Bearing from pier to cerberus/target is ",round(Pos_Pier_bearing_to_target,digits=1))
println("")

target_pos_estimate = triangulation(Pos_Shore,
                                    Pos_Shore_bearing_to_target,
                                    Pos_Pier,
                                    Pos_Pier_bearing_to_target)

print("The estimated position of target is ")
Crosspoint_println(target_pos_estimate)
println("In coordinate format of [Easting value,Northing value]")

println("")
println("Now we do the same thing again using measurement")
println("")

#
# Now use Measurement to measure the uncertainties
#    help?>  ±
#    "±" can be typed by \pm<tab>

using Measurements

Shore = (325335.78 ± 5.0,5795975.09 ± 5.0)
Pier  = (325164.29 ± 5.0,5795832.74 ± 5.0)
Shore_bearing_to_target = 263.3 ± 0.5
Pier_bearing_to_target  = 305.8 ± 0.5

target_pos_estimate2 = triangulation(Shore,
                                     Shore_bearing_to_target,
                                     Pier,
                                     Pier_bearing_to_target)

value_two_decimal_places(m)       = rstrip(  @sprintf("%-16.2f",Measurements.value(m))  )
uncertainty_two_decimal_places(m) = rstrip(  @sprintf("%-16.2f",Measurements.uncertainty(m))  )

print("The estimated position of target using Measurements package is ")
print("[")
print(value_two_decimal_places(target_pos_estimate2[1]))
print(" ± ")
print(uncertainty_two_decimal_places(target_pos_estimate2[1]))
print(",")
print(value_two_decimal_places(target_pos_estimate2[2]))
print(" ± ")
print(uncertainty_two_decimal_places(target_pos_estimate2[2]))
println("]")
println("")

dict_uncertainty_in_easting  = Measurements.uncertainty_components(target_pos_estimate2[1])
dict_uncertainty_in_northing = Measurements.uncertainty_components(target_pos_estimate2[2])

println("Uncertainty in the easting value")
for key in keys(dict_uncertainty_in_easting)
    println(key," => ",dict_uncertainty_in_easting[key])
end
println("")

println("Uncertainty in the northing value")
for key in keys(dict_uncertainty_in_northing)
    println(key," => ",dict_uncertainty_in_northing[key])
end

Now we display the screen output of the program

program start
The coordinates for :
Shore is [325335.78,5795975.09]
Pier is [325164.29,5795832.74]

Bearing from shore to cerebus/target is 263.3
Bearing from pier to cerebus/target is 305.8

The estimated position of target is [325018.58,5795937.83]
In coordinate format of [Easting value,Northing value]

Now we do the same thing again using measurement

The estimated position of target using Measurements
package is [325018.58 ± 10.32,5795937.83 ± 5.04]

Uncertainty in the easting value
(5.79583274e6, 5.0, 0x0000000000000004) => 5.961637475221173
(263.3, 0.5, 0x0000000000000005) => 3.345995971168122
(305.8, 0.5, 0x0000000000000006) => 2.304707563592274
(5.79597509e6, 5.0, 0x0000000000000002) => 5.961637475221173
(325164.29, 5.0, 0x0000000000000003) => 4.299668552553106
(325335.78, 5.0, 0x0000000000000001) => 0.7003314474468948

Uncertainty in the northing value
(5.79583274e6, 5.0, 0x0000000000000004) => 0.7003314474468948
(263.3, 0.5, 0x0000000000000005) => 2.4132084035663866
(5.79597509e6, 5.0, 0x0000000000000002) => 4.299668552553106
(305.8, 0.5, 0x0000000000000006) => 0.27074091483518714
(325335.78, 5.0, 0x0000000000000001) => 0.5050949698748158
(325164.29, 5.0, 0x0000000000000003) => 0.5050949698748158
julia>

Next we will talk about each section of the Julia program to explain how it is accomplished.

Crosspoint_print()

Now we will talk about the function Crosspoint_print() and Crosspoingg_println()

using Printf

#=  Crosspoint_print(x)
    Print coordinates to two decimal places
=#
function Crosspoint_print(x)
    if length(x) == 1
        print(rstrip(@sprintf("%-10.2f", x)))
        return
    end
    local flag = false
    print("[")
    for ele in x
        print(flag==true ? "," : "")
        print(rstrip(@sprintf("%-10.2f", ele)))
        flag = true
    end
    print("]")
end

Crosspoint_println(x) = begin Crosspoint_print(x); println(); end

Crosspoint_print() will printout a floating point number (or a tuple or an array of floating point number) to two decimal places.

For example:

julia> Crosspoint_print(pi)
3.14
julia> Crosspoint_print([pi,2*pi])
[3.14,6.28]
julia> Crosspoint_print([pi,2*pi,pi^2])
[3.14,6.28,9.87]
julia> Crosspoint_print((pi,2*pi,pi^2))
[3.14,6.28,9.87]

The purposes of Crosspoint_print() is to easily print out the coordinates. For example:

julia> Pos_Shore = (325335.78,5795975.09)  
(325335.78, 5.79597509e6)

julia> print(  Pos_Shore  )
(325335.78, 5.79597509e6)

julia> Crosspoint_print(  Pos_Shore  )
[325335.78,5795975.09]

intersect2lines() and triangulation()

Now we will talk about the function intersect2lines() and triangulation()

#=  intersect2lines(A,a,B,b)

    Return the intersection of two lines from
    two known points A,B using sin_degree() and
    cos_degree() only. a and b are the bearings
    from the unknown position to the two known points A,B
=#
function intersect2lines(A,a,B,b)
    den = sind(a - b)    # sind() function is Sine in degrees, cosd() is similar
    tempX = B[1]*cosd(b)*sind(a) - (A[1]*cosd(a) + (B[2] - A[2])*sind(a))*sind(b)
    tempY = A[2]*cosd(b)*sind(a) - cosd(a)*((A[1] - B[1])*cosd(b) + B[2]*sind(b))
    newX = tempX/den
    newY = tempY/den
    return [newX,newY]
end

#=  triangulation(A,a,B,b)

    Return the position of an unknown location
    by triagulation from two known points A,B and
    the bearing from those two known points A,B to
    the unknown location
=#
triangulation(A,a,B,b) = intersect2lines(A,a+180.0,B,b+180.0)

In trigonometry and geometry, triangulation is the process of determining the location of a point by forming triangles to the point from two (or more) known points.

So the question is "How do we develop triangulation?"

We start from FIRST PRINCIPLES

Next we assume we knew the equation of the line (aka Y = M * X + C ) and work out the location of the unknown location.

Now to find the value of X and Y, both X & Y at point P , all we have to do is to find the values of

So now we know that

M1 = Cot(bearing to flagpole A)
M2 = Cot(bearing to flagpole B)
C1 = Y_flagpoleA - Cot(bearing to flagpole A) * X_flagpoleA
C2 = Y_flagpoleB - Cot(bearing to flagpole B) * X_flagpoleB

Next we put together everything that we know

The value for X is

The value for Y is

And so in SUMMARY, the location of unknown location is

Which is what is in our function intersect2lines(A,a,B,b)

function intersect2lines(A,a,B,b)
    den = sind(a - b)    # sind() function is Sine in degrees, cosd() is similar
    tempX = B[1]*cosd(b)*sind(a) - (A[1]*cosd(a) + (B[2] - A[2])*sind(a))*sind(b)
    tempY = A[2]*cosd(b)*sind(a) - cosd(a)*((A[1] - B[1])*cosd(b) + B[2]*sind(b))
    newX = tempX/den
    newY = tempY/den
    return [newX,newY]
end

Regarding the function triangulation()

The difference between intersect2lines() and triangulation() is very simple.

The function intersect2lines() assumes that you are in location P and you have the following two bearings:

bearing a from location P to Flagpole A
bearing b from location P to Flagpole B

Whereas the function triangulation() assumes that you are
at the two Flagpoles and you have the following two bearings:

bearing a from Flagpole A to location P
bearing b from Flagpole B to location P

So you see the only difference is 180 degrees so

triangulation(A,a,B,b) = intersect2lines(A,a+180.0,B,b+180.0)

Location of Cerberus

So with the Julia code

target_pos_estimate = triangulation(Pos_Shore,
                                    Pos_Shore_bearing_to_target,
                                    Pos_Pier,
                                    Pos_Pier_bearing_to_target)

print("The estimated position of target is ")
Crosspoint_println(target_pos_estimate)
println("In coordinate format of [Easting value,Northing value]")

We can now get the location of the Cerberus

The estimated position of target is [325018.58,5795937.83]
In coordinate format of [Easting value,Northing value]

Once you got the location of Cerberus in UTM format, you can convert it back to Latitude and Longitude using this website

https://awsm-tools.com/utm-to-lat-long

Now we do the same thing again but this time with Measurements

Measurements is a Julia package that allows us to include (the amount of) uncertainties in our input values.

To input uncertainties, we type '\' , 'p' , 'm' , <TAB> on the REPL or a Julia Editor.

#
# Now use Measurement to measure the uncertainties
#    help?>  ±
#    "±" can be typed by \pm<tab>

using Measurements

Shore = (325335.78 ± 5.0,5795975.09 ± 5.0)
Pier  = (325164.29 ± 5.0,5795832.74 ± 5.0)
Shore_bearing_to_target = 263.3 ± 0.5
Pier_bearing_to_target  = 305.8 ± 0.5

For example: Shore_bearing_to_target = 263.3 ± 0.5

The variable "Shore_bearing_to_target" is a measurement.
The measurement "Shore_bearing_to_target" has two components.

The first component is called "value"
The second component is called "uncertainty"

This function returns the "value" of a measurement

Measurements.value()

This other function returns the "uncertainty" of a measurement

Measurements.uncertainty()

I have the following two functions to return the numeric value of the two above components in a string

value_two_decimal_places(m)       = rstrip(  @sprintf("%-16.2f",Measurements.value(m))  )
uncertainty_two_decimal_places(m) = rstrip(  @sprintf("%-16.2f",Measurements.uncertainty(m))  )

The location with uncertainties

So here is the location with uncertainties

target_pos_estimate2 = triangulation(Shore,
                                     Shore_bearing_to_target,
                                     Pier,
                                     Pier_bearing_to_target)

print("The estimated position of target using Measurements package is ")
print("[")
print(value_two_decimal_places(target_pos_estimate2[1]))
print(" ± ")
print(uncertainty_two_decimal_places(target_pos_estimate2[1]))
print(",")
print(value_two_decimal_places(target_pos_estimate2[2]))
print(" ± ")
print(uncertainty_two_decimal_places(target_pos_estimate2[2]))
println("]")
println("")

Here is the output

Now we do the same thing again using measurement

The estimated position of target using Measurements
package is [325018.58 ± 10.32,5795937.83 ± 5.04]

Relative contributions to the uncertainties

dict_uncertainty_in_easting  = Measurements.uncertainty_components(target_pos_estimate2[1])
dict_uncertainty_in_northing = Measurements.uncertainty_components(target_pos_estimate2[2])

println("Uncertainty in the easting value")
for key in keys(dict_uncertainty_in_easting)
    println(key," => ",dict_uncertainty_in_easting[key])
end
println("")

println("Uncertainty in the northing value")
for key in keys(dict_uncertainty_in_northing)
    println(key," => ",dict_uncertainty_in_northing[key])
end

Consider the uncertainties for the "northing" of Cerberus

5795937.83 ± 5.04

The uncertainty for the northing is 5.04 but where does it comes from?

It came from the following input measurements

Shore Easting
Shore Northing
Shore_bearing_to_target
Pier Easting
Pier Northing
Pier_bearing_to_target

There is a way to determine the "relative" contribution to the final uncertainty.

Uncertainty in the northing value
(5.79583274e6, 5.0, 0x0000000000000004) => 0.7003314474468948
(263.3, 0.5, 0x0000000000000005) => 2.4132084035663866
(5.79597509e6, 5.0, 0x0000000000000002) => 4.299668552553106
(305.8, 0.5, 0x0000000000000006) => 0.27074091483518714
(325335.78, 5.0, 0x0000000000000001) => 0.5050949698748158
(325164.29, 5.0, 0x0000000000000003) => 0.5050949698748158

Here we see that the value 5795975.09 (Which is Shore Northing) contributes 4.29966 of relative weight to the final uncertainty of the northing value of Cerberus and this INPUT parameter is THE BIGGEST contribution.

This means that if you reduce the uncertainty for Shore Northing, you will get the biggest bang per buck for lowering the final uncertainty for "northing" of the Cerberus.

This mean that if you want to get the biggest bang per buck for lowering the final uncertainty for "northing" of the Cerberus, you should put all your effort into reducing the uncertainty for Shore Northing.

Convolution in the form of a binary infix operator

Steven Siew — Fri, 16 Dec 2022 15:21:09 +0000

Basic introduction to Convolution

Here is a quick mathematic problem:

What is the result of the following polynomial multiplication?

(1 + 2*X + 3*X^2) * (4 + 5*X + 6*X^2 + 7*X^3 + 8*X^4 +9*X^5)

well the answer is (using mathematica)

4 + 13*X + 28*X^2 + 34*X^3 + 40*X^4 + 46*X^5 + 42*X^6 + 27*X^7

You can obtain the answer by using the convolution function in Julia. Unfortunately, the convolution function is NOT in base Julia but is found in the DSP.jl package

So after you have added DSP to your manifest using package manager, you can do the following.

julia> import DSP: conv

julia> conv([1,2,3],[4,5,6,7,8,9])
8-element Vector{Int64}:
  4
 13
 28
 34
 40
 46
 42
 27

You can perform the convolution as a binary infix function if you define the operator \odot as the convolution function.

# help?>  ⊙
# "⊙" can be typed by \odot<tab>

Now you can do this

julia> ⊙(x,y) = conv(x,y)
⊙ (generic function with 1 method)

julia> [1,2,3] ⊙ [4,5,6,7,8,9]
8-element Vector{Int64}:
  4
 13
 28
 34
 40
 46
 42
 27

Link to paper about convolution in Julia

https://arxiv.org/pdf/1612.08825.pdf

Link to discussion about user defined binary infix operator

3blue1brown explaination of what convolution is

Here is the link to the youtube video about convolution by 3blue1brown. He explains what a convolution really is.

How to use convolution to predict the outcome of an election for a political party

Suppose the Lion political party has three candidates (Adam,Bob and Charles) in three different electorates in an election such that

Adam has 80% chance of winning his seat
Bob has 70% chance of winning his seat
Charles has 60% chance of winning his seat

What are the chances of the Lion political party winning 0,1,2,3 seats?

julia> adam = [(1.0 - 0.8),0.8]
2-element Vector{Float64}:
 0.19999999999999996
 0.8

julia> bob = [(1.0 - 0.7),0.7]
2-element Vector{Float64}:
 0.30000000000000004
 0.7

julia> charles = [(1.0-0.6),0.6]
2-element Vector{Float64}:
 0.4
 0.6

julia> adam ⊙ bob ⊙ charles
4-element Vector{Float64}:
 0.02400000000000002
 0.18800000000000006
 0.45199999999999996
 0.33599999999999997

The answer is

2.4 % chance of winning no seats
18.8 % chance of winning only one seats
45.2 % chance of winning only two seats
33.6 % chance of winning all three seats

Using convolution to do ordinary multiplication

julia> 1234 * 5678
7006652

'''
(4 + 3 * 10^1 + 2 * 10^2 + 1 * 10^3) * (8 + 7 * 10^1 + 6 * 10^2 + 5 * 10^3)
'''

julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
 32
 52
 61
 60
 34
 16
  5

julia> for k = 2:length(result)
           result[k] = result[k] + (result[k-1] ÷ 10)
           result[k-1] = result[k-1] % 10
       end

julia> result
7-element Vector{Int64}:
 2
 5
 6
 6
 0
 0
 7

julia> result = reverse(result)
7-element Vector{Int64}:
 7
 0
 0
 6
 6
 5
 2

julia> string.(result) |> join
"7006652"

You may think that it is a bit silly to use convolution to do multiplication but if you do not have access to BigInt then your only choice left is to use convolution to multiple two gigantic integers together.

Implementing convolution using Fast Fourier Transform

We can implement convolution using fast fourier transform as this article below shows

julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
 32
 52
 61
 60
 34
 16
  5

# Now we do the same thing using FFT

julia> using FFTW

julia> a=Float64.([4,3,2,1])
4-element Vector{Float64}:
 4.0
 3.0
 2.0
 1.0

julia> b=Float64.([8,7,6,5])
4-element Vector{Float64}:
 8.0
 7.0
 6.0
 5.0

julia> len_final = length(a) + length(b) - 1
7

julia> extra_padding_for_a = len_final - length(a)
3

# Pad a with some extra paddings of zeros
julia> aa = vcat(a,[0,0,0])
7-element Vector{Float64}:
 4.0
 3.0
 2.0
 1.0
 0.0
 0.0
 0.0

julia> extra_padding_for_b = len_final - length(b)
3

# Pad b with some extra paddings of zeros
julia> bb = vcat(b,[0,0,0])
7-element Vector{Float64}:
 8.0
 7.0
 6.0
 5.0
 0.0
 0.0
 0.0

julia> a_fft = fft(aa)
7-element Vector{ComplexF64}:
               10.0 + 0.0im
  4.524458669761152 - 4.729234010885294im
 2.1539892641849523 - 1.275184775842325im
 2.3215520660538953 - 0.7129161645984384im
 2.3215520660538953 + 0.7129161645984384im
 2.1539892641849523 + 1.275184775842325im
  4.524458669761152 + 4.729234010885294im

julia> b_fft = fft(bb)
7-element Vector{ComplexF64}:
              26.0 + 0.0im
 6.524458669761152 - 13.491806545954942im
 4.153989264184952 - 0.31203553822726793im
 4.321552066053895 - 3.2208368399238467im
 4.321552066053895 + 3.2208368399238467im
 4.153989264184952 + 0.31203553822726793im
 6.524458669761152 + 13.491806545954942im

julia> product_fft = a_fft .* b_fft
7-element Vector{ComplexF64}:
             260.0 + 0.0im
 -34.2862667915158 - 91.89881294123597im
  8.54974531072476 - 5.969225068086821im
 7.736521480791037 - 10.558244744191306im
 7.736521480791037 + 10.558244744191306im
  8.54974531072476 + 5.969225068086821im
 -34.2862667915158 + 91.89881294123597im

julia> ff = ifft(product_fft)
7-element Vector{ComplexF64}:
              32.0 + 0.0im
              52.0 + 0.0im
              61.0 + 0.0im
 60.00000000000001 + 0.0im
 34.00000000000001 + 0.0im
              16.0 + 0.0im
               5.0 + 0.0im

julia> final_result = real.(ff)
7-element Vector{Float64}:
 32.0
 52.0
 61.0
 60.00000000000001
 34.00000000000001
 16.0
  5.0

# compare with
julia> result = [4,3,2,1] ⊙ [8,7,6,5]
7-element Vector{Int64}:
 32
 52
 61
 60
 34
 16
  5

Now put them together to perform integer multiplication using FFT

function fft_multiplication(str_a,str_b)
    local arr_a, arr_b, len_final
    local extra_padding_for_a
    arr_a = Float64.( reverse(map(x->parse(Int64,x),split(str_a,""))) )
    arr_b = Float64.( reverse(map(x->parse(Int64,x),split(str_b,""))) )
    len_final = length(arr_a) + length(arr_b) - 1
    extra_padding_for_a = len_final - length(arr_a)
    extra_padding_for_b = len_final - length(arr_b)
    arr_a = vcat(arr_a,[ 0.0 for k in 1:extra_padding_for_a ])
    arr_b = vcat(arr_b,[ 0.0 for k in 1:extra_padding_for_b ])
    local a_fft = fft(arr_a)
    local b_fft = fft(arr_b)
    local product_fft = a_fft .* b_fft
    local complex_temp = ifft(product_fft)
    local result = real.(complex_temp)
    result = Int64.(  map(x->round(x,RoundNearest),result)  )
    local k
    for k = 2:length(result)
           result[k] = result[k] + (result[k-1] ÷ 10)
           result[k-1] = result[k-1] % 10
    end
    result = reverse(result)
    result = join(string.(result))
    return result
end

julia> println( fft_multiplication("12345","67890")  )
838102050

# compare with
julia> 12345 * 67890
838102050

Finally the algorithm for a very naive implementation of convolution

How can we talk about convolution without a naive implementation of convolution. This is a very inefficient implementation of an algorithm to calculate convolution.

function naive_convolution(data, kernel)
    local len_data = length(data)
    local len_kernel = length(kernel)
    local len_ans = len_data + len_kernel - 1
    local typeofdata = eltype(data)
    local ans = zeros(typeofdata, len_ans)
    local j,k,acc,loc

    # Next we work out the value for each of the element of the solution
    for j = 1:len_ans
        acc = zero(typeofdata)
        for k = 1:len_kernel
            loc = j - (k - 1)
            if !( loc < 1 || loc > len_data )
                acc += data[loc] * kernel[k]
            end
        end
        ans[j] = acc
    end
    return ans
end

Convolution is the mathematical equivalence of looking for your slippers

Convolution and the Fourier Transform explained visually

What is the Frequency of a Fast Fourier Transform Output?

Where are the magnitude and phase in the output of the FFT?

Why is the output of FFT symmetrical?

Because there are two identical frequencies. One is positive and travel forward in time and the other is negative and travel backwards in time.

Okay, okay. Please do not take me seriously. I am just joking.

Another way to understand Convolution by calculating amount of smoke in the air during an event with fireworks

How to plot Logarithmic Graphs

Steven Siew — Sat, 17 Sep 2022 08:16:49 +0000

Here is a quick julia source code for how to plot Logarithmic graphs

using Plots

cases = [10, 31, 47, 93, 109, 109, 171, 222, 266, 348, 399, 415, 430, 553, 619, 702, 827, 918,
919, 919, 1033, 1110, 1247, 1395, 1561, 1648, 1658, 1999, 2106, 2312, 2477, 2748, 2762, 2766, 2938,
3443, 3725, 4019, 4492, 4499, 4504, 4771, 5401, 5611, 6182, 6640, 6651, 6673, 6947, 7667, 7913,
8271, 9482, 9491, 9507, 9860, 10570, 11094, 11930, 13012, 13027, 13053, 13523, 14607, 15668,
16243, 17236, 17263, 17307, 18263, 19314, 20933, 21488, 22863, 22953, 23018, 24009, 25544, 26264,
27178, 28828, 28848, 28856, 31110, 32373, 33820, 34431, 35809, 35916, 35998, 36949, 38714,
40022, 40855, 42039, 42244, 42342, 43691, 44832, 45911, 47114, 48304, 48439, 48446, 49814,
50813, 51726, 52318, 53241, 53522, 53640, 54186, 55269, 56350, 56686, 57682, 57839, 57839,
58046, 59597, 60291, 60835]


function getLogTicks(vector,count=1)
    function poststr(expnum)
        n = Int64(expnum) ÷ 3
        if n == 0
            return ""
        elseif n == 1
            return "k"
        else
            return "×\$10^{$(n*3)}\$"
        end
    end
    min = ceil(log10(minimum(vector)))
    max = floor(log10(maximum(vector)))
    major = [ [k*n for k = 1:count] for n in 10 .^ collect(min:max) ]
    major = hcat(major...)
    major = vec(major)
    major = filter(x->x<=maximum(vector),major)
    #majorText = ["\$10^{$(round(Int64,i))}\$" for i=min:max]
    majorText = ["$(Int64(j*10^(i%3)))$(poststr(i))" for i=min:max for j=1:count]
    majorText = majorText[1:length(major)]
    minor = [j*10^i for i=(min-1):(max+1) for j=(count+1):9]
    minor = minor[findall(minimum(vector) .<= minor .<= maximum(vector))]
    ([major; minor], [majorText; fill("", length(minor))])
end

#
# Now Plot the cases 
# using LOGSCALE for the Y-Axis
# We can't take log(0) since it is Minus Infinity
# So we assume any value less than one is one

yticks_scale = getLogTicks( [ k<1.0 ? 1.0 : k    for k in cases ], 6 );

plot(0:length(cases)-1,
    [ k<1.0 ? 1.0 : k    for k in cases ],
    title="Monkeypox Cases",
    linewidth=2, ylabel="Cases",
    yaxis=:log10, yticks=yticks_scale,
    xlabel="days",
    legend=:topleft,label="Monkeypox",
    markershape=:circle,size=(800,700)) |> display

And here is the result.

The secret is the function getLogTicks which takes two arguements. The first is the numeric data in a vector format and the second argument is number of text labels on the y-axis.