diff --git a/uncertainties/core.py b/uncertainties/core.py index f5e14ffb..b0df3395 100644 --- a/uncertainties/core.py +++ b/uncertainties/core.py @@ -74,39 +74,27 @@ def correlated_values(nom_values, covariance_mat, tags=None): """ - Return numbers with uncertainties (AffineScalarFunc objects) - that correctly reproduce the given covariance matrix, and have - the given (float) values as their nominal value. + Return numbers with uncertainties (AffineScalarFunc objects) that have nominal + values given by nom_values and are correlated according to covariance_mat. - The correlated_values_norm() function returns the same result, - but takes a correlation matrix instead of a covariance matrix. + nom_values -- length (N,) sequence of (real) nominal values for the returned numbers + with uncertainties. - The list of values and the covariance matrix must have the - same length, and the matrix must be a square (symmetric) one. + covariance_mat -- (N, N) array representing the target covariance matrix for the + returned numbers with uncertainties. This matrix must be positive semi-definite. - The numbers with uncertainties returned depend on newly - created, independent variables (Variable objects). + tags -- optional length (N,) sequence of strings to use as tags for the variables on + which the returned numbers with uncertainties depend. - nom_values -- sequence with the nominal (real) values of the - numbers with uncertainties to be returned. - - covariance_mat -- full covariance matrix of the returned numbers with - uncertainties. For example, the first element of this matrix is the - variance of the first number with uncertainty. This matrix must be a - NumPy array-like (list of lists, NumPy array, etc.). - - tags -- if 'tags' is not None, it must list the tag of each new - independent variable. - - This function raises NotImplementedError if numpy cannot be - imported. + This function raises NotImplementedError if numpy cannot be imported. """ if numpy is None: msg = ( - "uncertainties was not able to import numpy so " - "correlated_values is unavailable." + "uncertainties was not able to import numpy so correlated_values is " + "unavailable." ) raise NotImplementedError(msg) + covariance_mat = numpy.array(covariance_mat) if not numpy.all(numpy.isreal(covariance_mat)): @@ -123,55 +111,64 @@ def correlated_values(nom_values, covariance_mat, tags=None): if tags is None: tags = [None] * n_dims - ufloat_atoms = numpy.array([Variable(0, 1, tags[dim]) for dim in range(n_dims)]) + X = numpy.array([Variable(0, 1, tags[dim]) for dim in range(n_dims)]) try: """ - Covariance matrices for linearly independent random variables are - symmetric and positive-definite so they can be decomposed sa + Covariance matrices for linearly independent random variables are symmetric and + positive-definite so they can be Cholesky decomposed as + C = L * L.T - with L a lower triangular matrix. - Let R be a vector of independent random variables with zero mean and - unity variance. Then consider - Y = L * R + with L a lower triangular matrix. Let X be a vector of independent random + variables with zero mean and unity variance. Then consider + + Y = L * X + and - Cov(Y) = E[Y * Y.T] = E[L * R * R.T * L.T] = L * E[R * R.t] * L.T - = L * Cov(R) * L.T = L * I * L.T = L * L.T = C - where Cov(R) = I because the random variables in V are independent with - unity variance. So Y defined as above has covariance C. + + Cov(Y) = E[Y * Y.T] = E[L * X * X.T * L.T] = L * E[X * X.t] * L.T + = L * Cov(X) * L.T = L * I * L.T = L * L.T = C + + where Cov(X) = I because the random variables in X are independent with + unity variance. So Y = L * X has covariance C. """ L = numpy.linalg.cholesky(covariance_mat) - Y = L @ ufloat_atoms + Y = L @ X except numpy.linalg.LinAlgError: """" - If two random variables are linearly dependent, e.g. x and y=2*x, then - their covariance matrix will be degenerate. In this case, a Cholesky - decomposition is not possible, but an eigenvalue decomposition is. Even - in this case, covariance matrices are symmetric, so they can be - decomposed as + If two random variables are linearly dependent, e.g. x and y=2*x, then their + covariance matrix will be degenerate. In this case, a Cholesky decomposition is + not possible, but an eigenvalue decomposition is. Even in this case, covariance + matrices are symmetric, so they can be decomposed as - C = U V U^T + C = U D U^T + + with U orthogonal and D diagonal with non-negative (though possibly + zero-valued) entries. Let S = sqrt(D) and + + Y = U * S * X - with U orthogonal and V diagonal with non-negative (though possibly - zero-valued) entries. Let S = sqrt(V) and - Y = U * S * R Then - Cov(Y) = E[Y * Y.T] = E[U * S * R * R.T * S.T * U.T] - = U * S * E[R * R.T] * S.T * U.T + + Cov(Y) = E[Y * Y.T] = E[U * S * X * X.T * S.T * U.T] + = U * S * E[X * X.T] * S.T * U.T = U * S * I * S.T * U.T - = U * S * S.T * U.T = U * V * U.T + = U * S * S.T * U.T = U * D * U.T = C - So Y defined as above has covariance C. + + So Y = U * S * X defined as above has covariance C. """ - eig_vals, eig_vecs = numpy.linalg.eigh(covariance_mat) + eig_vals, U = numpy.linalg.eigh(covariance_mat) """ Eigenvalues may be close to zero but negative. We clip these to zero. """ eig_vals = numpy.clip(eig_vals, a_min=0, a_max=None) - std_devs = numpy.diag(numpy.sqrt(eig_vals)) - Y = numpy.transpose(eig_vecs @ std_devs @ ufloat_atoms) + if numpy.any(eig_vals < 0): + raise ValueError("covariance_mat must be positive semi-definite.") + S = numpy.diag(numpy.sqrt(eig_vals)) + Y = numpy.transpose(U @ S @ X) result = numpy.array(nom_values) + Y return result @@ -179,29 +176,22 @@ def correlated_values(nom_values, covariance_mat, tags=None): def correlated_values_norm(values_with_std_dev, correlation_mat, tags=None): """ - Return correlated values like correlated_values(), but takes - instead as input: - - - nominal (float) values along with their standard deviation, and - - a correlation matrix (i.e. a normalized covariance matrix). + Return numbers with uncertainties (AffineScalarFunc objects) that have nominal + values and standard deviations given by values_with_std_dev and whose correlation + matrix is given by correlation_mat. - values_with_std_dev -- sequence of (nominal value, standard - deviation) pairs. The returned, correlated values have these - nominal values and standard deviations. + values_with_std_dev -- Length (N,) sequence of tuples of (nominal_value, std_dev) + corresponding to the nominal values and standard deviations of the returned numbers + with uncertainty. std_devs must be non-negative. - correlation_mat -- correlation matrix between the given values, except - that any value with a 0 standard deviation must have its correlations - set to 0, with a diagonal element set to an arbitrary value (something - close to 0-1 is recommended, for a better numerical precision). When - no value has a 0 variance, this is the covariance matrix normalized by - standard deviations, and thus a symmetric matrix with ones on its - diagonal. This matrix must be an NumPy array-like (list of lists, - NumPy array, etc.). + correlation_mat -- (N, N) array representing the target normalized correlation + matrix between the returned values with uncertainty. This matrix must be positive + semi-definite. - tags -- like for correlated_values(). + tags -- optional length (N,) sequence of strings to use as tags for the variables on + which the returned numbers with uncertainties depend. - This function raises NotImplementedError if numpy cannot be - imported. + This function raises NotImplementedError if numpy cannot be imported. """ if numpy is None: msg = ( @@ -211,11 +201,21 @@ def correlated_values_norm(values_with_std_dev, correlation_mat, tags=None): raise NotImplementedError(msg) nominal_values, std_devs = tuple(zip(*values_with_std_dev)) + if any(numpy.array(std_devs) < 0): + raise ValueError("All std_devs must be non-negative.") + + """ + The entries Corr[i, j] of the correlation matrix are related to the entries + Cov[i, j] by + Cov[i, j] = Corr[i, j] * s[i] * s[j] + + where s[i] is the standard deviation of random variable i. + """ outer_std_devs = numpy.outer(std_devs, std_devs) cov_mat = correlation_mat * outer_std_devs - return correlated_values(nominal_values, cov_mat) + return correlated_values(nominal_values, cov_mat, tags) def correlation_matrix(nums_with_uncert):