docstrings

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,85 +111,87 @@ 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
 
 
 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):