Covariance Matrix#

Historical context#

Generally, the covariance matrix is a matrix that describes the relationship between the elements of a multivariate normal distribution. The covariance matrix is a square matrix of size [n x n] where n is the number of elements in the multivariate normal distribution. The covariance matrix is a symmetric matrix and is defined as follows:

\[\begin{split}\Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} & \cdots & \sigma_{1n} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} & \cdots & \sigma_{2n} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} & \cdots & \sigma_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sigma_{n1} & \sigma_{n2} & \sigma_{n3} & \cdots & \sigma_{nn} \end{bmatrix}\end{split}\]

The covariance matrix for the material balance sequence is of particular importance for several statistical tests, including SITMUF and GEMUF. It’s impossible to know the exact matrix in practice, but it can be estimated by knowing expected sensor performance.

Recall that the muf sequence is defined as follows:

\[ \boldsymbol{\text{muf}} = \{ \text{muf}_0, \text{muf}_1, ... \text{muf}_n \} \]

With

\[ \text{muf}_i = \sum_{l \in l_0}\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}I_{t,l} - \sum_{l \in l_1}\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}O_{t,l} - \sum_{l \in l_2}(C_{i,l} - C_{i-1, l}) \]

It’s generally assumed that since the error models are normally distributed, individual muf values (i.e., mbp$_i$) and the muf sequence (i.e., muf) will also be normally distributed. Consequently, the muf sequence can be thought of as a multivariate normal distribution such that:

\[ \boldsymbol{\text{muf}} \sim \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma}) \]

The covariance matrix contains the covariance between different material balances in the sequence. For example, consider the entry $\sigma_{2n}^2$ of the covariance matrix below. This term is the variance between material balance $n$ and $2$.

(1)#\[\begin{split} \begin{aligned} \boldsymbol{\Sigma} & = \begin{bmatrix} {\sigma}_{11}^2 & {\sigma}_{12}^2 & \ldots & {\sigma}_{1n}^2 \\ {\sigma}_{21}^2 & {\sigma}_{22}^2 & \ldots & {\sigma}_{2n}^2 \\ \vdots & \vdots & \ddots & \vdots \\ {\sigma}_{n1}^2 & {\sigma}_{n2}^2 & \ldots & {\sigma}_{nn}^2 \\ \end{bmatrix} = \begin{bmatrix} \boldsymbol{\Sigma_{i-1}} & \boldsymbol{\sigma_{i-1}} \\ \boldsymbol{\sigma_{i-1}}^T & \sigma_{i,i} \\ \end{bmatrix} \end{aligned} \end{split}\]

The covariance matrix itself is often calculated using relative standard deviations, similar to the calculation for $\sigma$muf. In fact, the diagonal terms (i.e., $\Sigma_{1,1}, \Sigma_{2,2}, ... $) are the variance of the material balance (the covariance of material balance $i$ with itself is the variance). There’s two key expressions; one for the covariance diagonals (i.e., $\sigma_{x,x}$) and one for the covariance off diagonals (i.e., $\sigma_{x,x^{\prime}}$). These are calculated by using various variance and covariance rules and propagating the terms.

Recall the expression that was derived from $\sigma$muf; which is used for the diagonals for the covariance matrix:

\[\begin{split} \sigma_{i,i}^2 & \approx \sum_{l \in l_0} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}I_{l,t}\right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) + \sum_{l \in l_2} \left( (C_{i,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) \\ & + \sum_{l \in l_2} \left( (C_{i-1,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) \\ & + \sum_{l \in l_1} \left( \left( \int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i} O_{l,t} \right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) \end{split}\]

The off-diagonal is calculated in a similar manner, but has more terms. The off-diagonal is the covariance between material balance $i$ and $j$.

\[\begin{split} \sigma_{i,j}^2 \approx \text{cov} \bigg(& \sum_{l \in l_0}\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}I_{t,l} - \sum_{l \in l_1}\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}O_{t,l} - \sum_{l \in l_2}(C_{i,l} - C_{i-1, l}), \\ & \sum_{l \in l_0}\int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}I_{t,l} - \sum_{l \in l_1}\int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}O_{t,l} - \sum_{l \in l_2}(C_{j,l} - C_{j-1, l}) \bigg)\\ \end{split}\]

Following the same rules used to derive $\sigma$muf leads to the expression for the covariance off diagonal:

