SHRIMP: Sub WtdAvCorr

SQUID 2.50 Sub: WtdAvCorr

[The following is a generalised function, used and reused by SQUID 2.50 on at least a couple of different occasions during a 'typical' data-reduction. As such, the inputs and outputs of the functions have generalised names. The identities of the specific inputs and outputs are defined when the function is called.]

The function calculates the weighted average of a single variable, where the errors of index-adjacent values are correlated. This calculation requires a variance-covariance matrix, which can either by defined as an input, or constructed from an input sigma-rho matrix, as required.

Usage

WtdAvCorr(Values, VarCov, N, MeanVal, SigmaMeanVal, MSWD, Prob, SigRho, Bad)

Mandatory inputs

Values: Vector of length N containing values of a single variable, for which the errors of index-adjacent values are correlated.

VarCov: Array of size N x N defining EITHER the variance-covariance matrix OR the sigma-rho matrix associated with Values.

N: Integer defining the number of data-points submitted to the calculation.

MeanVal: The error-weighted mean of Values.

SigmaMeanVal: The 1-sigma absolute error of the MeanVal.

MSWD: The classically calculated Mean Square of Weighted Deviates associated with the MeanVal.

Prob: The classically calculated probability of fit associated with the MSWD and (N - 1) degrees of freedom.

SigRho: Boolean, specifying whether the input VarCov array is a sigma-rho matrix (TRUE) or a variance-covariance matrix (FALSE).

Bad: Boolean for the purpose of error-trapping during calculations.

---

Definition of variables

Values of type Boolean
Bad, SigRho

Values of type Integer
i, j, N

Values of type Double
MeanVal, SigmaMeanVal, MSWD, Prob, Numer, Denom, OmI, OMij,

Variables of type Variant (produced by matrix transposition/multiplication operations)
Omega, SumWtdResids, TransUnwtdResids, TempMat

Arrays of size N x 1 (addressed [1 to N, 1 to 1] ) comprising values of type Double
UnwtdResids

Arrays of size N x N (addressed [1 to N, 1 to N] ) comprising values of type Double
OmegaInv

---

WtdAvCorr

The first step is to construct the variance-covariance matrix. If the input VarCov is a variance-covariance matrix, its elements are written directly to OmegaInv. If the input VarCov is a sigma-rho matrix, the sigmas and rhos are converted to variances and covariances before writing to Omega Inv. Finally, Bad is set to TRUE, and reverts to FALSE if the remaining calculations in WtdAvCorr are completed without incident.

For i = 1 To N  
  For j = 1 To N  
    If SigRho = False  
      OmI = VarCov[i, j]  
    ElseIf i = j --convert sigmas to variances  
      OmI = VarCov[i, i] ^ 2  
    Else  --convert rhos to covariances  
      OmI = VarCov[i, j] * VarCov[j, j] * VarCov[i, i]  
    End If  

  OmegaInv[i, j] = OmI  
  Next j  
Next i  

Bad = TRUE

Omega is generated by inverting the matrix OmegaInv. Ludwig does this using Excel's MINVERSE function (https://msdn.microsoft.com/en-us/library/microsoft.office.interop.excel.worksheetfunction.minverse.aspx).

Omega = MINVERSE(OmegaInv)  

If [Omega gives error]  
  Exit Sub  
End If

Use the elements of Omega to calculate the sums for Numer and Denom:

For i = 1 To N  
  For j = 1 To N  
    If N > 1 --because if N = 1 the matrix inverse has only one dimension  
      OMij = Omega[i, j]      
    Else  
      OMij = Omega[1]  
    End If  
    
    Numer = Numer + (Values[i] + Values[j]) * OMij  
    Denom = Denom + OMij  
  Next j  
Next i  

Denominator error-trapping, evaluation of MeanVal and SigmaMeanVal, and generation of row-vector UnwtdResids for the purpose of calculating MSWD and Prob:

If Denom <= 0  
  Exit Sub  
End If  
  
MeanVal = Numer / Denom / 2
SigmaMeanVal = sqrt(1 / Denom)

For i = 1 To N
  UnwtdResids[i, 1] = Values[i] - MeanVal
Next i

TransUnwtdResids, TempMat, and SumWtdResids are generated from Omega via matrix transposition and multiplication operations. Ludwig used the Excel-implemented functions TRANSPOSE (trivial: https://msdn.microsoft.com/en-us/library/office/ff196261.aspx) and MMULT (https://msdn.microsoft.com/en-us/library/office/ff835849.aspx), and the (legacy) Excel function FDIST (https://support.office.com/en-US/article/FDIST-function-ECF76FBA-B3F1-4E7D-A57E-6A5B7460B786) to calculate Prob:

TransUnwtdResids = TRANSPOSE(UnwtdResids)
TempMat = MMULT(TransUnwtdResids, Omega)
SumWtdResids = MMULT(TempMat, UnwtdResids)

If N = 1 --added "If N = 1" part to avoid divide-by-zero
  MSWD = 0  
  Prob = 0  
Else  
  MSWD = SumWtdResids[1] / (N - 1)  
  Prob = FDIST(MSWD, (N - 1), 1E+9)  
End If  
  
Bad = FALSE  
  
End Sub  

---