cma.utilities.math

eigendecomposition of a symmetric matrix, much slower than numpy.linalg.eigh, return (EVals, Basis), the eigenvalues and an orthonormal basis of the corresponding eigenvectors, where

Basis[i]

the i-th row of Basis

columns of Basis, [Basis[j][i] for j in range(len(Basis))]

the i-th eigenvector with eigenvalue EVals[i]

return geometric standard deviation of vals.

The gsd is invariant under linear scaling and independent of the choice of the log-exp base.

kwargs are passed to np.std, in particular ddof.

Hessian estimate for f at x0

remove np type of a number, return int or float.

Return the unchanged argument iff isinstance(float_, int).

Remark that np.float64 is an instance of float however np.int64 is not an instance if int (for good reasons). That is, isinstance(np.array([1])[0], int) is False while isinstance(np.array([1.])[0], float) is True, as of numpy 2.x in 2025.

make print([np.float64(3.3)]) look nicer

rolling average without biasing boundary effects.

The first entries give the average over all first values (until the window width is reached).

If w is not an integer, expontential smoothing with weights proportionate to (1 - 1/w)**i summing to one is executed, thereby putting about 1 - exp(-1) ≈ 0.63 of the weight sum on the last w entries.

Details: the average is mainly based on np.convolve, whereas exponential smoothing is for the time being numerically inefficient and scales quadratically with the length of x.

return an approximation of the inverse cdf value of a standard normal distribution,

assuming 0 < p <= 1/2. This is the "sigma" for which the tail distribution has probability p (AKA quantile or percentile point function).

For p=1/2 we have sigma = 0. The approximation 0.79 + 1.49 x sqrt(-log(p)) + alpha x sqrt(p), where alpha=0.637... is such that p=1/2 maps to sigma=0, has a sigma error within [-0.02, 0.029] for p >= 1e-12 and for p >= 1e-22 with an additional correction as applied. The error is less than 0.0121 for all p <= 0.1. The relative sigma error for p->1/2 and hence sigma->0 is < 0.11.

The input p may be an np.array.

The following is the Python code to assess the accuracy:

%pylab
import scipy
from scipy import stats
import cma

normal_ppf = cma.utilities.math.normal_ppf

pp = np.logspace(-15, np.log10(0.5), 5400)  # target tail probabilities

if 1 < 3:  # sigma vs p
    figure(53)
    gcf().clear()
    grid(True, which='both')
    xlabel('tail probability')
    ylabel('sigma')
    semilogx(pp, normal_ppf(pp), label='approximation')
    semilogx(pp, scipy.stats.norm.ppf(pp), label='true')
    legend()
    # semilogx(pp, -sqrt(7 * log(1/pp/4 + pp)) / 2)
    # semilogx(pp, -scipy.stats.norm.ppf(pp/2))
    # semilogx(pp, sqrt(7 * log((1/pp + pp) / 2)) / 2)
if 1 < 3:  # Delta sigma vs p
    figure(54)
    gcf().clear()
    grid(True, which='both')
    xlabel('tail probability')
    ylabel('sigma difference (error)')
    semilogx(pp, normal_ppf(pp) - scipy.stats.norm.ppf(pp),
            label='absolute')
    abs_error = max(np.abs(normal_ppf(pp) - scipy.stats.norm.ppf(pp)))
    semilogx(pp, (normal_ppf(pp) - scipy.stats.norm.ppf(pp)) / np.maximum(1e-12, -scipy.stats.norm.ppf(pp)),
            label='relative')
    rel_error = max((normal_ppf(pp) - scipy.stats.norm.ppf(pp))
                      / np.maximum(1e-12, -scipy.stats.norm.ppf(pp)))
    ylim(-max(np.abs(ylim())), max(np.abs(ylim())))
    text(xlim()[0], ylim()[0] * 0.98,
        ' the largest absolute and relative errors are {} and {}'
        .format(np.round(abs_error, 6), np.round(rel_error, 6)))
    legend()
if 11 < 3:
    figure(55)
    gcf().clear()
    grid(True, which='both')
    xlabel('tail probability')
    ylabel('sigma ratio')
    semilogx(pp, normal_ppf(pp) / scipy.stats.norm.ppf(pp))
if 11 < 3:
    figure(56)
    gcf().clear()
    ylabel('probability ratio')
    xlabel('sigma')
    plot(stats.norm.ppf(pp), stats.norm.cdf(normal_ppf(pp)) / pp)
    grid(True, which='both')
if 1 < 3:  # true p of sigma vs p input
    figure(56)
    gcf().clear()
    ylabel('true probability ratio')
    xlabel('tail probability (input)')
    semilogx(pp, pp / stats.norm.cdf(normal_ppf(pp)))
    grid(True, which='both')