\[\begin{split} \sigma_{i,j}^2 \approx & \sum_{l \in l_2} \left( \left( C_{i,l}C_{j,l} + C_{i-1,l}C_{j-1,l} \right) * (\delta_{S,l})^2 \right) \\ - & \sum_{l \in l_2} \left( \left( C_{i,l}C_{j-1,l} \right) * \left( (\delta_{S,l})^2 + P(j-1 == i)*(\delta_{R,l})^2 \right) \right)\\ - & \sum_{l \in l_2} \left( \left( C_{i-1,l}C_{j,l} \right) * \left( (\delta_{S,l})^2 + P(i-1 == j)*(\delta_{R,l})^2 \right) \right)\\ + & \sum_{l \in l_0} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_j}I_{t,l}\right) \left( \int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}I_{t,l} \right) (\delta_{S,l})^2 \right) \\ + & \sum_{l \in l_1} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_j}O_{t,l}\right) \left( \int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}O_{t,l} \right) (\delta_{S,l})^2 \right) \\ \end{split}\]

Where

\[\begin{split} [P] \equiv \begin{cases} 0 & P \text{ is false} \\ 1 & P \text{ is true} \\ \end{cases}\end{split}\]

Implementation#

The covariance matrix is a $NxN$ matrix at the N-th material balance period. Since MAPIT adopts the “NRTA” style calculation, all $NxN$ entries must be updated at each balance period which simulates the arrival of new information. The MAPIT SITMUF calculation is not well vectorized as the $NxN$ must be resized and calculated at each balance. The calculation starts by looping over balance periods and each entry in the covariance matrix at balance $P$:

AuxFunctions.py

    for currentMB in range(1, int(totalMBPs)):
        for j in range(0, currentMB):

The variables for the different times have a different meaning than in the expressions that were defined above. This is for legacy purposes and to improve alignment with the papers. The following table describes the mapping between the derived expressions and associated code:

Model Component	Code Expression
Balance $i$	`I`
Balance $j$	`IPrime`

For simplicity, the diagonal and off-diagonal terms are broken into multiple components.

Diagonal terms#

\[\begin{split} \sigma_{i,i}^2 & \approx \colorbox{#5755aa}{$\sum_{l \in l_0} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}I_{l,t}\right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) $} \\ & + \colorbox{#0074CC}{$\sum_{l \in l_1} \left( \left( \int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i} O_{l,t} \right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} \\ & + \colorbox{#A2782A}{$\sum_{l \in l_2} \left( (C_{i,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} \\ & + \colorbox{#3E8E87}{$\sum_{l \in l_2} \left( (C_{i-1,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} - \colorbox{#CC3300}{$\sum_{l \in l_2} \left( 2C_{i-1,l}C_{i,l}(\delta_{S,l})^2 \right)$} \\\end{split}\]

\[\begin{split} \sigma_{i,i}^2 & \approx \colorbox{#9392C8}{$\sum_{l \in l_0} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i}I_{l,t}\right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right) $} \\ & + \colorbox{#5CB8FF}{$\sum_{l \in l_1} \left( \left( \int_{t=\text{MBP}_{i-1}}^{\text{MBP}_i} O_{l,t} \right)^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} \\ & + \colorbox{#DEBC7D}{$\sum_{l \in l_2} \left( (C_{i,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} \\ & + \colorbox{#8ECDC8}{$\sum_{l \in l_2} \left( (C_{i-1,l})^2 * ((\delta_{R,l})^2 + (\delta_{S,l})^2) \right)$} - \colorbox{#FF855C}{$\sum_{l \in l_2} \left( 2C_{i-1,l}C_{i,l}(\delta_{S,l})^2 \right)$} \\\end{split}\]

Term 1 AuxFunctions.py

Term 1
AuxFunctions.py

                for k in range(len(inputAppliedError)):

                    logicalInterval = np.logical_and(processedInputTimes[k] >= IPrevious_time,processedInputTimes[k] <= I_time).reshape((-1,))  #select the indices for the relevant time

                    if inputTypes[k] == 'continuous':
                        term1 += trapSum(logicalInterval,processedInputTimes[k],inputAppliedError[k]) **2 * (ErrorMatrix[k, 0]**2 + ErrorMatrix[k, 1]**2)
                    elif inputTypes[k] == 'discrete':
                        term1 += (inputAppliedError[k][:, logicalInterval].sum(axis=1)**2 * (ErrorMatrix[k, 0]**2 + ErrorMatrix[k, 1]**2 ))
                    else:
                        raise Exception("inputTypes[j] is not 'continuous' or 'discrete'")                      

Term 2
AuxFunctions.py

                for k in range(len(outputAppliedError)):

                    logicalInterval = np.logical_and(processedOutputTimes[k] >= IPrevious_time,processedOutputTimes[k] <= I_time).reshape((-1,))
                    locMatrixRow = k + len(inputAppliedError) + len(inventoryAppliedError)                      

                    if outputTypes[k] == 'continuous':
                        term2 += trapSum(logicalInterval,processedOutputTimes[k],outputAppliedError[k])**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)
                    elif outputTypes[k] == 'discrete':
                        term2 += (outputAppliedError[k][:, logicalInterval].sum(axis=1)**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2 ))
                    else:
                        raise Exception("outputTypes[j] is not 'continuous' or 'discrete'")

Term 3
AuxFunctions.py

                for k in range(len(inventoryAppliedError)):
                    locMatrixRow = k+len(inputAppliedError)

                    startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - IPrevious_time).argmin()
                    endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                    term3 += inventoryAppliedError[k][:,endIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)


                if j != 0:
                    for k in range(len(inventoryAppliedError)):
                        locMatrixRow = k + len(inputAppliedError)
                        startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) -IPrevious_time).argmin()
                        endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                        term4 += inventoryAppliedError[k][:,startIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)
                        term5 += inventoryAppliedError[k][:,startIdx] * inventoryAppliedError[k][:,endIdx] * ErrorMatrix[locMatrixRow, 1]**2

                covmatrix[:,j,j] = term1 + term2 + term3 + term4 - 2 * term5

