This particular function is organized to take in all of the X data points and all of

the Y data points and operate on these vectors one data point,

one element at a time.

So, the first thing that we do is we enter a loop that iterates

through all of these input data points one at a time and makes

updates to the capacity elements one at a time based on the data

measured from the initial point in time up to the Kth point in time.

So, at the beginning of the loop,

we first update the unscaled recursive parameter values as you learned about,

and then we update the scaled parameter values as well.

The code then computes updates to the total capacity estimates for all

of the data up to and including

this particular timestamp in the loop using each one of the four methods.

We start by computing a capacity estimate using the weighted least-squares method.

This is the simplest of all of the methods where capacity is computed as c2

divided by c1 and the Hessian is computed as two times c1.

So, the capacity and the variance are both

stored in the first column of the output matrices in the row

corresponding to this particular iteration through

the data vectors that were provided as input to this function.

Second, we implement the weighted total least-squares method.

Remember that this method is not recursive

and so instead of using the recursive parameters,

we need to use all of the input data up to and including the data measured at

this iteration inclusive for updating the total capacity.

I first make a new variable X which

contains all of the X data up until this point in time,

I make a new variable Y,

and I make variables SX and SY for the variances on X and Y. I also make a vector G,

which is the forgetting factor raised to different powers,

such that the weighting of the present sample is equal to gamma to the zero or one,

and the weighting of the previous sample is equal to gamma to the one or just gamma,

and an earlier samples are simply gamma to

greater powers which result in smaller values for the waiting constants.

Remember that we need to use a numeric solution to find

the weighted total least-squares answer and we need to initialize this numeric solution.

Here I choose to initialize it to the estimate of capacity that

was produced by the weighted least-square solution just a moment ago.

Then, we enter a Newton-Raphson update where we repeatedly

compute the Jacobian J and the Hessian and H of the cost function

and the new capacity is estimated as

the previous capacity minus J divided by H. Within about five iterations,

this Newton-Raphson method should converge to

a good solution providing that the input data are consistent.

So, we store the capacity estimate,

and the variance estimate,

and the output matrices.

Sometimes, I find that the initial estimates of

the weighted total least-squares method do not converge,

which goes against what we believe to be guaranteed by the theory of the problem.

I believe that is true because the input data are

sometimes not consistent if there's a lot of noise on the data.

So, the signs of X and Y should both be positive if I remember correctly and maybe

one is negative because of the numeric noise that's added.

So, these cases, the WTLS method may not converge.

In those cases, I assign an output of NaN,

which stands for "not a number" to

the capacity estimate and I assign values of zero to the variance,

so that the function that is calling this routine can then understand,

don't use that value of total capacity,

wait until later until the capacity is reliable.