The approximation is by construction exact for p=1/2. For small probabilities (sigma < -7), the probability is quite sensitive: roughly speaking, a Delta sigma of 0.05 / 0.25 / 1 changes the probability by a factor of 1.4 / 3 / 50, respectively.

n iid dim-dimensional vectors with length norm_(vector).

The vectors are uniformly distributed on a hypersphere surface.

CMA-ES diverges with popsize 100 in 15-D without option 'CSA_clip_length_value': [0,0].

>>> from cma.utilities.math import randhss
>>> dim = 3
>>> assert dim - 1e-7 < sum(randhss(1, dim)[0]**2) < dim + 1e-7

n iid vectors uniformly distributed on the hypersphere surface with mixing in of normal distribution, which can be beneficial in smaller dimension.

test (success) frequency nsucc1 of n1 against nsucc2 of n2.

By default, n2 = n1, hence n2 can be omitted in this case. The frequency data nsucc1, n1, nsucc2, n2 and the kwargs are passed to scipy.stats.contingency.chi2_contingency for a chi-square test of independence.

Return the probability that both success "frequencies" stem from the same underlying distribution (H0). More specifically, return the probability (upper bound) to observe the given or a greater discrepancy between the two success frequency data, given the data stem from the same distribution (formally, are independent of the "row" index (1, 2), hence "test of independence").

Example: when testing 3 versus 8 successes observed in 25 trials respectively, we find that p=0.17 under H0:

import cma
cma.utilities.math.testchisquare(3, 25, 8)

0.17207161769462614

Details: testing the successes nsucc1, nsucc2 is equivalent to testing the respective failures n1 - nsucc1, n2 - nsucc2.

With n1 = n2 = 15, we test nsucc1 = [0, 1, ...] against several nsucc2:

import cma
t = cma.utilities.math.testchisquare  # scipy needs to be installed

n = 15
print('nsucc   nfail   p')
for p1 in range(6):
    for delta in [4, 5, 6]:
        p2 = p1 + (p1 > 1) + delta
        print(' {} {} (= {} {}) {:.2}'.format(
                p1, p2, n-p1, n-p2, t(p1, n, p2)))
    print(' ')

nsucc   nfail   p
 0 4 (= 15 11) 0.11
 0 5 (= 15 10) 0.05
 0 6 (= 15 9) 0.022

 1 5 (= 14 10) 0.17
 1 6 (= 14 9) 0.084
 1 7 (= 14 8) 0.039

 2 7 (= 13 8) 0.11
 2 8 (= 13 7) 0.053
 2 9 (= 13 6) 0.023

 3 8 (= 12 7) 0.13
 3 9 (= 12 6) 0.062
 3 10 (= 12 5) 0.027

 4 9 (= 11 6) 0.14
 4 10 (= 11 5) 0.067
 4 11 (= 11 4) 0.028

 5 10 (= 10 5) 0.14
 5 11 (= 10 4) 0.067
 5 12 (= 10 3) 0.027

Testing 1 (of n=15) versus 6 yields p=0.084 (4-th row), testing 1 versus 7 yields p=0.039. The shown p-values are monotonuously increasing(!) with increasing n: with n = 7, the mentioned p-values are smaller (about 0.033 and 0.007, respectively), with n -> infinity they reach about 0.13 and 0.077, respectively.

Remark: under H0, p~U[0,1], hence E ln(p) = -1 and the geometric average p equals 1/e.

return p-value computed with mannwhitneyu from scipy.stats.

testranksum(data1, data2) returns the probability that data1 and data2, or more discrepant data, are generated while P(d1 < d2) = P(d1 > d2) = 1/2 is true (AKA H0).

This function calls scipy.stats.mannwhitneyu passing data and kwargs, setting by default alternative='two-sided' and method='exact' if len(data1) * len(data2) < 55000 else 'auto'. Pass method='auto' when execution speed is an issue.

The minimal data sizes to possibly get p < 1% two-sided are:

import cma
ranksum = cma.utilities.math.testranksum
[(n, 1e-4 * int(0.5 + 1e4 *  # round and remove np type
     ranksum(range(n[0]), range(n[0], n[0] + n[1]))))
        for n in [[5, 5], [4, 6], [3, 9], [2, 19], [1, 200]]]
# n1, n2, p-value
[([5, 5], 0.0079),
 ([4, 6], 0.0095),
 ([3, 9], 0.0091),
 ([2, 19], 0.0095),
 ([1, 200], 0.01)]

The probabilities are exactly doubled with 'alterative='less' passed as argument.

Details: the above threshold for method='exact' can be changed by assigning cma.utilities.math._testranksum_method_exact_threshold.

Remark: under H0, p~U[0,1], hence E ln(p) = -1 and the geometric average p equals 1/e.

change C in place into a correlation matrix, AKA whitening