Term 4
AuxFunctions.py

                for k in range(len(inventoryAppliedError)):
                    locMatrixRow = k+len(inputAppliedError)

                    startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - IPrevious_time).argmin()
                    endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                    term3 += inventoryAppliedError[k][:,endIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)


                if j != 0:
                    for k in range(len(inventoryAppliedError)):
                        locMatrixRow = k + len(inputAppliedError)
                        startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) -IPrevious_time).argmin()
                        endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                        term4 += inventoryAppliedError[k][:,startIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)
                        term5 += inventoryAppliedError[k][:,startIdx] * inventoryAppliedError[k][:,endIdx] * ErrorMatrix[locMatrixRow, 1]**2

                covmatrix[:,j,j] = term1 + term2 + term3 + term4 - 2 * term5

Term 5
AuxFunctions.py

The factor 2 is included later when the terms are added and assigned to the covariance matrix.#

                for k in range(len(inventoryAppliedError)):
                    locMatrixRow = k+len(inputAppliedError)

                    startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - IPrevious_time).argmin()
                    endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                    term3 += inventoryAppliedError[k][:,endIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)


                if j != 0:
                    for k in range(len(inventoryAppliedError)):
                        locMatrixRow = k + len(inputAppliedError)
                        startIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) -IPrevious_time).argmin()
                        endIdx = np.abs(processedInventoryTimes[k].reshape((-1,)) - I_time).argmin()

                        term4 += inventoryAppliedError[k][:,startIdx]**2 * (ErrorMatrix[locMatrixRow, 0]**2 + ErrorMatrix[locMatrixRow, 1]**2)
                        term5 += inventoryAppliedError[k][:,startIdx] * inventoryAppliedError[k][:,endIdx] * ErrorMatrix[locMatrixRow, 1]**2

                covmatrix[:,j,j] = term1 + term2 + term3 + term4 - 2 * term5

Off-diagonal terms#

\[\begin{split}\sigma_{i,j}^2 \approx & \colorbox{#5755aa}{$ \sum_{l \in l_0} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_j}I_{t,l}\right) \left( \int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}I_{t,l} \right) (\delta_{S,l})^2 \right) $}\\ + & \colorbox{#0074CC}{$ \sum_{l \in l_1} \left( \left(\int_{t=\text{MBP}_{i-1}}^{\text{MBP}_j}O_{t,l}\right) \left( \int_{t=\text{MBP}_{j-1}}^{\text{MBP}_j}O_{t,l} \right) (\delta_{S,l})^2 \right) $} \\ + & \colorbox{#A2782A}{$\sum_{l \in l_2} \left( \left( C_{i,l}C_{j,l} + C_{i-1,l}C_{j-1,l} \right) * (\delta_{S,l})^2 \right)$} \\ - & \colorbox{#3E8E87}{$\sum_{l \in l_2} \left( \left( C_{i,l}C_{j-1,l} \right) * \left( (\delta_{S,l})^2 + P(j-1 == i)*(\delta_{R,l})^2 \right) \right) $}\\ - & \colorbox{#CC3300}{$\sum_{l \in l_2} \left( \left( C_{i-1,l}C_{j,l} \right) * \left( (\delta_{S,l})^2 + P(i-1 == j)*(\delta_{R,l})^2 \right) \right) $}\\\end{split}\]