Created
February 17, 2020 12:49
-
-
Save dominikmuller/6dca46d78cf41f2d972e2cba896b8e10 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from numpy import matrix | |
| from numpy import zeros | |
| from numpy import ones | |
| from numpy import sqrt | |
| import numpy | |
| import logging as log | |
| class DimensionError(numpy.linalg.LinAlgError): | |
| pass | |
| class BLUE: | |
| """A simple implementation of the Best Linear Unbiased Estimator, based on | |
| "How to combine correlated estimates of a single physical quantity" | |
| in NIM A, 270, pp 110-117""" | |
| def __init__(self): | |
| # List to store uncertainty sources and there cov matrices | |
| self._vlist = [] | |
| def add_measurement(self, m_list): | |
| """Method to add the measured values. m_list has to give proper | |
| len() and must be convertible to numpy.matrix objects. The length | |
| of m_list defines the dimensionality of the average. Everything | |
| else has to match this!""" | |
| self._n_meas = len(m_list) | |
| self._values = matrix(m_list) | |
| # Now that we know the dimensions, initialise utility and storage | |
| self._unit = matrix(ones(len(m_list))).T | |
| self._ematrix = matrix(zeros((len(m_list), len(m_list)))) | |
| def add_source(self, v_matrix, name=''): | |
| """Add a source of uncertainty represented by the covariance | |
| matrix of the measurements for this source. If name is specified, | |
| the contributions of the uncertainty on the overall uncertainty | |
| can be accessed.""" | |
| # matrix.n_dim returns a tuple of the dimensions, check that each | |
| # is identical to the number of measurements | |
| uncmatrix = matrix(v_matrix) | |
| if False in [self._n_meas == x for x in uncmatrix.shape]: | |
| raise DimensionError( | |
| 'Wrong shape. Expected {0}x{0}, got {1}x{2}'.format( | |
| self._n_meas, *uncmatrix.shape)) | |
| if name: | |
| self._vlist.append((name, uncmatrix)) | |
| self._ematrix += uncmatrix | |
| def do(self): | |
| log.debug('Total covariance matrix: \n {}'.format(self._ematrix)) | |
| self.weights = self._ematrix.I * self._unit | |
| self.weights /= self._unit.T * self._ematrix.I * self._unit | |
| log.debug('Weights: {}'.format(self.weights.T)) | |
| self.mean = (self._values*self.weights).max() | |
| log.debug('Mean: {}'.format(self.mean)) | |
| self.error = sqrt(self.weights.T*self._ematrix*self.weights).max() | |
| log.debug('Total uncertainty: {}'.format(self.error)) | |
| log.debug('Splitting errors based on input:') | |
| for name, mat in self._vlist: | |
| tmp_error = sqrt((self.weights.T*mat*self.weights).max()) | |
| setattr(self, name, tmp_error) | |
| log.debug('{} : {}'.format(name, tmp_error)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment