-
Notifications
You must be signed in to change notification settings - Fork 16
SHRIMP: Sub SimpleWtdAv
This subroutine calculates a simple, inverse-variance weighted mean from vectors X and Sigma (where Sigma is 1sigma ABSOLUTE) comprising double-precision values. This subroutine is not intended for independent implementation if it can be avoided; instead, this documentation is supplied purely in order to make clear the functionality implied by the invocation of Ludwig's "SimpleWtdAv" function within the parent subroutine FindCoherentGroup.
Squid3 ought to invoke an existing 'Ludwig Library' function (such as WeightedAverage?) in order to perform the same arithmetic as SimpleWtdAv. The documentation below is intended to assist in determining exactly which function should be used for that purpose, and what its arguments should be.
SimpleWtdAv(N, X, Sigma, MeanVal, SigmaMean, MSWD, Probfit)
N: Integer specifying the number of data-rows passed to this subroutine.
X: Vector of length N, containing double-precision ages/values, for all the data-rows passed to this subroutine.
Sigma: Vector of length N, containing double-precision uncertainties (1sigma ABSOLUTE) associated with the ages/values, for all the data-rows passed to this subroutine.
MeanVal: Double-precision value of the weighted mean defined by X and Sigma.
SigmaMean: Double-precision 1sigma ABSOLUTE uncertainty on the value of the weighted mean as defined by non-zero X and Sigma.
MSWD: Double-precision value of the MSWD associated with the weighted mean as defined by non-zero X and Sigma.
Probfit: Double-precision value of the probability of equivalence associated with the weighted mean as defined by non-zero X and Sigma.
Values of type Integer
df, i, Nn
Values of type Double
m, s, sw, Sx, w
For i = 1 To N
If ( IsNumeric( X[i] ) = TRUE ) AND ( X[i] <> 0 )
Nn = 1 + Nn
w = 1 / Sigma[i] ^ 2
Sx = Sx + ( w * X[i] )
sw = sw + w
End If
Next i
If Nn = 0
MeanVal = 0
Exit Sub
ElseIf Nn = 1
MeanVal = Sx / w
m = 0
SigmaMean = SQRT( 1 / sw )
Probfit = 1
Else
df = Nn - 1
MeanVal = Sx / sw
SigmaMean = SQRT( 1 / sw )
s = 0
For i = 1 To N
If ( IsNumeric( X[i] ) = TRUE ) AND ( X[i] <> 0 )
s = s + ( ( X[i] - MeanVal ) / Sigma[i] ) ^ 2
End If
Next i
m = s / df
End If
MSWD = m
If df > 0
Probfit = FDIST(m, df, 1e9) --FDIST is an Excel function; see below
End If
End Sub
I think we have covered Excel's FDIST function before, but if not, there is plenty of documentation online e.g. https://www.excelfunctions.net/excel-fdist-function.html . (Note that SQUID 2.50 uses the 'inaccurate' 2003-era function FDIST, rather than the improved F.DIST.RT function available in Excel 2010 and newer versions.)