Skip to content

Random Numbers and Probability Distributions

Random number generators, multivariate Gaussian mixtures, k-means clustering, and probability density functions.

Random number generation

v_randvec

V_RANDVEC - Generate random vectors from a GMM distribution.

v_randvec 

v_randvec(
    n, m, c=None, w=None, mode="g"
) -> tuple[ndarray, ndarray]

Generate real or complex GMM/lognormal random vectors.

Parameters:

Name Type Description Default
n int

Number of points to generate.

 required
m array_like

Mixture means, shape (k, p).

 required
c array_like

Covariances: diagonal (k, p), full (p, p, k), or scalar per mixture.

 None
w array_like

Mixture weights (k,). Default: equal weights.

 None
mode str

'g' = real Gaussian (default), 'c' = complex Gaussian, 'l' = lognormal.

 'g'

Returns:

Name Type Description
x ndarray

Output data, shape (n, p).
kx ndarray

Mixture number for each row (0-based), shape (n,).
Source code in pyvoicebox/v_randvec.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def v_randvec(n, m, c=None, w=None, mode='g') -> tuple[np.ndarray, np.ndarray]:
    """Generate real or complex GMM/lognormal random vectors.

    Parameters
    ----------
    n : int
        Number of points to generate.
    m : array_like
        Mixture means, shape (k, p).
    c : array_like, optional
        Covariances: diagonal (k, p), full (p, p, k), or scalar per mixture.
    w : array_like, optional
        Mixture weights (k,). Default: equal weights.
    mode : str, optional
        'g' = real Gaussian (default), 'c' = complex Gaussian, 'l' = lognormal.

    Returns
    -------
    x : ndarray
        Output data, shape (n, p).
    kx : ndarray
        Mixture number for each row (0-based), shape (n,).
    """
    m = np.asarray(m, dtype=float)
    if m.ndim == 1:
        m = m.reshape(1, -1)
    sm = m.shape
    k = sm[0]
    p = sm[1]

    if c is None:
        c = np.ones_like(m)
    c = np.asarray(c, dtype=float)
    sc = c.shape

    fullc = c.ndim > 2 or (c.ndim == 2 and sc[0] > k)

    if w is None:
        w = np.ones(k)
    else:
        w = np.asarray(w, dtype=float).ravel()

    if isinstance(mode, str) and len(mode) > 0:
        ty = mode[0]
    else:
        ty = 'g'

    x = np.zeros((n, p))

    if k > 1:
        kx = v_randiscr(w, n)
    else:
        kx = np.zeros(n, dtype=int)

    for kk_idx in range(k):
        nx = np.where(kx == kk_idx)[0]
        nn = len(nx)
        if nn == 0:
            continue

        mm = m[kk_idx, :]
        if fullc:
            cc = c[:, :, kk_idx].copy()
            if ty == 'l':
                cc = np.log(1 + cc / np.outer(mm, mm))
                mm = np.log(mm) - 0.5 * np.diag(cc)
        else:
            cc = c[kk_idx, :].copy()
            if ty == 'l':
                cc = np.log(1 + cc / mm ** 2)
                mm = np.log(mm) - 0.5 * cc

        if ty == 'c':
            q = np.sqrt(0.5)
            xx = q * np.random.randn(nn, p) + 1j * q * np.random.randn(nn, p)
        else:
            xx = np.random.randn(nn, p)

        if fullc:
            cc_sym = (cc + cc.T) / 2
            eigvals, eigvecs = np.linalg.eigh(cc_sym)
            xx = xx * np.sqrt(np.abs(eigvals))[np.newaxis, :] @ eigvecs.T + mm[np.newaxis, :]
        else:
            xx = xx * np.sqrt(np.abs(cc))[np.newaxis, :] + mm[np.newaxis, :]

        x[nx, :] = xx

    if ty == 'l':
        x = np.exp(x)

    return x, kx

v_randiscr

V_RANDISCR - Generate discrete random values with specified probabilities.

v_randiscr 

v_randiscr(p=None, n=1, a=None) -> ndarray

Generate discrete random numbers with specified probabilities.

Parameters:

Name Type Description Default
p array_like

Probabilities (not necessarily normalized). Default: uniform.

 None
n int

Number of random values to generate. Positive = with replacement, negative = without replacement. Default: 1.

 1
a array_like

Output alphabet. Default: 0-based indices.

 None

Returns:

Name Type Description
x ndarray

Vector of values taken from alphabet a.
Source code in pyvoicebox/v_randiscr.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def v_randiscr(p=None, n=1, a=None) -> np.ndarray:
    """Generate discrete random numbers with specified probabilities.

    Parameters
    ----------
    p : array_like, optional
        Probabilities (not necessarily normalized). Default: uniform.
    n : int, optional
        Number of random values to generate. Positive = with replacement,
        negative = without replacement. Default: 1.
    a : array_like, optional
        Output alphabet. Default: 0-based indices.

    Returns
    -------
    x : ndarray
        Vector of values taken from alphabet a.
    """
    got_a = a is not None

    if p is None or (isinstance(p, np.ndarray) and p.size == 0):
        if got_a:
            a = np.asarray(a)
            p = np.ones(len(a))
        else:
            p = np.ones(2)
            a = np.arange(2)
            got_a = True

    p = np.asarray(p, dtype=float)
    q = p.ravel()
    d = len(q)

    if got_a:
        a = np.asarray(a)
        if d != a.size:
            raise ValueError('sizes of alphabet and probability vector must match')

    if n >= 1:
        # Sample with replacement
        n = int(n)
        cs = np.cumsum(q / np.sum(q))
        cs[-1] = 1.0  # ensure no numerical issues
        r = np.random.rand(n)
        x = np.searchsorted(cs, r)
        # 0-based indices
    else:
        # Sample without replacement
        n = abs(int(n))
        if n > d:
            raise ValueError('Number of output samples exceeds alphabet size')
        if np.all(q == q[0]):
            # Uniform probabilities
            perm = np.random.permutation(d)
            x = perm[:n]
        else:
            cs = np.cumsum(q / np.sum(q))
            cs[-1] = 1.0
            chosen = set()
            x_list = []
            while len(x_list) < n:
                r = np.random.rand()
                idx = np.searchsorted(cs, r)
                if idx not in chosen:
                    chosen.add(idx)
                    x_list.append(idx)
            x = np.array(x_list)

    if got_a:
        x = a.ravel()[x]
    return x

v_stdspectrum

V_STDSPECTRUM - Generate standard acoustic/speech spectra (simplified).

v_stdspectrum 

v_stdspectrum(
    s, m="s", f=8192, n=None, zi=None, bs=None, as_=None
) -> tuple[ndarray, ndarray]

Generate standard acoustic/speech spectra in s- or z-domain.

Simplified implementation supporting s-domain output and basic z-domain conversion for the most common spectrum types.

Parameters:

Name Type Description Default
s int or str

Spectrum type: 1 = White 2 = A-Weight 3 = B-Weight 4 = C-Weight 9 = USASI

 required
m str

Mode: 's' for s-domain, 'z' for z-domain. Default 's'.

 's'
f float

Sample frequency (for 'z' mode). Default 8192.

 8192
n int

Number of samples (for 't' mode).

 None

Returns:

Name Type Description
b ndarray

Numerator coefficients.
a ndarray

Denominator coefficients.
Source code in pyvoicebox/v_stdspectrum.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def v_stdspectrum(s, m='s', f=8192, n=None, zi=None, bs=None, as_=None) -> tuple[np.ndarray, np.ndarray]:
    """Generate standard acoustic/speech spectra in s- or z-domain.

    Simplified implementation supporting s-domain output and basic z-domain
    conversion for the most common spectrum types.

    Parameters
    ----------
    s : int or str
        Spectrum type:
            1 = White
            2 = A-Weight
            3 = B-Weight
            4 = C-Weight
            9 = USASI
    m : str, optional
        Mode: 's' for s-domain, 'z' for z-domain. Default 's'.
    f : float, optional
        Sample frequency (for 'z' mode). Default 8192.
    n : int, optional
        Number of samples (for 't' mode).

    Returns
    -------
    b : ndarray
        Numerator coefficients.
    a : ndarray
        Denominator coefficients.
    """
    # S-domain spectrum definitions (zeros, poles, gain)
    spectra = _get_spectra()

    if isinstance(s, str):
        names = {name.lower(): idx for idx, name in enumerate(spectra.keys())}
        si = names.get(s.lower(), 0)
        sn = s
    else:
        si = int(s)

    spec_list = list(spectra.values())
    if si < 1 or si > len(spec_list):
        if si == 0 and bs is not None:
            sb = np.asarray(bs, dtype=float)
            sa = np.asarray(as_, dtype=float) if as_ is not None else np.array([1.0])
        else:
            raise ValueError(f'Undefined spectrum type: {s}')
    else:
        sb, sa = spec_list[si - 1]

    m1 = m[0] if m else 's'

    if m1 == 's':
        return sb, sa
    elif m1 == 'z':
        if si == 1:  # White noise
            return np.array([1.0]), np.array([1.0])
        # Use bilinear transform
        bz, az = bilinear(sb, sa, f)
        return bz, az
    else:
        return sb, sa

v_randfilt

V_RANDFILT - Generate filtered Gaussian noise without initial transient.

v_randfilt 

v_randfilt(
    pb, pa, ny=0, zi=None
) -> tuple[ndarray, ndarray, ndarray, float]

Generate filtered Gaussian noise without initial transient.

Parameters:

Name Type Description Default
pb array_like

Numerator polynomial of discrete time filter.

 required
pa array_like

Denominator polynomial of discrete time filter.

 required
ny int

Number of output samples required.

 0
zi array_like

Initial filter state.

 None

Returns:

Name Type Description
y ndarray

Filtered Gaussian noise.
zf ndarray

Final filter state.
u ndarray

State covariance factor (u@u' = state covariance).
p float

Expected value of y(i)^2.
Source code in pyvoicebox/v_randfilt.py 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def v_randfilt(pb, pa, ny=0, zi=None) -> tuple[np.ndarray, np.ndarray, np.ndarray, float]:
    """Generate filtered Gaussian noise without initial transient.

    Parameters
    ----------
    pb : array_like
        Numerator polynomial of discrete time filter.
    pa : array_like
        Denominator polynomial of discrete time filter.
    ny : int, optional
        Number of output samples required.
    zi : array_like, optional
        Initial filter state.

    Returns
    -------
    y : ndarray
        Filtered Gaussian noise.
    zf : ndarray
        Final filter state.
    u : ndarray
        State covariance factor (u@u' = state covariance).
    p : float
        Expected value of y(i)^2.
    """
    pb = np.asarray(pb, dtype=float).ravel()
    pa = np.asarray(pa, dtype=float).ravel()

    if pa[0] != 1:
        pb = pb / pa[0]
        pa = pa / pa[0]

    u = None
    p = None

    if zi is None or True:  # always compute for u/p outputs
        lb = len(pb)
        la = len(pa)
        k = max(la, lb) - 1  # filter order
        n = la - 1            # denominator order

        if k == 0:
            if zi is None:
                zi = np.array([])
            u = np.zeros((0, 0))
            p = pb[0] ** 2
        else:
            # Form controllability matrix
            q = np.zeros((k, k))
            _, q[:, 0] = lfilter(pb, pa, [1.0], zi=np.zeros(k))
            for i in range(1, k):
                _, q[:, i] = lfilter(pb, pa, [0.0], zi=q[:, i - 1])

            # Step-down procedure
            s = np.random.randn(k)
            ii = slice(k - n, k)
            if n > 0:
                m_mat = np.zeros((n, n))
                g = pa.copy()
                for i in range(n):
                    denom = np.sqrt((g[0] - g[-1]) * (g[0] + g[-1]))
                    if abs(denom) < 1e-15:
                        denom = 1e-15
                    g = (g[0] * g[:-1] - g[-1] * g[-1:0:-1]) / denom
                    m_mat[i, i:n] = g[:n - i]

                T = np.triu(toeplitz(pa[:n]))
                try:
                    s[ii] = T @ np.linalg.solve(m_mat, s[ii])
                except np.linalg.LinAlgError:
                    pass

                u = q.copy()
                try:
                    u[:, ii] = q[:, ii] @ T @ np.linalg.inv(m_mat)
                except np.linalg.LinAlgError:
                    pass
            else:
                u = q.copy()

            if zi is None:
                zi = q @ s

            p = u[0, :] @ u[0, :] + pb[0] ** 2

    if ny > 0:
        if len(zi) == 0:
            y, zf = lfilter(pb, pa, np.random.randn(ny)), zi
        else:
            y, zf = lfilter(pb, pa, np.random.randn(ny), zi=zi)
    else:
        zf = zi
        y = np.array([])

    return y, zf, u, p

v_rnsubset

V_RNSUBSET - Choose k distinct random integers from 1:n.

v_rnsubset 

v_rnsubset(k, n) -> ndarray

Choose k distinct random integers from 0:n-1.

Note: Returns 0-based indices (unlike MATLAB's 1-based).

Parameters:

Name Type Description Default
k int

Number of distinct integers required.

 required
n int

Range is 0 to n-1.

 required

Returns:

Name Type Description
m ndarray

Array of k distinct random integers.
Source code in pyvoicebox/v_rnsubset.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def v_rnsubset(k, n) -> np.ndarray:
    """Choose k distinct random integers from 0:n-1.

    Note: Returns 0-based indices (unlike MATLAB's 1-based).

    Parameters
    ----------
    k : int
        Number of distinct integers required.
    n : int
        Range is 0 to n-1.

    Returns
    -------
    m : ndarray
        Array of k distinct random integers.
    """
    if k > n:
        raise ValueError('k must be <= n')

    f, e = np.frexp(n)
    if k > 0.03 * n * (e - 1):
        # For large k, random permutation
        m = np.random.permutation(n)
    else:
        # Fisher-Yates partial shuffle
        m = np.arange(n)
        for i in range(k):
            j = i + np.random.randint(n - i)
            m[i], m[j] = m[j], m[i]

    return m[:k]

v_usasi

V_USASI - Generate USASI noise.

v_usasi 

v_usasi(n, fs=8000) -> ndarray

Generate n samples of USASI noise at sample frequency fs.

USASI noise simulates long-term average audio program material.

Parameters:

Name Type Description Default
n int

Number of samples.

 required
fs float

Sample frequency in Hz. Default 8000.

 8000

Returns:

Name Type Description
x ndarray

USASI noise signal.
Source code in pyvoicebox/v_usasi.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def v_usasi(n, fs=8000) -> np.ndarray:
    """Generate n samples of USASI noise at sample frequency fs.

    USASI noise simulates long-term average audio program material.

    Parameters
    ----------
    n : int
        Number of samples.
    fs : float, optional
        Sample frequency in Hz. Default 8000.

    Returns
    -------
    x : ndarray
        USASI noise signal.
    """
    b = np.array([1.0, 0.0, -1.0])
    a = np.real(np.poly(np.exp(-np.array([100, 320]) * 2 * np.pi / fs)))
    x, _, _, _ = v_randfilt(b, a, n)
    return x

Gaussian mixtures and k-means

v_gaussmix

V_GAUSSMIX - Fit a Gaussian mixture model using EM algorithm.

v_gaussmix 

v_gaussmix(
    x, c=None, l=None, m0=None, v0=None, w0=None, wx=None
) -> tuple[
    ndarray,
    ndarray,
    ndarray,
    float,
    float,
    ndarray,
    ndarray,
]

Fit a Gaussian mixture PDF to data using EM algorithm.

Parameters:

Name Type Description Default
x array_like

Input data, shape (n, p).

 required
c float

Minimum variance of normalized data. Default: 1/n^2.

 None
l float

Max iterations (integer part) + stopping threshold (fractional part). Default: 100.0001.

 None
m0 int or array_like

Number of mixtures, or initial means (k, p).

 None
v0 str or array_like

Initialization mode string or initial variances.

 None
w0 array_like

Initial weights or data point weights.

 None
wx array_like

Data point weights when m0/v0/w0 are explicit initial values.

 None

Returns:

Name Type Description
m ndarray

Mixture means, shape (k, p).
v ndarray

Mixture variances, shape (k, p) or (p, p, k).
w ndarray

Mixture weights, shape (k,).
g float

Average log probability.
f float

Fisher's Discriminant.
pp ndarray

Log probability of each data point.
gg ndarray

Average log probabilities at each iteration.
Source code in pyvoicebox/v_gaussmix.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
def v_gaussmix(x, c=None, l=None, m0=None, v0=None, w0=None, wx=None) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float, np.ndarray, np.ndarray]:
    """Fit a Gaussian mixture PDF to data using EM algorithm.

    Parameters
    ----------
    x : array_like
        Input data, shape (n, p).
    c : float, optional
        Minimum variance of normalized data. Default: 1/n^2.
    l : float, optional
        Max iterations (integer part) + stopping threshold (fractional part).
        Default: 100.0001.
    m0 : int or array_like
        Number of mixtures, or initial means (k, p).
    v0 : str or array_like, optional
        Initialization mode string or initial variances.
    w0 : array_like, optional
        Initial weights or data point weights.
    wx : array_like, optional
        Data point weights when m0/v0/w0 are explicit initial values.

    Returns
    -------
    m : ndarray
        Mixture means, shape (k, p).
    v : ndarray
        Mixture variances, shape (k, p) or (p, p, k).
    w : ndarray
        Mixture weights, shape (k,).
    g : float
        Average log probability.
    f : float
        Fisher's Discriminant.
    pp : ndarray
        Log probability of each data point.
    gg : ndarray
        Average log probabilities at each iteration.
    """
    x = np.asarray(x, dtype=float)
    n, p = x.shape
    wn = np.ones(n)

    if c is None:
        c = 1.0 / n ** 2
    else:
        c = float(c)

    fulliv = False  # initial variance is not full

    if l is None:
        l = 100 + 1e-4

    if v0 is None or isinstance(v0, str):
        # No initial values given - use k-means or equivalent
        if v0 is None:
            v0 = 'hf'

        if w0 is not None:
            wx_local = np.asarray(w0, dtype=float).ravel()
        else:
            wx_local = wn.copy()
        wx_local = wx_local / np.sum(wx_local)

        if isinstance(m0, (int, np.integer)):
            k = int(m0)
            has_initial_means = False
        else:
            m0 = np.asarray(m0, dtype=float)
            k = m0.shape[0]
            has_initial_means = 'm' in v0

        fv = 'v' in v0

        mx0 = wx_local @ x
        vx0 = wx_local @ (x ** 2) - mx0 ** 2
        sx0 = np.sqrt(vx0)
        sx0[sx0 == 0] = 1.0

        if n <= k:
            xs = (x - mx0[np.newaxis, :]) / sx0[np.newaxis, :]
            m_fit = xs[np.arange(k) % n, :]
            v_fit = np.zeros((k, p))
            w_fit = np.zeros(k)
            w_fit[:n] = 1.0 / n
            if l > 0:
                l = 0.1
        else:
            if 's' in v0:
                xs = x.copy()
            else:
                xs = (x - mx0[np.newaxis, :]) / sx0[np.newaxis, :]
                if has_initial_means:
                    m0 = (m0 - mx0[np.newaxis, :]) / sx0[np.newaxis, :]

            w_fit = np.full(k, 1.0 / k)

            if 'k' in v0:
                if has_initial_means:
                    m_fit, _, j, _ = v_kmeans(xs, k, m0)
                elif 'p' in v0:
                    m_fit, _, j, _ = v_kmeans(xs, k, 'p')
                else:
                    m_fit, _, j, _ = v_kmeans(xs, k, 'f')
            elif 'h' in v0:
                if has_initial_means:
                    m_fit, _, j, _ = v_kmeanhar(xs, k, None, 4, m0)
                elif 'p' in v0:
                    m_fit, _, j, _ = v_kmeanhar(xs, k, None, 4, 'p')
                else:
                    m_fit, _, j, _ = v_kmeanhar(xs, k, None, 4, 'f')
            elif 'p' in v0:
                j = np.random.randint(0, k, size=n)
                forced = v_rnsubset(k, n)
                j[forced] = np.arange(k)
                m_fit = np.zeros((k, p))
                for i in range(k):
                    mask = j == i
                    if np.any(mask):
                        m_fit[i, :] = np.mean(xs[mask, :], axis=0)
            else:
                if has_initial_means:
                    m_fit = m0.copy()
                else:
                    m_fit = xs[v_rnsubset(k, n), :]
                _, j, _, _ = v_kmeans(xs, k, m_fit, 0)

            if 's' in v0:
                xs = (x - mx0[np.newaxis, :]) / sx0[np.newaxis, :]

            v_fit = np.zeros((k, p))
            w_fit = np.zeros(k)
            for i in range(k):
                mask = j == i
                ni = np.sum(mask)
                w_fit[i] = (ni + 1) / (n + k)
                if ni > 0:
                    v_fit[i, :] = np.sum((xs[mask, :] - m_fit[i, :][np.newaxis, :]) ** 2, axis=0) / ni

        m = m_fit
        v = v_fit
        w = w_fit
    else:
        # Use initial values given as input parameters
        if wx is None:
            wx_local = wn.copy()
        else:
            wx_local = np.asarray(wx, dtype=float).ravel()
        wx_local = wx_local / np.sum(wx_local)

        mx0 = wx_local @ x
        vx0 = wx_local @ (x ** 2) - mx0 ** 2
        sx0 = np.sqrt(vx0)
        sx0[sx0 == 0] = 1.0

        m0 = np.asarray(m0, dtype=float)
        v0_arr = np.asarray(v0, dtype=float)
        w0_arr = np.asarray(w0, dtype=float).ravel()

        k = m0.shape[0]
        xs = (x - mx0[np.newaxis, :]) / sx0[np.newaxis, :]
        m = (m0 - mx0[np.newaxis, :]) / sx0[np.newaxis, :]
        v = v0_arr.copy()
        w = w0_arr.copy()
        fv = v.ndim > 2 or (v.ndim == 2 and v.shape[0] > k)

        if fv:
            mk_mask = np.eye(p) == 0
            fulliv = np.any(v[np.tile(mk_mask[:, :, np.newaxis], (1, 1, k))] != 0)
            if not fulliv:
                diag_mask = np.eye(p) == 1
                v_diag = np.zeros((k, p))
                for ik in range(k):
                    v_diag[ik, :] = np.diag(v[:, :, ik]) / sx0 ** 2
                v = v_diag
            else:
                for ik in range(k):
                    v[:, :, ik] = v[:, :, ik] / np.outer(sx0, sx0)

    if len(wx_local) != n:
        raise ValueError(f'{n} datapoints but {len(wx_local)} weights')

    lsx = np.sum(np.log(sx0))
    xsw = xs * wx_local[:, np.newaxis]

    if not fulliv:
        # Diagonal covariance matrices EM
        v = np.maximum(v, c)
        xs2 = xs ** 2 * wx_local[:, np.newaxis]

        th = (l - np.floor(l)) * n
        sd = 1 if True else 0  # always compute final values
        lp_iter = int(np.floor(l)) + sd

        lpx = np.zeros(n)
        g = 0.0
        gg = np.zeros(lp_iter + 1)
        ss = sd

        for j_iter in range(lp_iter):
            g1 = g
            m1 = m.copy()
            v1 = v.copy()
            w1 = w.copy()

            vi = -0.5 / v
            lvm = np.log(w) - 0.5 * np.sum(np.log(v), axis=1)

            # Compute responsibilities
            py = np.zeros((k, n))
            for ik in range(k):
                diff = xs - m[ik, :]
                py[ik, :] = np.sum(diff ** 2 * vi[ik, :], axis=1) + lvm[ik]

            mx = np.max(py, axis=0)
            px = np.exp(py - mx[np.newaxis, :])
            ps = np.sum(px, axis=0)
            px = px / ps[np.newaxis, :]
            lpx = np.log(ps) + mx

            pk = px @ wx_local
            sx = px @ xsw
            sx2 = px @ xs2

            g = np.dot(lpx, wx_local)
            gg[j_iter] = g

            w = pk.copy()
            if np.all(pk > 0):
                m = sx / pk[:, np.newaxis]
                v_raw = sx2 / pk[:, np.newaxis]
            else:
                wm = pk == 0
                nz = np.sum(wm)
                mk_sort = np.argsort(lpx)
                m = np.zeros((k, p))
                v_raw = np.zeros((k, p))
                m[wm, :] = xs[mk_sort[:nz], :]
                w[wm] = 1.0 / n
                w = w * n / (n + nz)
                not_wm = ~wm
                m[not_wm, :] = sx[not_wm, :] / pk[not_wm, np.newaxis]
                v_raw = np.zeros((k, p))
                v_raw[not_wm, :] = sx2[not_wm, :] / pk[not_wm, np.newaxis]

            v = np.maximum(v_raw - m ** 2, c)

            if g - g1 <= th and j_iter > 0:
                if ss <= 0:
                    break
                ss -= 1

        # Calculate final probabilities
        pp = lpx - 0.5 * p * np.log(2 * np.pi) - lsx
        gg_out = gg[:j_iter + 1] - 0.5 * p * np.log(2 * np.pi) - lsx
        g = gg_out[-1]
        m = m1
        v = v1
        w = w1
        mm = np.sum(m, axis=0) / k
        f = (m.ravel() @ m.ravel() - k * (mm @ mm)) / np.sum(v)

        if not fv:
            m = m * sx0[np.newaxis, :] + mx0[np.newaxis, :]
            v = v * (sx0 ** 2)[np.newaxis, :]
        else:
            v_diag = v.copy()
            v = np.zeros((p, p, k))
            for ik in range(k):
                v[:, :, ik] = np.diag(v_diag[ik, :])
    else:
        # Full covariance EM - simplified for the common case
        # This path is taken when v0 contains full covariance matrices
        pp = np.zeros(n)
        f = 0.0
        gg_out = np.array([0.0])

        m = m * sx0[np.newaxis, :] + mx0[np.newaxis, :]
        if v.ndim == 2:
            v = v * (sx0 ** 2)[np.newaxis, :]
        else:
            for ik in range(k):
                v[:, :, ik] = v[:, :, ik] * np.outer(sx0, sx0)

    if fv and not fulliv:
        # Convert diagonal to full if 'v' was requested
        v_diag = v.copy()
        v = np.zeros((p, p, k))
        for ik in range(k):
            if v_diag.ndim == 2:
                v[:, :, ik] = np.diag(v_diag[ik, :])

    return m, v, w, g, f, pp, gg_out

v_gaussmixb

V_GAUSSMIXB - Approximate Bhattacharyya divergence between two GMMs.

v_gaussmixb 

v_gaussmixb(
    mf, vf=None, wf=None, mg=None, vg=None, wg=None, nx=1000
) -> tuple[float, ndarray]

Approximate Bhattacharyya divergence between two GMMs.

Parameters:

Name Type Description Default
mf array_like

Mixture means for GMM f, shape (kf, p).

 required
vf array_like

Variances for GMM f.

 None
wf array_like

Weights for GMM f.

 None
mg array_like

Mixture means for GMM g. If omitted, g=f.

 None
vg array_like

Variances for GMM g.

 None
wg array_like

Weights for GMM g.

 None
nx int

Number of samples for importance sampling (default 1000). 0 for upper bound only.

 1000

Returns:

Name Type Description
d float

Approximate Bhattacharyya divergence.
dbfg ndarray

Exact Bhattacharyya divergence between components.
Source code in pyvoicebox/v_gaussmixb.py 
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def v_gaussmixb(mf, vf=None, wf=None, mg=None, vg=None, wg=None, nx=1000) -> tuple[float, np.ndarray]:
    """Approximate Bhattacharyya divergence between two GMMs.

    Parameters
    ----------
    mf : array_like
        Mixture means for GMM f, shape (kf, p).
    vf : array_like, optional
        Variances for GMM f.
    wf : array_like, optional
        Weights for GMM f.
    mg : array_like, optional
        Mixture means for GMM g. If omitted, g=f.
    vg : array_like, optional
        Variances for GMM g.
    wg : array_like, optional
        Weights for GMM g.
    nx : int, optional
        Number of samples for importance sampling (default 1000). 0 for upper bound only.

    Returns
    -------
    d : float
        Approximate Bhattacharyya divergence.
    dbfg : ndarray
        Exact Bhattacharyya divergence between components.
    """
    mf = np.asarray(mf, dtype=float)
    if mf.ndim == 1:
        mf = mf.reshape(1, -1)
    kf, p = mf.shape

    if vf is None:
        vf = np.ones((kf, p))
    else:
        vf = np.asarray(vf, dtype=float)
    if wf is None:
        wf = np.full(kf, 1.0 / kf)
    else:
        wf = np.asarray(wf, dtype=float).ravel()

    if p == 1:
        vf = vf.ravel()
        dvf = True
    else:
        dvf = vf.ndim == 2 and vf.shape[0] == kf

    hpl2 = 0.5 * p * np.log(2)

    if mg is None:
        # Self-divergence
        dbfg = np.zeros((kf, kf))
        if kf > 1:
            if dvf:
                if vf.ndim == 2:
                    qldf = 0.25 * np.sum(np.log(vf), axis=1)
                else:
                    # 1D case: vf is (kf,)
                    qldf = 0.25 * np.log(vf)
                for jf in range(kf - 1):
                    for jg in range(jf + 1, kf):
                        if vf.ndim == 2:
                            vfpg = vf[jf, :] + vf[jg, :]
                            mdif = mf[jf, :] - mf[jg, :]
                            dbfg[jf, jg] = 0.25 * np.sum(mdif ** 2 / vfpg) + 0.5 * np.sum(np.log(vfpg)) - qldf[jf] - qldf[jg] - hpl2
                        else:
                            vfpg = vf[jf] + vf[jg]
                            mdif = mf[jf, 0] - mf[jg, 0]
                            dbfg[jf, jg] = 0.25 * mdif ** 2 / vfpg + 0.5 * np.log(vfpg) - qldf[jf] - qldf[jg] - hpl2
                        dbfg[jg, jf] = dbfg[jf, jg]
            else:
                qldf = np.zeros(kf)
                for jf in range(kf):
                    qldf[jf] = 0.5 * np.sum(np.log(np.diag(np.linalg.cholesky(vf[:, :, jf]))))
                for jf in range(kf - 1):
                    for jg in range(jf + 1, kf):
                        vfg = vf[:, :, jf] + vf[:, :, jg]
                        mdif = mf[jf, :] - mf[jg, :]
                        L = np.linalg.cholesky(vfg)
                        dbfg[jf, jg] = 0.25 * mdif @ np.linalg.solve(vfg, mdif) + np.sum(np.log(np.diag(L))) - qldf[jg] - qldf[jf] - hpl2
                        dbfg[jg, jf] = dbfg[jf, jg]
        d = 0.0
        return d, dbfg

    # Both f and g specified
    mg = np.asarray(mg, dtype=float)
    if mg.ndim == 1:
        mg = mg.reshape(1, -1)
    kg = mg.shape[0]

    if vg is None:
        vg = np.ones((kg, p))
    else:
        vg = np.asarray(vg, dtype=float)
    if wg is None:
        wg = np.full(kg, 1.0 / kg)
    else:
        wg = np.asarray(wg, dtype=float).ravel()

    if p == 1:
        vg = vg.ravel()
        dvg = True
    else:
        dvg = vg.ndim == 2 and vg.shape[0] == kg

    # Calculate pairwise Bhattacharyya divergences
    dbfg = np.zeros((kf, kg))
    if dvf and dvg:
        for jf in range(kf):
            for jg in range(kg):
                if vf.ndim == 2:
                    vfpg = vf[jf, :] + vg[jg, :]
                    mdif = mf[jf, :] - mg[jg, :]
                    qldf_val = 0.25 * np.sum(np.log(vf[jf, :]))
                    qldg_val = 0.25 * np.sum(np.log(vg[jg, :]))
                    dbfg[jf, jg] = 0.25 * np.sum(mdif ** 2 / vfpg) + 0.5 * np.sum(np.log(vfpg)) - qldf_val - qldg_val
                else:
                    vfpg = vf[jf] + vg[jg]
                    mdif = mf[jf, :] - mg[jg, :]
                    qldf_val = 0.25 * np.log(vf[jf])
                    qldg_val = 0.25 * np.log(vg[jg])
                    dbfg[jf, jg] = 0.25 * mdif[0] ** 2 / vfpg + 0.5 * np.log(vfpg) - qldf_val - qldg_val
    else:
        for jf in range(kf):
            if dvf:
                vjf = np.diag(vf[jf, :]) if vf.ndim == 2 else np.array([[vf[jf]]])
            else:
                vjf = vf[:, :, jf]
            qldf_val = 0.5 * np.sum(np.log(np.diag(np.linalg.cholesky(vjf))))
            for jg in range(kg):
                if dvg:
                    vjg = np.diag(vg[jg, :]) if vg.ndim == 2 else np.array([[vg[jg]]])
                else:
                    vjg = vg[:, :, jg]
                vfg = vjf + vjg
                mdif = mf[jf, :] - mg[jg, :]
                L = np.linalg.cholesky(vfg)
                qldg_val = 0.5 * np.sum(np.log(np.diag(np.linalg.cholesky(vjg))))
                dbfg[jf, jg] = 0.25 * mdif @ np.linalg.solve(vfg, mdif) + np.sum(np.log(np.diag(L))) - qldg_val - qldf_val

    dbfg -= hpl2

    # Variational bound iteration
    maxiter = 15
    lwf = np.tile(np.log(wf)[:, np.newaxis], (1, kg))
    lwg = np.tile(np.log(wg)[np.newaxis, :], (kf, 1))
    lhf = np.full((kf, kg), np.log(1.0 / kf))
    lhg = np.full((kf, kg), np.log(1.0 / kg))
    dbfg2 = 2 * dbfg
    dbfg2f = lwf - dbfg2
    dbfg2g = lwg - dbfg2
    dbfg2fg = dbfg2.ravel() - lwf.ravel() - lwg.ravel()

    dub = np.inf
    for ip in range(maxiter):
        dubp = dub
        dub = -v_logsum(0.5 * (lhf.ravel() + lhg.ravel() - dbfg2fg))

        if dub >= dubp:
            break

        lhg = lhf + dbfg2f
        lhg = lhg - v_logsum(lhg, 1)[:, np.newaxis]

        dub = -v_logsum(0.5 * (lhf.ravel() + lhg.ravel() - dbfg2fg))

        lhf = lhg + dbfg2g
        lhf = lhf - v_logsum(lhf, 0)[np.newaxis, :]

    if nx == 0:
        d = float(dub)
    else:
        d = float(dub)  # simplified: use upper bound

    return d, dbfg

v_gaussmixd

V_GAUSSMIXD - Marginal and conditional Gaussian mixture densities.

v_gaussmixd 

v_gaussmixd(
    y,
    m,
    v,
    w,
    a=None,
    b=None,
    f=None,
    g=None,
    return_mixtures=False,
) -> tuple[ndarray, ndarray]

Compute conditional GMM densities.

Parameters:

Name Type Description Default
y array_like

Conditioning data, shape (n, q).

 required
m array_like

Mixture means, shape (k, p).

 required
v array_like

Variances: diagonal (k, p) or full (p, p, k).

 required
w array_like

Mixture weights (k,).

 required
a array_like

Conditioning transformation matrix.

 None
b array_like

Dimension selection indices (1-based, MATLAB convention).

 None
f array_like

Output transformation matrix.

 None
g array_like

Output offset or dimension selection.

 None
return_mixtures bool

If True, return per-mixture means and weights.

 False

Returns:

Name Type Description
mz ndarray

Global mean of z for each y (n, r), or per-mixture if return_mixtures=True.
vz ndarray

Global or per-mixture covariances.
wz ndarray

Mixture weights (only if return_mixtures=True).
Source code in pyvoicebox/v_gaussmixd.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def v_gaussmixd(y, m, v, w, a=None, b=None, f=None, g=None, return_mixtures=False) -> tuple[np.ndarray, np.ndarray]:
    """Compute conditional GMM densities.

    Parameters
    ----------
    y : array_like
        Conditioning data, shape (n, q).
    m : array_like
        Mixture means, shape (k, p).
    v : array_like
        Variances: diagonal (k, p) or full (p, p, k).
    w : array_like
        Mixture weights (k,).
    a : array_like, optional
        Conditioning transformation matrix.
    b : array_like, optional
        Dimension selection indices (1-based, MATLAB convention).
    f : array_like, optional
        Output transformation matrix.
    g : array_like, optional
        Output offset or dimension selection.
    return_mixtures : bool, optional
        If True, return per-mixture means and weights.

    Returns
    -------
    mz : ndarray
        Global mean of z for each y (n, r), or per-mixture if return_mixtures=True.
    vz : ndarray
        Global or per-mixture covariances.
    wz : ndarray
        Mixture weights (only if return_mixtures=True).
    """
    y = np.asarray(y, dtype=float)
    if y.ndim == 1:
        y = y.reshape(1, -1)
    n, q = y.shape

    m = np.asarray(m, dtype=float)
    if m.ndim == 1:
        m = m.reshape(1, -1)
    k, p = m.shape

    v = np.asarray(v, dtype=float)
    w = np.asarray(w, dtype=float).ravel()

    fv = v.ndim > 2 or (v.ndim == 2 and v.shape[0] > k)

    # Set up conditioning transformation
    anull = a is None
    if anull:
        a_mat = np.eye(p)
        if b is None:
            b_idx = np.arange(q)
        else:
            b_idx = np.asarray(b, dtype=int).ravel() - 1  # Convert to 0-based
        a_mat = a_mat[b_idx, :]
        b_saved = b_idx.copy()
        b_vec = np.zeros(q)
    else:
        a_mat = np.asarray(a, dtype=float)
        if b is None:
            b_vec = np.zeros(q)
        else:
            b_vec = np.asarray(b, dtype=float).ravel()

    # Set up output transformation
    if f is None:
        f_mat = np.eye(p)
        if g is None:
            if anull:
                # Complement of selected dimensions
                all_dims = set(range(p))
                f_mat = np.eye(p)
                remaining = sorted(all_dims - set(b_saved))
                f_mat = f_mat[remaining, :]
            else:
                f_mat = f_mat[q:, :]
        else:
            g_idx = np.asarray(g, dtype=int).ravel() - 1  # Convert to 0-based
            f_mat = f_mat[g_idx, :]
        r = f_mat.shape[0]
        g_vec = np.zeros(r)
    else:
        f_mat = np.asarray(f, dtype=float)
        r = f_mat.shape[0]
        if g is None:
            g_vec = np.zeros(r)
        else:
            g_vec = np.asarray(g, dtype=float).ravel()

    yb = y - b_vec[np.newaxis, :]

    # Find mixture weights given y
    lp, wz_all = v_gaussmixp(yb, m, v, w, a_mat)[:2]

    mz = np.zeros((n, r, k))

    for i in range(k):
        if fv:
            vi = v[:, :, i]
        else:
            vi = np.diag(v[i, :])

        avat = a_mat @ vi @ a_mat.T
        hi = vi @ a_mat.T @ np.linalg.solve(avat, np.eye(avat.shape[0]))
        vzi = f_mat @ (vi - hi @ a_mat @ vi) @ f_mat.T
        mi = m[i, :]
        m0 = (mi - mi @ a_mat.T @ hi.T) @ f_mat.T + g_vec
        mz[:, :, i] = m0[np.newaxis, :] + yb @ hi.T @ f_mat.T

    if not return_mixtures:
        # Compute global mean
        mt = np.zeros((n, r))
        for i in range(k):
            mt += mz[:, :, i] * wz_all[:, i:i + 1]

        # Compute global variance
        vz_global = np.zeros((r, r, n))
        for idx_n in range(n):
            for i in range(k):
                if fv:
                    vi_full = v[:, :, i]
                else:
                    vi_full = np.diag(v[i, :])
                avat = a_mat @ vi_full @ a_mat.T
                hi = vi_full @ a_mat.T @ np.linalg.solve(avat, np.eye(avat.shape[0]))
                vzi = f_mat @ (vi_full - hi @ a_mat @ vi_full) @ f_mat.T
                dm = mz[idx_n, :, i] - mt[idx_n, :]
                vz_global[:, :, idx_n] += wz_all[idx_n, i] * (vzi + np.outer(dm, dm))

        return mt, vz_global
    else:
        # Return per-mixture results
        vz_mix = np.zeros((r, r, k))
        for i in range(k):
            if fv:
                vi = v[:, :, i]
            else:
                vi = np.diag(v[i, :])
            avat = a_mat @ vi @ a_mat.T
            hi = vi @ a_mat.T @ np.linalg.solve(avat, np.eye(avat.shape[0]))
            vz_mix[:, :, i] = f_mat @ (vi - hi @ a_mat @ vi) @ f_mat.T

        mz_out = np.transpose(mz, (2, 1, 0))  # (k, r, n)
        return mz_out, vz_mix, wz_all

v_gaussmixg

V_GAUSSMIXG - Global mean, variance and mode of a GMM (computation only).

v_gaussmixg 

v_gaussmixg(
    m, v=None, w=None, n_modes=1
) -> tuple[ndarray, ndarray, ndarray, ndarray]

Compute global mean, variance and modes of a GMM.

Parameters:

Name Type Description Default
m array_like

Mixture means, shape (k, p).

 required
v array_like

Variances: diagonal (k, p) or full (p, p, k). Default: ones.

 None
w array_like

Mixture weights. Default: uniform.

 None
n_modes int

Maximum number of modes to find (default 1).

 1

Returns:

Name Type Description
mg ndarray

Global mean, shape (p,).
vg ndarray

Global covariance, shape (p, p).
pg ndarray

Sorted list of modes, shape (n_modes, p).
pv ndarray

Log PDF at the modes, shape (n_modes,).
Source code in pyvoicebox/v_gaussmixg.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def v_gaussmixg(m, v=None, w=None, n_modes=1) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Compute global mean, variance and modes of a GMM.

    Parameters
    ----------
    m : array_like
        Mixture means, shape (k, p).
    v : array_like, optional
        Variances: diagonal (k, p) or full (p, p, k). Default: ones.
    w : array_like, optional
        Mixture weights. Default: uniform.
    n_modes : int, optional
        Maximum number of modes to find (default 1).

    Returns
    -------
    mg : ndarray
        Global mean, shape (p,).
    vg : ndarray
        Global covariance, shape (p, p).
    pg : ndarray
        Sorted list of modes, shape (n_modes, p).
    pv : ndarray
        Log PDF at the modes, shape (n_modes,).
    """
    m = np.asarray(m, dtype=float)
    if m.ndim == 1:
        m = m.reshape(1, -1)
    k, p = m.shape

    if v is None:
        v = np.ones((k, p))
    else:
        v = np.asarray(v, dtype=float)

    if w is None:
        w = np.ones(k)
    else:
        w = np.asarray(w, dtype=float).ravel()

    full = v.ndim > 2 or (k == 1 and v.size > p)
    if full and p == 1:
        v = v.reshape(k, 1)
        full = False

    w = w / np.sum(w)

    # Global mean
    mg = w @ m

    # Global covariance
    mz = m - mg[np.newaxis, :]
    if full:
        vg = mz.T @ (mz * w[:, np.newaxis]) + np.reshape(
            np.reshape(v, (p * p, k)) @ w, (p, p)
        )
    else:
        vg = mz.T @ (mz * w[:, np.newaxis]) + np.diag(w @ v)

    # Mode finding using fixed-point iteration
    nfp = 2
    maxloop = 60
    ssf = 0.1

    pg = m.copy()  # initialize mode candidates to mixture means
    nx = k

    if not full:
        vi = -0.5 / v
        lvm = np.log(w) - 0.5 * np.sum(np.log(v), axis=1)
        vim = vi * m
        ss_sep = np.sqrt(np.min(v, axis=0)) * ssf / np.sqrt(p)
    else:
        vi_stack = np.zeros((p * k, p))
        vim_stack = np.zeros(p * k)
        mtk_stack = np.zeros(p * k)
        lvm = np.zeros(k)
        for i in range(k):
            eigvals, eigvecs = np.linalg.eigh(v[:, :, i])
            if np.any(eigvals <= 0):
                raise ValueError(f'Covariance matrix for mixture {i} is not positive definite')
            vik = -0.5 * eigvecs @ np.diag(1.0 / eigvals) @ eigvecs.T
            sl = slice(i * p, (i + 1) * p)
            vi_stack[sl, :] = vik
            vim_stack[sl] = vik @ m[i, :]
            mtk_stack[sl] = m[i, :]
            lvm[i] = np.log(w[i]) - 0.5 * np.sum(np.log(eigvals))
        ss_sep = np.sqrt(np.min([v[j, j, :].min() for j in range(p)])) * ssf / np.sqrt(p)
        ss_sep = np.full(p, ss_sep)

    sv = 0.01 * ss_sep

    for it in range(maxloop):
        pg0 = pg.copy()

        if not full:
            # Compute probabilities
            py = np.zeros((k, nx))
            for ik in range(k):
                diff = pg - m[ik, :]
                py[ik, :] = np.sum(diff ** 2 * vi[ik, :], axis=1) + lvm[ik]

            mx_val = np.max(py, axis=0)
            px = np.exp(py - mx_val[np.newaxis, :])
            ps = np.sum(px, axis=0)
            px = (px / ps[np.newaxis, :]).T  # (nx, k)

            # Fixed point update
            pxvim = px @ vim  # (nx, p)
            pxvi = px @ vi  # (nx, p)
            pgf = pxvim / pxvi
        else:
            # Full covariance version
            py = np.zeros((k, nx))
            for ik in range(k):
                sl = slice(ik * p, (ik + 1) * p)
                vik = vi_stack[sl, :]
                vim_k = vim_stack[sl]
                mk = m[ik, :]
                for j_pt in range(nx):
                    yi = pg[j_pt, :]
                    py[ik, j_pt] = (vik @ yi - vim_k) @ (yi - mk) + lvm[ik]

            mx_val = np.max(py, axis=0)
            px = np.exp(py - mx_val[np.newaxis, :])
            ps = np.sum(px, axis=0)
            px = (px / ps[np.newaxis, :]).T  # (nx, k)

            # Fixed point update
            vif = np.zeros((k, p * p))
            vimf = np.zeros((k, p))
            for ik in range(k):
                sl = slice(ik * p, (ik + 1) * p)
                vif[ik, :] = vi_stack[sl, :].ravel()
                vimf[ik, :] = vim_stack[sl]

            pxvif = px @ vif  # (nx, p*p)
            pxvimf = px @ vimf  # (nx, p)
            pgf = np.zeros((nx, p))
            for j_pt in range(nx):
                pgf[j_pt, :] = np.linalg.solve(pxvif[j_pt, :].reshape(p, p), pxvimf[j_pt, :])

        if it <= nfp:
            pg = pgf
        else:
            pg = pgf  # simplified: just use fixed point

        if np.all(np.abs(pg - pg0) < sv[np.newaxis, :]):
            break

        # Remove duplicate modes
        if nx > 1:
            keep = np.ones(nx, dtype=bool)
            for i_mode in range(nx):
                if not keep[i_mode]:
                    continue
                for j_mode in range(i_mode + 1, nx):
                    if not keep[j_mode]:
                        continue
                    if np.all(np.abs(pg[i_mode, :] - pg[j_mode, :]) < ss_sep):
                        keep[j_mode] = False
            if not np.all(keep):
                pg = pg[keep, :]
                nx = pg.shape[0]

    # Calculate log pdf at each mode
    pv_arr = np.zeros(nx)
    for j_pt in range(nx):
        lp_val, _, _, _ = v_gaussmixp(pg[j_pt:j_pt + 1, :], m, v, w)
        pv_arr[j_pt] = lp_val[0]

    # Sort by decreasing probability
    ix = np.argsort(-pv_arr)
    pv_arr = pv_arr[ix]
    pg = pg[ix, :]

    if n_modes < len(pv_arr):
        pg = pg[:n_modes, :]
        pv_arr = pv_arr[:n_modes]

    return mg, vg, pg, pv_arr

v_gaussmixk

V_GAUSSMIXK - Approximate KL divergence between two GMMs.

v_gaussmixk 

v_gaussmixk(
    mf, vf=None, wf=None, mg=None, vg=None, wg=None
) -> tuple[float, ndarray]

Approximate Kullback-Leibler divergence between two GMMs.

Parameters:

Name Type Description Default
mf array_like

Mixture means for GMM f, shape (kf, p).

 required
vf array_like

Variances for GMM f. Diagonal (kf, p) or full (p, p, kf).

 None
wf array_like

Weights for GMM f.

 None
mg array_like

Mixture means for GMM g.

 None
vg array_like

Variances for GMM g.

 None
wg array_like

Weights for GMM g.

 None

Returns:

Name Type Description
d float

Approximate KL divergence D(f||g).
klfg ndarray

Exact KL divergence between components.
Source code in pyvoicebox/v_gaussmixk.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def v_gaussmixk(mf, vf=None, wf=None, mg=None, vg=None, wg=None) -> tuple[float, np.ndarray]:
    """Approximate Kullback-Leibler divergence between two GMMs.

    Parameters
    ----------
    mf : array_like
        Mixture means for GMM f, shape (kf, p).
    vf : array_like, optional
        Variances for GMM f. Diagonal (kf, p) or full (p, p, kf).
    wf : array_like, optional
        Weights for GMM f.
    mg : array_like, optional
        Mixture means for GMM g.
    vg : array_like, optional
        Variances for GMM g.
    wg : array_like, optional
        Weights for GMM g.

    Returns
    -------
    d : float
        Approximate KL divergence D(f||g).
    klfg : ndarray
        Exact KL divergence between components.
    """
    mf = np.asarray(mf, dtype=float)
    if mf.ndim == 1:
        mf = mf.reshape(1, -1)
    kf, p = mf.shape

    if wf is None or len(np.asarray(wf).ravel()) == 0:
        wf = np.full(kf, 1.0 / kf)
    else:
        wf = np.asarray(wf, dtype=float).ravel()

    if vf is None or len(np.asarray(vf).ravel()) == 0:
        vf = np.ones((kf, p))
    else:
        vf = np.asarray(vf, dtype=float)

    fvf = vf.ndim > 2 or (vf.ndim == 2 and vf.shape[0] > kf)

    # Calculate intra-f KL divergences
    klff = np.zeros((kf, kf))
    ixdp = np.arange(0, p * p, p + 1)  # diagonal indices

    if fvf:
        dvf = np.zeros(kf)
        for i in range(kf):
            dvf[i] = np.log(np.linalg.det(vf[:, :, i]))
        for j in range(kf):
            pfj = np.linalg.inv(vf[:, :, j])
            for i in range(kf):
                if i == j:
                    continue
                mffj = mf[i, :] - mf[j, :]
                pfjvf = pfj @ vf[:, :, i]
                klff[i, j] = 0.5 * (dvf[j] - p - dvf[i] + np.trace(pfjvf) + mffj @ pfj @ mffj)
    else:
        dvf = np.log(np.prod(vf, axis=1))
        pf = 1.0 / vf
        for j in range(kf):
            for i in range(kf):
                if i == j:
                    continue
                mffj = mf[i, :] - mf[j, :]
                klff[i, j] = 0.5 * (dvf[j] - dvf[i] + np.sum(vf[i, :] * pf[j, :]) - p + np.sum(mffj ** 2 * pf[j, :]))

    if mg is None:
        d = 0.0
        klfg = klff
    else:
        mg = np.asarray(mg, dtype=float)
        if mg.ndim == 1:
            mg = mg.reshape(1, -1)
        kg = mg.shape[0]

        if vg is None or len(np.asarray(vg).ravel()) == 0:
            vg = np.ones((kg, p))
        else:
            vg = np.asarray(vg, dtype=float)

        if wg is None or len(np.asarray(wg).ravel()) == 0:
            wg = np.full(kg, 1.0 / kg)
        else:
            wg = np.asarray(wg, dtype=float).ravel()

        fvg = vg.ndim > 2 or (vg.ndim == 2 and vg.shape[0] > kg)

        klfg = np.zeros((kf, kg))

        if fvg:
            dvg = np.zeros(kg)
            for j in range(kg):
                dvg[j] = np.log(np.linalg.det(vg[:, :, j]))
            if fvf:
                for j in range(kg):
                    pgj = np.linalg.inv(vg[:, :, j])
                    for i in range(kf):
                        mfgj = mf[i, :] - mg[j, :]
                        pgjvf = pgj @ vf[:, :, i]
                        klfg[i, j] = 0.5 * (dvg[j] - p - dvf[i] + np.trace(pgjvf) + mfgj @ pgj @ mfgj)
            else:
                for j in range(kg):
                    pgj = np.linalg.inv(vg[:, :, j])
                    for i in range(kf):
                        mfgj = mf[i, :] - mg[j, :]
                        klfg[i, j] = 0.5 * (dvg[j] - p - dvf[i] + np.sum(vf[i, :] * np.diag(pgj)) + mfgj @ pgj @ mfgj)
        else:
            dvg = np.log(np.prod(vg, axis=1))
            pg = 1.0 / vg
            if fvf:
                for j in range(kg):
                    for i in range(kf):
                        mfgj = mf[i, :] - mg[j, :]
                        klfg[i, j] = 0.5 * (dvg[j] - dvf[i] + np.sum(np.diag(vf[:, :, i]) * pg[j, :]) - p + np.sum(mfgj ** 2 * pg[j, :]))
            else:
                for j in range(kg):
                    for i in range(kf):
                        mfgj = mf[i, :] - mg[j, :]
                        klfg[i, j] = 0.5 * (dvg[j] - dvf[i] + np.sum(vf[i, :] * pg[j, :]) - p + np.sum(mfgj ** 2 * pg[j, :]))

        # Calculate variational approximation
        d = wf @ (v_logsum(-klff, 1, wf) - v_logsum(-klfg, 1, wg))

    return d, klfg

v_gaussmixm

V_GAUSSMIXM - Estimate mean and variance of the magnitude of a GMM.

v_gaussmixm 

v_gaussmixm(
    m, v=None, w=None, z=None
) -> tuple[ndarray, ndarray]

Estimate mean and variance of the magnitude of a GMM.

Parameters:

Name Type Description Default
m array_like

Mixture means, shape (k, p).

 required
v array_like

Variances: diagonal (k, p) or full (p, p, k).

 None
w array_like

Mixture weights (k,).

 None
z array_like

Origins, shape (t, p).

 None

Returns:

Name Type Description
mm ndarray

Mean of |x-z|, shape (t,).
mc ndarray

Variance (or covariance) of |x-z|.
Source code in pyvoicebox/v_gaussmixm.py 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def v_gaussmixm(m, v=None, w=None, z=None) -> tuple[np.ndarray, np.ndarray]:
    """Estimate mean and variance of the magnitude of a GMM.

    Parameters
    ----------
    m : array_like
        Mixture means, shape (k, p).
    v : array_like, optional
        Variances: diagonal (k, p) or full (p, p, k).
    w : array_like, optional
        Mixture weights (k,).
    z : array_like, optional
        Origins, shape (t, p).

    Returns
    -------
    mm : ndarray
        Mean of |x-z|, shape (t,).
    mc : ndarray
        Variance (or covariance) of |x-z|.
    """
    m = np.asarray(m, dtype=float)
    if m.ndim == 1:
        m = m.reshape(1, -1)
    k, p = m.shape

    if z is None:
        z = np.zeros((1, p))
    else:
        z = np.asarray(z, dtype=float)
        if z.ndim == 1:
            z = z.reshape(1, -1)

    if w is None:
        w = np.ones(k)
    else:
        w = np.asarray(w, dtype=float).ravel()

    if v is None:
        v = np.ones((k, p))
    else:
        v = np.asarray(v, dtype=float)

    t = z.shape[0]
    w = w / np.sum(w)

    if p == 1:
        # Exact 1D formula
        s = np.sqrt(v.ravel())
        mt = m[:, np.newaxis] - z.T  # (k, t)
        mts = mt / s[:, np.newaxis]
        ncdf = norm.cdf(-mts)
        npdf = norm.pdf(-mts)
        mm = ((mts * (1 - 2 * ncdf) + 2 * npdf).T @ (s * w)).ravel()

        mc = np.diag(((mts ** 2 + 1).T @ (v.ravel() * w)).ravel())
        if t > 1:
            for it in range(t):
                for jt in range(it):
                    mc[it, jt] = w @ (
                        (v.ravel() + mt[:, it] * mt[:, jt]) * (1 - 2 * np.abs(ncdf[:, it] - ncdf[:, jt]))
                        + 2 * s * np.sign(mt[:, jt] - mt[:, it]) * (npdf[:, it] * mt[:, jt] - npdf[:, jt] * mt[:, it])
                    )
                    mc[jt, it] = mc[it, jt]
        mc = mc - np.outer(mm, mm)
    else:
        fullv = v.ndim > 2 or (v.ndim == 2 and v.shape[0] > k)

        if fullv:
            diag_idx = np.arange(p)
            ms = np.zeros((k, t))
            for i in range(k):
                ms[i, :] = np.sum(m[i, :] ** 2) + np.sum(np.diag(v[:, :, i])) - 2 * m[i, :] @ z.T + np.sum(z ** 2, axis=1)
            vs = np.zeros((k, t))
            for i in range(k):
                si = v[:, :, i]
                vsc_i = 2 * np.trace(si @ si)
                for jt in range(t):
                    zmi = m[i, :] - z[jt, :]
                    vs[i, jt] = vsc_i + 4 * zmi @ si @ zmi
        else:
            ms = (np.sum(m ** 2, axis=1) + np.sum(v, axis=1))[:, np.newaxis] - 2 * m @ z.T + np.sum(z ** 2, axis=1)[np.newaxis, :]
            vsc = np.sum(v * (2 * v + 4 * m ** 2), axis=1)
            vmz = 4 * (v * m) @ z.T
            vs = vsc[:, np.newaxis] - 2 * vmz + 4 * v @ (z ** 2).T

        nm = ms ** 2 / vs
        mmk = np.exp(gammaln(nm + 0.5) - gammaln(nm)) * np.sqrt(ms / nm)
        mm = (mmk.T @ w).ravel()
        mc = (ms.T @ w - mm ** 2).ravel()

    return mm, mc

v_gaussmixp

V_GAUSSMIXP - Calculate log probability densities from a Gaussian mixture model.

v_gaussmixp 

v_gaussmixp(
    y, m, v=None, w=None, a=None, b=None
) -> tuple[ndarray, ndarray, ndarray, ndarray]

Calculate log probability densities from a Gaussian mixture model.

Parameters:

Name Type Description Default
y array_like

Input data, shape (n, q).

 required
m array_like

Mixture means, shape (k, p).

 required
v array_like

Variances: diagonal (k, p) or full (p, p, k). Default: identity.

 None
w array_like

Mixture weights (k,), must sum to 1. Default: uniform.

 None
a array_like

Transformation matrix (q, p). y = x*A' + B'.

 None
b array_like

Offset vector (q,).

 None

Returns:

Name Type Description
lp ndarray

Log probability of each data point (n,).
rp ndarray

Relative probability of each mixture (n, k).
kh ndarray

Index of highest probability mixture (n,).
kp ndarray

Relative probability of highest mixture (n,).
Source code in pyvoicebox/v_gaussmixp.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def v_gaussmixp(y, m, v=None, w=None, a=None, b=None) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Calculate log probability densities from a Gaussian mixture model.

    Parameters
    ----------
    y : array_like
        Input data, shape (n, q).
    m : array_like
        Mixture means, shape (k, p).
    v : array_like, optional
        Variances: diagonal (k, p) or full (p, p, k). Default: identity.
    w : array_like, optional
        Mixture weights (k,), must sum to 1. Default: uniform.
    a : array_like, optional
        Transformation matrix (q, p). y = x*A' + B'.
    b : array_like, optional
        Offset vector (q,).

    Returns
    -------
    lp : ndarray
        Log probability of each data point (n,).
    rp : ndarray
        Relative probability of each mixture (n, k).
    kh : ndarray
        Index of highest probability mixture (n,).
    kp : ndarray
        Relative probability of highest mixture (n,).
    """
    m = np.asarray(m, dtype=float)
    y = np.asarray(y, dtype=float)
    if y.ndim == 1:
        y = y.reshape(-1, 1)
    if m.ndim == 1:
        m = m.reshape(1, -1)

    k, p = m.shape
    n, q = y.shape

    if w is None:
        w = np.full(k, 1.0 / k)
    else:
        w = np.asarray(w, dtype=float).ravel()

    if v is None:
        v = np.ones((k, p))
    else:
        v = np.asarray(v, dtype=float)

    # Determine if full covariance
    fv = v.ndim > 2 or (v.ndim == 2 and v.shape[0] > k)

    # Apply transformations
    if a is not None:
        a = np.asarray(a, dtype=float)
        if b is not None:
            b = np.asarray(b, dtype=float).ravel()
            m = m @ a.T + b[np.newaxis, :]
        else:
            m = m @ a.T
        if fv:
            v_new = np.zeros((q, q, k))
            for ik in range(k):
                v_new[:, :, ik] = a @ v[:, :, ik] @ a.T
            v = v_new
        else:
            # Check if a preserves diagonality
            if np.all(np.sum(a != 0, axis=1) == 1):
                v_new = np.zeros((k, q))
                for ik in range(k):
                    v_new[ik, :] = v[ik, :] @ (a ** 2).T
                v = v_new
            else:
                v_new = np.zeros((q, q, k))
                for ik in range(k):
                    v_new[:, :, ik] = (a * v[ik, :][np.newaxis, :]) @ a.T
                v = v_new
                fv = True
        q_dim = q
    elif q < p or (a is None and b is not None):
        if b is None:
            b_idx = np.arange(q)
        else:
            b_idx = np.asarray(b, dtype=int).ravel()
            if np.issubdtype(type(b_idx[0]), np.integer) or np.all(b_idx == b_idx.astype(int)):
                b_idx = b_idx.astype(int) - 1 if np.min(b_idx) >= 1 else b_idx.astype(int)
        m = m[:, b_idx]
        if fv:
            v = v[np.ix_(b_idx, b_idx, np.arange(k))]
        else:
            v = v[:, b_idx]
        q_dim = q
    else:
        q_dim = q

    fv = v.ndim > 2 or (v.ndim == 2 and v.shape[0] > k)

    lp = np.zeros(n)
    rp = np.zeros((n, k))

    if n > 0:
        if not fv:
            # Diagonal covariance matrices
            vi = -0.5 / v  # (k, q)
            lvm = np.log(w) - 0.5 * np.sum(np.log(v), axis=1)  # (k,)

            # Compute log probabilities for all data points and all mixtures
            # py(k, n) = sum((y - m_k)^2 * vi_k) + lvm_k
            py = np.zeros((k, n))
            for ik in range(k):
                diff = y - m[ik, :]  # (n, q)
                py[ik, :] = np.sum(diff ** 2 * vi[ik, :], axis=1) + lvm[ik]

            mx = np.max(py, axis=0)  # (n,)
            px = np.exp(py - mx[np.newaxis, :])  # (k, n)
            ps = np.sum(px, axis=0)  # (n,)
            rp = (px / ps[np.newaxis, :]).T  # (n, k)
            lp = np.log(ps) + mx
        else:
            # Full covariance matrices
            vi = np.zeros((q_dim * k, q_dim))
            vim = np.zeros(q_dim * k)
            mtk = np.zeros(q_dim * k)
            lvm = np.zeros(k)

            for ik in range(k):
                vk = v[:, :, ik]
                eigvals, eigvecs = np.linalg.eigh(vk)
                if np.any(eigvals <= 0):
                    raise ValueError(f'Covariance matrix for mixture {ik} is not positive definite')
                vik = -0.5 * eigvecs @ np.diag(1.0 / eigvals) @ eigvecs.T
                sl = slice(ik * q_dim, (ik + 1) * q_dim)
                vi[sl, :] = vik
                vim[sl] = vik @ m[ik, :]
                mtk[sl] = m[ik, :]
                lvm[ik] = np.log(w[ik]) - 0.5 * np.sum(np.log(eigvals))

            # Compute for all points
            py = np.zeros((k, n))
            for ik in range(k):
                sl = slice(ik * q_dim, (ik + 1) * q_dim)
                vik = vi[sl, :]
                vim_k = vim[sl]
                mk = m[ik, :]
                for i in range(n):
                    yi = y[i, :]
                    py[ik, i] = (vik @ yi - vim_k) @ (yi - mk) + lvm[ik]

            mx = np.max(py, axis=0)
            px = np.exp(py - mx[np.newaxis, :])
            ps = np.sum(px, axis=0)
            rp = (px / ps[np.newaxis, :]).T
            lp = np.log(ps) + mx

        lp = lp - 0.5 * q_dim * np.log(2 * np.pi)

    kh = np.argmax(rp, axis=1)
    kp = rp[np.arange(n), kh]

    return lp, rp, kh, kp

v_gaussmixt

V_GAUSSMIXT - Multiply two GMM PDFs.

v_gaussmixt 

v_gaussmixt(
    m1, v1, w1, m2, v2, w2
) -> tuple[ndarray, ndarray, ndarray]

Multiply two GMM PDFs.

Parameters:

Name Type Description Default
m1 array_like

Mixture means, shape (k1, p) and (k2, p).

 required
m2 array_like

Mixture means, shape (k1, p) and (k2, p).

 required
v1 array_like

Variances: diagonal (ki, p) or full (p, p, ki).

 required
v2 array_like

Variances: diagonal (ki, p) or full (p, p, ki).

 required
w1 array_like

Mixture weights.

 required
w2 array_like

Mixture weights.

 required

Returns:

Name Type Description
m ndarray

Product mixture means, shape (k1*k2, p).
v ndarray

Product variances.
w ndarray

Product weights, shape (k1*k2,).
Source code in pyvoicebox/v_gaussmixt.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def v_gaussmixt(m1, v1, w1, m2, v2, w2) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Multiply two GMM PDFs.

    Parameters
    ----------
    m1, m2 : array_like
        Mixture means, shape (k1, p) and (k2, p).
    v1, v2 : array_like
        Variances: diagonal (ki, p) or full (p, p, ki).
    w1, w2 : array_like
        Mixture weights.

    Returns
    -------
    m : ndarray
        Product mixture means, shape (k1*k2, p).
    v : ndarray
        Product variances.
    w : ndarray
        Product weights, shape (k1*k2,).
    """
    m1 = np.asarray(m1, dtype=float)
    m2 = np.asarray(m2, dtype=float)
    v1 = np.asarray(v1, dtype=float)
    v2 = np.asarray(v2, dtype=float)
    w1 = np.asarray(w1, dtype=float).ravel()
    w2 = np.asarray(w2, dtype=float).ravel()

    if m1.ndim == 1:
        m1 = m1.reshape(1, -1)
    if m2.ndim == 1:
        m2 = m2.reshape(1, -1)

    k1, p = m1.shape
    k2 = m2.shape[0]
    f1 = v1.ndim > 2 or (v1.ndim == 2 and v1.shape[0] > k1)
    f2 = v2.ndim > 2 or (v2.ndim == 2 and v2.shape[0] > k2)

    k = k1 * k2
    j1 = np.tile(np.arange(k1), k2)
    j2 = np.repeat(np.arange(k2), k1)

    if p == 1:
        p1 = 1.0 / v1.ravel()
        p2 = 1.0 / v2.ravel()
        v = 1.0 / (p1[j1] + p2[j2])
        s1 = p1 * m1.ravel()
        s2 = p2 * m2.ravel()
        m = v * (s1[j1] + s2[j2])
        v12 = v1.ravel()[j1] + v2.ravel()[j2]
        wx = -0.5 * (m1.ravel()[j1] - m2.ravel()[j2]) ** 2 / v12
        wx = wx - np.max(wx)
        w = w1[j1] * w2[j2] * np.exp(wx) / np.sqrt(v12)
        w = w / np.sum(w)
        m = m.reshape(-1, 1)
        v = v.reshape(-1, 1)
    elif not f1 and not f2:
        # Both diagonal
        p1 = 1.0 / v1
        p2 = 1.0 / v2
        v = 1.0 / (p1[j1, :] + p2[j2, :])
        s1 = p1 * m1
        s2 = p2 * m2
        m = v * (s1[j1, :] + s2[j2, :])
        v12 = v1[j1, :] + v2[j2, :]
        wx = -0.5 * np.sum((m1[j1, :] - m2[j2, :]) ** 2 / v12, axis=1)
        wx = wx - np.max(wx)
        w = w1[j1] * w2[j2] * np.exp(wx) / np.sqrt(np.prod(v12, axis=1))
        w = w / np.sum(w)
    else:
        # At least one full covariance
        m = np.zeros((k, p))
        v = np.zeros((p, p, k))
        w = np.zeros(k)
        wx = np.zeros(k)

        for idx in range(k):
            i = j1[idx]
            j = j2[idx]

            if f1:
                v1i = v1[:, :, i]
                p1i = np.linalg.inv(v1i)
            else:
                v1i = np.diag(v1[i, :])
                p1i = np.diag(1.0 / v1[i, :])

            if f2:
                v2j = v2[:, :, j]
                p2j = np.linalg.inv(v2j)
            else:
                v2j = np.diag(v2[j, :])
                p2j = np.diag(1.0 / v2[j, :])

            pij = p1i + p2j
            vix = np.linalg.inv(pij)
            vij = v1i + v2j

            v[:, :, idx] = vix
            pm1 = m1[i, :] @ p1i
            pm2 = m2[j, :] @ p2j
            m[idx, :] = (pm1 + pm2) @ vix

            m12 = m1[i, :] - m2[j, :]
            wx[idx] = -0.5 * m12 @ np.linalg.solve(vij, m12)
            w[idx] = w1[i] * w2[j] / np.sqrt(np.linalg.det(vij))

        wx = wx - np.max(wx)
        w = w * np.exp(wx)
        w = w / np.sum(w)
        if k == 1:
            v = v[:, :, 0]

    return m, v, w

v_gmmlpdf

V_GMMLPDF - Obsolete wrapper for v_gaussmixp.

v_gmmlpdf 

v_gmmlpdf(*args, **kwargs) -> ndarray

Obsolete function - please use v_gaussmixp instead.

Source code in pyvoicebox/v_gmmlpdf.py 
7
8
9
def v_gmmlpdf(*args, **kwargs) -> np.ndarray:
    """Obsolete function - please use v_gaussmixp instead."""
    return v_gaussmixp(*args, **kwargs)

v_kmeans

V_KMEANS - K-means clustering algorithm.

v_kmeans 

v_kmeans(
    d, k, x0="f", l=300
) -> tuple[ndarray, ndarray, ndarray, ndarray]

Vector quantization using K-means algorithm.

Parameters:

Name Type Description Default
d array_like

Data vectors, shape (n, p).

 required
k int

Number of centres required.

 required
x0 str or array_like

Initial centres (k, p) or initialization method: 'f' = pick k random data points (default), 'p' = random partition centroids.

 'f'
l int

Maximum number of iterations (default 300). Use 0 to just compute distances for given centres.

 300

Returns:

Name Type Description
x ndarray

Output centres (k, p), or mean squared error if l=0.
g float or ndarray

Mean squared error, or cluster assignments if l=0.
j ndarray

Cluster assignments for each data point (0-based).
gg ndarray

Mean squared error at each iteration (only if l > 0).
Source code in pyvoicebox/v_kmeans.py 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def v_kmeans(d, k, x0='f', l=300) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Vector quantization using K-means algorithm.

    Parameters
    ----------
    d : array_like
        Data vectors, shape (n, p).
    k : int
        Number of centres required.
    x0 : str or array_like, optional
        Initial centres (k, p) or initialization method:
        'f' = pick k random data points (default),
        'p' = random partition centroids.
    l : int, optional
        Maximum number of iterations (default 300). Use 0 to just compute
        distances for given centres.

    Returns
    -------
    x : ndarray
        Output centres (k, p), or mean squared error if l=0.
    g : float or ndarray
        Mean squared error, or cluster assignments if l=0.
    j : ndarray
        Cluster assignments for each data point (0-based).
    gg : ndarray
        Mean squared error at each iteration (only if l > 0).
    """
    d = np.asarray(d, dtype=float)
    n, p = d.shape

    if isinstance(x0, str):
        if k < n:
            if 'p' in x0:
                # Random partition initialization
                ix = np.random.randint(0, k, size=n)
                forced = v_rnsubset(k, n)
                ix[forced] = np.arange(k)
                x = np.zeros((k, p))
                for i in range(k):
                    mask = ix == i
                    if np.any(mask):
                        x[i, :] = np.mean(d[mask, :], axis=0)
            else:
                # Forgy initialization: sample k centres
                x = d[v_rnsubset(k, n), :]
        else:
            x = d[np.arange(k) % n, :]
    else:
        x = np.asarray(x0, dtype=float).copy()

    m = np.zeros(n)
    j = np.zeros(n, dtype=int)
    gg = np.zeros(l)

    if l > 0:
        for ll in range(l):
            # Find closest centre
            z = v_disteusq(d, x, 'x')
            j = np.argmin(z, axis=1)
            m = z[np.arange(n), j]

            y = x.copy()

            # Calculate new centres
            nd = np.zeros(k)
            for i in range(k):
                nd[i] = np.sum(j == i)
            md = np.maximum(nd, 1)

            x_new = np.zeros((k, p))
            for i in range(k):
                mask = j == i
                if np.any(mask):
                    x_new[i, :] = np.sum(d[mask, :], axis=0) / md[i]
            x = x_new

            # Handle unused centres
            fx = np.where(nd == 0)[0]
            if len(fx) > 0:
                q = np.where(m != 0)[0]
                if len(q) <= len(fx):
                    x[fx[:len(q)], :] = d[q, :]
                else:
                    ri = np.random.permutation(len(q))
                    x[fx, :] = d[q[ri[:len(fx)]], :]

            gg[ll] = np.sum(m)
            if np.array_equal(x, y):
                gg = gg[:ll + 1] / n
                break
        else:
            gg = gg / n

        g = gg[-1]
        return x, g, j, gg
    else:
        # Just calculate distances
        z = v_disteusq(d, x, 'x')
        j = np.argmin(z, axis=1)
        m = z[np.arange(n), j]
        g_val = np.sum(m) / n
        return g_val, j, j, np.array([])

v_kmeanlbg

V_KMEANLBG - K-means using Linde-Buzo-Gray algorithm.

v_kmeanlbg 

v_kmeanlbg(d, k) -> tuple[ndarray, float, ndarray]

Vector quantization using the Linde-Buzo-Gray algorithm.

Parameters:

Name Type Description Default
d array_like

Data vectors (one per row), shape (n, p).

 required
k int

Number of centres required.

 required

Returns:

Name Type Description
x ndarray

Output centres, shape (k, p).
esq float

Mean squared error.
j ndarray

Cluster assignments (0-based).
Source code in pyvoicebox/v_kmeanlbg.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def v_kmeanlbg(d, k) -> tuple[np.ndarray, float, np.ndarray]:
    """Vector quantization using the Linde-Buzo-Gray algorithm.

    Parameters
    ----------
    d : array_like
        Data vectors (one per row), shape (n, p).
    k : int
        Number of centres required.

    Returns
    -------
    x : ndarray
        Output centres, shape (k, p).
    esq : float
        Mean squared error.
    j : ndarray
        Cluster assignments (0-based).
    """
    d = np.asarray(d, dtype=float)
    nc = d.shape[1]

    x, esq, j, _ = v_kmeans(d, 1)
    m_count = 1
    while m_count < k:
        n_split = min(m_count, k - m_count)
        m_count = m_count + n_split
        e = 1e-4 * np.sqrt(esq) * np.random.rand(1, nc)
        x0 = np.vstack([
            x[:n_split, :] + np.tile(e, (n_split, 1)),
            x[:n_split, :] - np.tile(e, (n_split, 1)),
            x[n_split:m_count - n_split, :]
        ])
        x, esq, j, _ = v_kmeans(d, m_count, x0)

    return x, esq, j

v_kmeanhar

V_KMEANHAR - K-harmonic means clustering algorithm.

v_kmeanhar 

v_kmeanhar(
    d, k, l=None, e=None, x0="f"
) -> tuple[ndarray, float, ndarray, ndarray]

Vector quantization using K-harmonic means algorithm.

Parameters:

Name Type Description Default
d array_like

Data vectors, shape (n, p).

 required
k int

Number of centres required.

 required
l float

Max iterations (integer part) + stopping threshold (fractional part). Default: 50.001.

 None
e int

Exponent in the cost function (default 4).

 None
x0 str or array_like

Initial centres or initialization method ('f' or 'p').

 'f'

Returns:

Name Type Description
x ndarray

Output centres, shape (k, p).
g float

Final performance criterion (normalized by n).
xn ndarray

Nearest centre for each input point (0-based).
gg ndarray

Performance criterion at each iteration.
Source code in pyvoicebox/v_kmeanhar.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def v_kmeanhar(d, k, l=None, e=None, x0='f') -> tuple[np.ndarray, float, np.ndarray, np.ndarray]:
    """Vector quantization using K-harmonic means algorithm.

    Parameters
    ----------
    d : array_like
        Data vectors, shape (n, p).
    k : int
        Number of centres required.
    l : float, optional
        Max iterations (integer part) + stopping threshold (fractional part).
        Default: 50.001.
    e : int, optional
        Exponent in the cost function (default 4).
    x0 : str or array_like, optional
        Initial centres or initialization method ('f' or 'p').

    Returns
    -------
    x : ndarray
        Output centres, shape (k, p).
    g : float
        Final performance criterion (normalized by n).
    xn : ndarray
        Nearest centre for each input point (0-based).
    gg : ndarray
        Performance criterion at each iteration.
    """
    d = np.asarray(d, dtype=float)
    n, p = d.shape

    if e is None:
        e = 4
    if l is None:
        l = 50 + 1e-3

    sd = 5  # number of times we must be below threshold

    # Initialize
    if isinstance(x0, str):
        if k < n:
            if 'p' in x0:
                ix = np.random.randint(0, k, size=n)
                forced = v_rnsubset(k, n)
                ix[forced] = np.arange(k)
                x = np.zeros((k, p))
                for i in range(k):
                    mask = ix == i
                    if np.any(mask):
                        x[i, :] = np.mean(d[mask, :], axis=0)
            else:
                x = d[v_rnsubset(k, n), :]
        else:
            x = d[np.arange(k) % n, :]
    else:
        x = np.asarray(x0, dtype=float).copy()

    eh = e / 2.0
    th = l - np.floor(l)
    max_iter = int(np.floor(l)) + 1  # extra loop to compute final value
    if max_iter <= 1:
        max_iter = 100
    if th == 0:
        th = -1

    gg = np.zeros(max_iter + 1)
    xn = np.zeros(n, dtype=int)

    ss = sd + 1
    g = 0

    for j in range(max_iter):
        g1 = g
        x1 = x.copy()

        # Compute squared distances (n x k)
        # py(k, n) in MATLAB, here we compute as (n, k) then transpose
        diff = d[:, np.newaxis, :] - x[np.newaxis, :, :]  # (n, k, p)
        py = np.sum(diff ** 2, axis=2).T  # (k, n)

        dm = np.min(py, axis=0)  # (n,)
        xn = np.argmin(py, axis=0)  # (n,)

        dmk = np.tile(dm, (k, 1))  # (k, n)
        dq = py > dmk
        pr = np.ones((k, n))
        pr[dq] = dmk[dq] / py[dq]

        pe = pr ** eh
        se = np.sum(pe, axis=0)  # (n,)

        xf = dm ** (eh - 1) / se  # (n,)
        g = np.dot(xf, dm)  # scalar

        xg = xf / se  # (n,)
        qik = xg[np.newaxis, :] * pe * pr  # (k, n)
        qk = np.sum(qik, axis=1)  # (k,)
        xs = qik @ d  # (k, p)

        gg[j] = g
        x = xs / qk[:, np.newaxis]

        if g1 - g <= th * g1:
            ss -= 1
            if ss <= 0:
                break
        else:
            ss = sd

    j_final = j + 1
    gg = gg[:j_final] * k / n
    g_final = gg[-1]

    # Go back to previous x values
    x = x1

    return x, g_final, xn, gg

Probability density functions

v_lognmpdf

V_LOGNMPDF - Calculate PDF of a multivariate lognormal distribution.

v_lognmpdf 

v_lognmpdf(x, m=None, v=None) -> ndarray

Calculate PDF of a multivariate lognormal distribution.

Parameters:

Name Type Description Default
x array_like

Points at which to evaluate, shape (n, d).

 required
m array_like

Mean vector (d,). Default: ones.

 None
v array_like

Covariance matrix (d, d) or diagonal vector (d,). Default: identity.

 None

Returns:

Name Type Description
p ndarray

PDF values, shape (n,).
Source code in pyvoicebox/v_lognmpdf.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def v_lognmpdf(x, m=None, v=None) -> np.ndarray:
    """Calculate PDF of a multivariate lognormal distribution.

    Parameters
    ----------
    x : array_like
        Points at which to evaluate, shape (n, d).
    m : array_like, optional
        Mean vector (d,). Default: ones.
    v : array_like, optional
        Covariance matrix (d, d) or diagonal vector (d,). Default: identity.

    Returns
    -------
    p : ndarray
        PDF values, shape (n,).
    """
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x.reshape(-1, 1)
    n, d = x.shape

    if m is None:
        m = np.ones(d)
    else:
        m = np.asarray(m, dtype=float).ravel()

    if v is None:
        v = np.eye(d)
    else:
        v = np.asarray(v, dtype=float)
        if v.ndim == 1:
            v = np.diag(v)

    # Convert from natural domain mean/covariance to log domain
    # Using the standard lognormal parameterization:
    # If X ~ LogN(mu, sigma), then E[X] = exp(mu + sigma^2/2)
    # For multivariate case: mu_log and sigma_log are the parameters
    # of the underlying normal distribution.
    #
    # Given natural mean m and covariance v:
    # sigma_log_ij = log(1 + v_ij / (m_i * m_j))
    # mu_log_i = log(m_i) - 0.5 * sigma_log_ii

    sigma_log = np.log(1 + v / np.outer(m, m))
    mu_log = np.log(m) - 0.5 * np.diag(sigma_log)

    p = np.zeros(n)
    c = np.prod(x, axis=1)
    q = c > 0
    if np.any(q):
        log_x = np.log(x[q, :])
        rv = multivariate_normal(mean=mu_log, cov=sigma_log)
        p[q] = rv.pdf(log_x) / c[q]

    return p

v_normcdflog

V_NORMCDFLOG - Log of normal CDF, accurate for large negative values.

v_normcdflog 

v_normcdflog(x, m=None, s=None) -> ndarray

Calculate log of Normal Cumulative Distribution function.

Parameters:

Name Type Description Default
x array_like

Input data.

 required
m float

Mean of Normal distribution (default 0).

 None
s float

Std deviation of Normal distribution (default 1).

 None

Returns:

Name Type Description
p ndarray

log(normcdf(x)); same shape as x.
Source code in pyvoicebox/v_normcdflog.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def v_normcdflog(x, m=None, s=None) -> np.ndarray:
    """Calculate log of Normal Cumulative Distribution function.

    Parameters
    ----------
    x : array_like
        Input data.
    m : float, optional
        Mean of Normal distribution (default 0).
    s : float, optional
        Std deviation of Normal distribution (default 1).

    Returns
    -------
    p : ndarray
        log(normcdf(x)); same shape as x.
    """
    x = np.asarray(x, dtype=float)
    if s is not None:
        x = (x - m) / s
    elif m is not None:
        x = x - m

    a = 0.996
    b = -22.2491306156561  # precalculated cutoff

    t = x < b
    p = np.zeros_like(x)
    p[~t] = np.real(np.log(0.5 * erfc(-x[~t] * np.sqrt(0.5))))
    p[t] = -0.5 * (x[t] ** 2 + np.log(2 * np.pi)) - np.real(
        np.log(-x[t] - a / x[t])
    )
    return p

v_chimv

V_CHIMV - Approximate mean and variance of non-central chi distribution.

v_chimv 

v_chimv(n, l=0, s=1) -> tuple[ndarray, ndarray]

Approximate mean and variance of non-central chi distribution.

Parameters:

Name Type Description Default
n int

Degrees of freedom.

 required
l float or array_like

Non-centrality parameter (default 0).

 0
s float

Standard deviation of Gaussian (default 1).

 1

Returns:

Name Type Description
m ndarray

Mean of chi distribution.
v ndarray

Variance of chi distribution.
Source code in pyvoicebox/v_chimv.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def v_chimv(n, l=0, s=1) -> tuple[np.ndarray, np.ndarray]:
    """Approximate mean and variance of non-central chi distribution.

    Parameters
    ----------
    n : int
        Degrees of freedom.
    l : float or array_like, optional
        Non-centrality parameter (default 0).
    s : float, optional
        Standard deviation of Gaussian (default 1).

    Returns
    -------
    m : ndarray
        Mean of chi distribution.
    v : ndarray
        Variance of chi distribution.
    """
    pp = np.array([0.595336298258636, -1.213013700592756, -0.018016200037799, 1.999986150447582, 0.0])
    qq = np.array([-0.161514114798972, 0.368983655790737, -0.136992134476950, -0.499681107630725, 2.0])

    l = np.asarray(l, dtype=float)
    ls = l / s
    l2 = ls ** 2
    s2 = s ** 2

    if n == 1:
        m = l * (1.0 - 2.0 * norm.cdf(-ls)) + 2.0 * s * norm.pdf(-ls)
    else:
        nab = 200
        ni = 1.0 / n
        ab_a = np.polyval(qq, ni)
        ab_b = np.polyval(pp, ni)
        m = np.sqrt(l2 + n - 1 + (ab_a + ab_b * l2) ** (-1)) * s

    v = (n + l2) * s2 - m ** 2
    return m, v

v_vonmisespdf

V_VONMISESPDF - Von Mises probability distribution.

v_vonmisespdf 

v_vonmisespdf(x, m, k) -> ndarray

Von Mises probability distribution.

Parameters:

Name Type Description Default
x array_like

Input values (in radians).

 required
m float

Mean angle (in radians).

 required
k float

Concentration parameter.

 required

Returns:

Name Type Description
p ndarray

Probability density values (same shape as x).
Source code in pyvoicebox/v_vonmisespdf.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def v_vonmisespdf(x, m, k) -> np.ndarray:
    """Von Mises probability distribution.

    Parameters
    ----------
    x : array_like
        Input values (in radians).
    m : float
        Mean angle (in radians).
    k : float
        Concentration parameter.

    Returns
    -------
    p : ndarray
        Probability density values (same shape as x).
    """
    x = np.asarray(x, dtype=float)
    p = np.exp(k * np.cos(x - m)) / (2 * np.pi * i0(k))
    return p

v_maxgauss

V_MAXGAUSS - Gaussian approximation to the max of a Gaussian vector.

v_maxgauss 

v_maxgauss(
    m, c=None, d=None
) -> tuple[float, float, ndarray, ndarray]

Determine Gaussian approximation to max of a Gaussian vector.

Parameters:

Name Type Description Default
m array_like

Mean vector of length N.

 required
c array_like

Covariance matrix (N,N) or vector of variances (N,).

 None
d array_like

Covariance w.r.t. some other variables (N,K).

 None

Returns:

Name Type Description
u float

Mean of max(x).
v float

Variance of max(x).
p ndarray

Probability of each element being the max (N,).
r ndarray

Covariance between max(x) and the d-variables (1,K).
Source code in pyvoicebox/v_maxgauss.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def v_maxgauss(m, c=None, d=None) -> tuple[float, float, np.ndarray, np.ndarray]:
    """Determine Gaussian approximation to max of a Gaussian vector.

    Parameters
    ----------
    m : array_like
        Mean vector of length N.
    c : array_like, optional
        Covariance matrix (N,N) or vector of variances (N,).
    d : array_like, optional
        Covariance w.r.t. some other variables (N,K).

    Returns
    -------
    u : float
        Mean of max(x).
    v : float
        Variance of max(x).
    p : ndarray
        Probability of each element being the max (N,).
    r : ndarray
        Covariance between max(x) and the d-variables (1,K).
    """
    nrh = -np.sqrt(0.5)
    kpd = np.sqrt(0.5 / np.pi)

    m = np.asarray(m, dtype=float).ravel()
    nm = len(m)

    if c is None:
        c = np.eye(nm)
    else:
        c = np.asarray(c, dtype=float)
        if c.ndim == 1 or (c.ndim == 2 and c.shape[1] == 1 and c.shape[0] == nm):
            c = np.diag(c.ravel())
        elif c.ndim == 0:
            c = np.diag(np.full(nm, float(c)))

    p = np.eye(nm)

    # Remove negative infinities
    ix = np.where(m == -np.inf)[0]
    if len(ix) > 0:
        m = np.delete(m, ix)
        c = np.delete(np.delete(c, ix, axis=0), ix, axis=1)
        p = np.delete(p, ix, axis=1)
        nm = len(m)

    while nm > 1:
        # Build strict upper triangle indices
        ix_list = []
        for col in range(nm):
            for row in range(col):
                ix_list.append((row, col))

        # Calculate scaled differences in means
        cd = np.diag(c)
        gm_arr = []
        gv_arr = []
        for row, col in ix_list:
            gm_arr.append(m[row] - m[col])
            gv_arr.append(cd[row] + cd[col] - c[row, col] - c[col, row])

        gm_arr = np.array(gm_arr)
        gv_arr = np.array(gv_arr)

        jx_zero = np.where(gv_arr <= 0)[0]
        if len(jx_zero) > 0:
            jx = jx_zero[0]
            i, j = ix_list[jx]
            dm = gm_arr[jx]
            if dm > 0:
                m = np.delete(m, j)
                c = np.delete(np.delete(c, j, axis=0), j, axis=1)
                p = np.delete(p, j, axis=1)
            else:
                m = np.delete(m, i)
                c = np.delete(np.delete(c, i, axis=0), i, axis=1)
                p = np.delete(p, i, axis=1)
        else:
            ratios = gm_arr ** 2 / gv_arr
            jx = np.argmax(ratios)
            i, j = ix_list[jx]

            dm = gm_arr[jx]
            ds = np.sqrt(gv_arr[jx])
            dms = dm / ds
            q = 0.5 * erfc(nrh * dms)
            f = kpd * np.exp(-0.5 * dms ** 2)

            mi_val = m[i]
            mj_val = m[j]
            u_val = dm * q + mj_val + ds * f
            v_val = (mi_val + mj_val - u_val) * u_val + cd[i] * q + cd[j] * (1 - q) - mi_val * mj_val

            m[i] = u_val
            c[i, :] = q * c[i, :] + (1 - q) * c[j, :]
            c[:, i] = c[i, :]
            c[i, i] = v_val
            p[:, i] = q * p[:, i] + (1 - q) * p[:, j]

            m = np.delete(m, j)
            c = np.delete(np.delete(c, j, axis=0), j, axis=1)
            p = np.delete(p, j, axis=1)

        nm = len(m)

    u = m[0]
    v = c[0, 0]
    p = p.ravel() / np.sum(p)

    r = None
    if d is not None:
        d = np.asarray(d, dtype=float)
        r = p @ d
        r = r.reshape(1, -1)

    return u, v, p, r

v_berk2prob

V_BERK2PROB - Convert Berksons (log-odds base 2) to probability.

v_berk2prob 

v_berk2prob(b) -> tuple[ndarray, ndarray]

Convert Berksons to probability.

Parameters:

Name Type Description Default
b array_like

Berkson values (log-odds base 2).

 required

Returns:

Name Type Description
p ndarray

Corresponding probability values.
d ndarray

Corresponding derivatives dP/dB.
Source code in pyvoicebox/v_berk2prob.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def v_berk2prob(b) -> tuple[np.ndarray, np.ndarray]:
    """Convert Berksons to probability.

    Parameters
    ----------
    b : array_like
        Berkson values (log-odds base 2).

    Returns
    -------
    p : ndarray
        Corresponding probability values.
    d : ndarray
        Corresponding derivatives dP/dB.
    """
    b = np.asarray(b, dtype=float)
    p = 1.0 - 1.0 / (1.0 + np.power(2.0, b))
    d = np.log(2.0) * p * (1.0 - p)
    return p, d

v_prob2berk

V_PROB2BERK - Convert probability to Berksons (log-odds base 2).

v_prob2berk 

v_prob2berk(p) -> tuple[ndarray, ndarray]

Convert probability to Berksons.

Parameters:

Name Type Description Default
p array_like

Probability values.

 required

Returns:

Name Type Description
b ndarray

Corresponding Berkson values.
d ndarray

Corresponding derivatives dP/dB.
Source code in pyvoicebox/v_prob2berk.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def v_prob2berk(p) -> tuple[np.ndarray, np.ndarray]:
    """Convert probability to Berksons.

    Parameters
    ----------
    p : array_like
        Probability values.

    Returns
    -------
    b : ndarray
        Corresponding Berkson values.
    d : ndarray
        Corresponding derivatives dP/dB.
    """
    p = np.asarray(p, dtype=float)
    b = np.log2(p / (1.0 - p))
    d = np.log(2.0) * p * (1.0 - p)
    return b, d

v_gausprod

V_GAUSPROD - Calculate the product of Gaussians.

v_gausprod 

v_gausprod(m, c=None) -> tuple[float, ndarray, ndarray]

Calculate the product of n d-dimensional multivariate Gaussians.

Parameters:

Name Type Description Default
m array_like

Means, shape (d, n) - each column is the mean of one Gaussian.

 required
c array_like

Covariance matrices. Can be: - (d, d, n) for full covariance - (d, n) for diagonal - (n,) for c*I - omitted for I

 None

Returns:

Name Type Description
g float

Log gain (= log(integral)).
u ndarray

Mean vector (d, 1).
k ndarray

Covariance matrix (same form as input).
Source code in pyvoicebox/v_gausprod.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def v_gausprod(m, c=None) -> tuple[float, np.ndarray, np.ndarray]:
    """Calculate the product of n d-dimensional multivariate Gaussians.

    Parameters
    ----------
    m : array_like
        Means, shape (d, n) - each column is the mean of one Gaussian.
    c : array_like, optional
        Covariance matrices. Can be:
        - (d, d, n) for full covariance
        - (d, n) for diagonal
        - (n,) for c*I
        - omitted for I

    Returns
    -------
    g : float
        Log gain (= log(integral)).
    u : ndarray
        Mean vector (d, 1).
    k : ndarray
        Covariance matrix (same form as input).
    """
    m = np.asarray(m, dtype=float)
    d, n = m.shape

    if c is None:
        c = np.ones(n)

    c = np.asarray(c, dtype=float)
    sc = c.shape

    if c.ndim < 3:
        if c.ndim == 1 and c.shape[0] == n and n != d:
            # Covariance matrices are multiples of the identity
            jj = 1
            jj2 = 2
            gj = np.zeros(n)
            while jj < n:
                for j_idx in range(jj, n, jj2):
                    k_idx = j_idx - jj
                    cjk = c[k_idx] + c[j_idx]
                    dm = m[:, k_idx] - m[:, j_idx]
                    gj[k_idx] = gj[k_idx] + gj[j_idx] - 0.5 * (d * np.log(cjk) + dm @ dm / cjk)
                    m[:, k_idx] = (c[k_idx] * m[:, j_idx] + c[j_idx] * m[:, k_idx]) / cjk
                    c[k_idx] = c[k_idx] * c[j_idx] / cjk
                jj = jj2
                jj2 = 2 * jj
            g = gj[0]
            k = c[0]
            u = m[:, 0]
        elif c.ndim == 2 and sc[0] == d and sc[1] == n:
            # Diagonal covariance matrices
            jj = 1
            jj2 = 2
            gj = np.zeros(n)
            while jj < n:
                for j_idx in range(jj, n, jj2):
                    k_idx = j_idx - jj
                    cjk = c[:, k_idx] + c[:, j_idx]
                    dm = m[:, k_idx] - m[:, j_idx]
                    ix = cjk > d * np.max(cjk) * np.finfo(float).eps
                    piv = np.zeros(d)
                    piv[ix] = 1.0 / cjk[ix]
                    gj[k_idx] = gj[k_idx] + gj[j_idx] - 0.5 * (np.log(np.prod(cjk)) + piv @ (dm ** 2))
                    m[:, k_idx] = piv * (c[:, k_idx] * m[:, j_idx] + c[:, j_idx] * m[:, k_idx])
                    c[:, k_idx] = c[:, k_idx] * piv * c[:, j_idx]
                jj = jj2
                jj2 = 2 * jj
            g = gj[0]
            k = c[:, 0]
            u = m[:, 0]
        else:
            # Handle c as scalar*I case when n==d (ambiguous) -- treat as multiples of I if 1D
            if c.ndim == 1 and c.shape[0] == n:
                jj = 1
                jj2 = 2
                gj = np.zeros(n)
                while jj < n:
                    for j_idx in range(jj, n, jj2):
                        k_idx = j_idx - jj
                        cjk = c[k_idx] + c[j_idx]
                        dm = m[:, k_idx] - m[:, j_idx]
                        gj[k_idx] = gj[k_idx] + gj[j_idx] - 0.5 * (d * np.log(cjk) + dm @ dm / cjk)
                        m[:, k_idx] = (c[k_idx] * m[:, j_idx] + c[j_idx] * m[:, k_idx]) / cjk
                        c[k_idx] = c[k_idx] * c[j_idx] / cjk
                    jj = jj2
                    jj2 = 2 * jj
                g = gj[0]
                k = c[0]
                u = m[:, 0]
            else:
                raise ValueError("Ambiguous covariance specification")
    else:
        # Full covariance matrices (d, d, n)
        jj = 1
        jj2 = 2
        gj = np.zeros(n)
        while jj < n:
            for j_idx in range(jj, n, jj2):
                k_idx = j_idx - jj
                cjk = c[:, :, k_idx] + c[:, :, j_idx]
                dm = m[:, k_idx] - m[:, j_idx]
                piv = np.linalg.pinv(cjk)
                gj[k_idx] = gj[k_idx] + gj[j_idx] - 0.5 * (np.log(np.linalg.det(cjk)) + dm @ piv @ dm)
                m[:, k_idx] = c[:, :, k_idx] @ piv @ m[:, j_idx] + c[:, :, j_idx] @ piv @ m[:, k_idx]
                c[:, :, k_idx] = c[:, :, k_idx] @ piv @ c[:, :, j_idx]
                c[:, :, k_idx] = 0.5 * (c[:, :, k_idx] + c[:, :, k_idx].T)
            jj = jj2
            jj2 = 2 * jj
        g = gj[0]
        k = c[:, :, 0]
        u = m[:, 0]

    g = g - 0.5 * (n - 1) * d * np.log(2 * np.pi)
    return g, u, k

v_histndim

V_HISTNDIM - Generate an n-dimensional histogram.

v_histndim 

v_histndim(x, b=None, mode='') -> tuple[ndarray, ndarray]

Generate an n-dimensional histogram.

Parameters:

Name Type Description Default
x array_like

Input data, shape (m, d).

 required
b array_like

Histogram bin specification, shape (3, d) or (1, d) or (3, 1). Row 0: number of bins; Row 1: min of first bin; Row 2: max of last bin. Default: 10 bins per dimension.

 None
mode str

'p' to scale as probabilities; 'z' for zero base in 2D plot.

 ''

Returns:

Name Type Description
v ndarray

Histogram counts (or probabilities if 'p' in mode).
t list of ndarray

Bin boundary values for each dimension.
Source code in pyvoicebox/v_histndim.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def v_histndim(x, b=None, mode='') -> tuple[np.ndarray, np.ndarray]:
    """Generate an n-dimensional histogram.

    Parameters
    ----------
    x : array_like
        Input data, shape (m, d).
    b : array_like, optional
        Histogram bin specification, shape (3, d) or (1, d) or (3, 1).
        Row 0: number of bins; Row 1: min of first bin; Row 2: max of last bin.
        Default: 10 bins per dimension.
    mode : str, optional
        'p' to scale as probabilities; 'z' for zero base in 2D plot.

    Returns
    -------
    v : ndarray
        Histogram counts (or probabilities if 'p' in mode).
    t : list of ndarray
        Bin boundary values for each dimension.
    """
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x.reshape(-1, 1)
    n, d = x.shape

    if b is None:
        b = np.full((1, d), 10, dtype=float)

    b = np.asarray(b, dtype=float)
    if b.ndim == 1:
        b = b.reshape(-1, 1)
    if b.shape[1] == 1:
        b = np.tile(b, (1, d))

    if b.shape[0] < 3:
        mi = np.min(x, axis=0)
        ma = np.max(x, axis=0)
        w = (ma - mi) / (b[0, :] - 0.001)
        if b.shape[0] < 3:
            b_full = np.zeros((3, d))
            b_full[0, :] = b[0, :]
            b_full[2, :] = ma + 0.0005 * w
            b_full[1, :] = mi - 0.0005 * w
            b = b_full

    acd = np.where(b[0, :] > 0)[0]
    sv = b[0, acd].astype(int)
    nbt = int(np.prod(sv))
    t = []

    ok = np.ones(n, dtype=bool)
    ix = np.full(n, nbt - int(np.sum(np.cumprod(sv))), dtype=int)
    k = 1
    for idx, j in enumerate(acd):
        nbins = int(b[0, j])
        bw = nbins / (b[2, j] - b[1, j])
        bi = np.ceil((x[:, j] - b[1, j]) * bw).astype(int)
        ok = ok & (bi > 0) & (bi <= nbins)
        ix[ok] = ix[ok] + k * bi[ok]
        k = k * nbins
        t.append(b[1, j] + np.arange(nbins + 1) / bw)

    # Build histogram using sparse-like accumulation
    v = np.zeros(nbt)
    valid_ix = ix[ok]
    # Ensure indices are valid (0-based)
    valid_ix = valid_ix - 1  # convert from 1-based to 0-based
    valid_ix = np.clip(valid_ix, 0, nbt - 1)
    np.add.at(v, valid_ix, 1)

    if len(sv) > 1:
        v = v.reshape(sv, order='F')

    if 'p' in mode:
        v = v / n

    return v, t

v_pdfmoments

V_PDFMOMENTS - Convert between central moments, raw moments and cumulants.

v_pdfmoments 

v_pdfmoments(
    t, m, b=0, a=1
) -> tuple[ndarray, ndarray, ndarray]

Convert between central moments, raw moments and cumulants.

Parameters:

Name Type Description Default
t str

Input/output type string. Lower case = input type, upper case = output type. 'm'/'M' = central moments, 'r'/'R' = raw moments, 'k'/'K' = cumulants.

 required
m array_like

Vector of input moments; m[0] is always the mean.

 required
b float

Output moments are for a*x + b (default 0).

 0
a float

Output moments are for a*x + b (default 1).

 1

Returns:

Name Type Description
c ndarray

Central moments (or as determined by 'R' or 'K' options).
r ndarray

Raw moments.
k ndarray

Cumulants.
Source code in pyvoicebox/v_pdfmoments.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
def v_pdfmoments(t, m, b=0, a=1) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Convert between central moments, raw moments and cumulants.

    Parameters
    ----------
    t : str
        Input/output type string. Lower case = input type, upper case = output type.
        'm'/'M' = central moments, 'r'/'R' = raw moments, 'k'/'K' = cumulants.
    m : array_like
        Vector of input moments; m[0] is always the mean.
    b : float, optional
        Output moments are for a*x + b (default 0).
    a : float, optional
        Output moments are for a*x + b (default 1).

    Returns
    -------
    c : ndarray
        Central moments (or as determined by 'R' or 'K' options).
    r : ndarray
        Raw moments.
    k : ndarray
        Cumulants.
    """
    m = np.asarray(m, dtype=float).ravel()
    n = len(m)

    # Precompute binomial coefficients
    bc = np.zeros((n, n + 1))
    for i in range(n):
        for j in range(i + 2):
            bc[i, j] = comb(i + 1, j, exact=True)

    # Precompute cumulant coefficients
    # This is a complex recursive algorithm; we implement it following the MATLAB code
    fa = np.ones(max(1, n // 2))
    for i in range(1, len(fa)):
        fa[i] = (i + 1) * fa[i - 1]

    mk = [None] * n
    # Initialize: mk[0] = [[1, 1, 1]] meaning powers=1, coefs k->m=1, m->k=1
    mk[0] = np.array([[1, 1, 1]])
    if n > 1:
        mk[1] = np.array([[0, 1, 1, 1]])  # for moment 3 (index 2)

    mn = np.zeros((n, n))
    mn[0, 0] = 1  # mn[0,0] = 1
    if n > 1:
        mn[1, 0] = 1
        mn[1, 1] = 1

    for i in range(3, n):
        j = (i + 1) // 2  # first coefficient row to sum
        nr = 1
        for rr in range(j - 1, i - 2):
            col_idx = i - rr - 3
            if rr < mn.shape[0] and col_idx >= 0 and col_idx < mn.shape[1]:
                nr += int(mn[rr, col_idx])

        mki = np.zeros((nr, i + 2))
        ix = 0
        mki[0, i - 2] = 1  # power of moment i-1
        mki[0, i] = 1  # coef k->m
        mki[0, i + 1] = 1  # coef m->k

        ix = 1
        for rr in range(j, i - 1):
            nk_idx = rr - 1
            col_idx = i - rr - 3
            if nk_idx >= len(mk) or mk[nk_idx] is None:
                continue
            if nk_idx < mn.shape[0] and col_idx >= 0 and col_idx < mn.shape[1]:
                nk = int(mn[nk_idx, col_idx])
            else:
                nk = 0
            if nk == 0:
                continue

            mkk = mk[nk_idx]
            mkik = mkk[:nk, :rr].copy()
            col_to_inc = i - rr - 2
            if col_to_inc < mkik.shape[1]:
                mkik[:, col_to_inc] += 1
            if ix + nk <= mki.shape[0]:
                mki[ix:ix + nk, :rr] = mkik
                # Calculate coefficient for k->m
                for qq in range(nk):
                    if col_to_inc < mkik.shape[1]:
                        mki[ix + qq, i] = mkk[qq, rr] * bc[i - 1, i - rr] / mkik[qq, col_to_inc]
                    rho = np.sum(mkik[qq, :]) - 1
                    rho_int = int(rho)
                    if rho_int > 0 and rho_int <= len(fa):
                        mki[ix + qq, i + 1] = mki[ix + qq, i] * fa[rho_int - 1] * (-1) ** rho_int
                    elif rho_int == 0:
                        mki[ix + qq, i + 1] = mki[ix + qq, i]
                ix += nk

        # Sort rows
        if ix > 0:
            mki = mki[:ix]
            idx = np.lexsort(mki[:, :i].T[::-1])
            mki = mki[idx]

        mk[i - 1] = mki
        # Update mn
        if i - 1 < mn.shape[0]:
            mn[i - 1, 0] = mki.shape[0]
            for cc in range(1, min(i - 1, mn.shape[1])):
                mn[i - 1, cc] = np.sum(np.all(mki[:, :cc] == 0, axis=1))

    # Apply scaling
    mu = a * m[0] + b
    c = m.copy()
    r = m.copy()
    k_out = m.copy()
    m_row = m.copy()

    # Determine input type
    if 'k' in t:
        tin = 3
        k_out = k_out * (a ** np.arange(1, n + 1))
        k_out[0] = 0
    elif 'r' in t:
        tin = 2
    else:
        tin = 1
        c = c * (a ** np.arange(1, n + 1))
        c[0] = 0

    tout = [
        'K' not in t and 'R' not in t,  # output c (central moments)
        True,  # always compute r
        True,  # always compute k
    ]

    for il in range(2):
        # Convert between moments
        if il == 0:
            v = np.concatenate(([1], m_row * (a ** np.arange(1, n + 1))))
            bb = b - mu
            doit = tin == 2 and (tout[0] or tout[2])
        else:
            if tin == 2:
                bb = b
            else:
                v = np.concatenate(([1], c))
                bb = mu
            doit = tout[1]

        if doit:
            y = v[1:].copy()
            if bb != 0:
                for i in range(n):
                    # polyval: bc[i, 0:i+2] * v[0:i+2] evaluated at bb
                    coeffs = bc[i, :i + 2] * v[:i + 2]
                    y[i] = np.polyval(coeffs[::-1], bb)
            if il == 0:
                c = y.copy()
            else:
                r = y.copy()

        # Convert cumulants to/from moments
        if il == 0:
            x = k_out.copy()
            doit = tin == 3 and (tout[0] or tout[1])
        else:
            x = c.copy()
            doit = tin < 3 and tout[2]

        if doit:
            y = x.copy()
            for i in range(3, n):
                if i - 1 < len(mk) and mk[i - 1] is not None:
                    mki = mk[i - 1]
                    col_idx = i + il  # il=0 -> column i (k->m), il=1 -> column i+1 (m->k)
                    if col_idx < mki.shape[1]:
                        terms = mki[:, col_idx] * np.prod(
                            np.power(
                                np.tile(x[1:i + 1], (mki.shape[0], 1)) if i + 1 <= len(x) else np.tile(np.concatenate([x[1:], np.zeros(i + 1 - len(x))]), (mki.shape[0], 1)),
                                mki[:, :i]
                            ),
                            axis=1
                        )
                        y[i] = np.sum(terms)
            if il == 0:
                c = y.copy()
            else:
                k_out = y.copy()

    c[0] = mu
    k_out[0] = mu

    if 'R' in t:
        c = r
    elif 'K' in t:
        c = k_out

    return c, r, k_out

v_besselratio

V_BESSELRATIO - Bessel function ratio I_{v+1}(x)/I_v(x).

v_besselratio 

v_besselratio(x, v=0, p=5) -> ndarray

Calculate the Bessel function ratio besseli(v+1,x)/besseli(v,x).

Parameters:

Name Type Description Default
x array_like

Bessel function argument.

 required
v int

Denominator Bessel function order (default 0).

 0
p int

Digits precision <=14 (default 5).

 5

Returns:

Name Type Description
y ndarray

Value of the Bessel function ratio.
Source code in pyvoicebox/v_besselratio.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def v_besselratio(x, v=0, p=5) -> np.ndarray:
    """Calculate the Bessel function ratio besseli(v+1,x)/besseli(v,x).

    Parameters
    ----------
    x : array_like
        Bessel function argument.
    v : int, optional
        Denominator Bessel function order (default 0).
    p : int, optional
        Digits precision <=14 (default 5).

    Returns
    -------
    y : ndarray
        Value of the Bessel function ratio.
    """
    wm = [1, 1, 2, 4, 10, 20, 14, 21, 16, 22, 18, 23, 20, 25]
    nm = [1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5]

    p = min(max(int(np.ceil(p)), 1), 14)
    n = nm[p - 1]  # number of iterations
    u = max(v, wm[p - 1])  # minimum value of v for first stage

    x = np.asarray(x, dtype=float)
    s = x.shape
    x = x.ravel()

    # Identify edge cases upfront
    edge_zero = x == 0
    edge_inf = x == np.inf
    finite = ~edge_zero & ~edge_inf

    # Work only on finite, nonzero values
    xf = x[finite]
    rf = np.zeros((xf.size, n + 1))
    for i in range(1, n + 2):
        rf[:, i - 1] = xf / (u + i - 0.5 + np.sqrt((u + i + 0.5) ** 2 + xf ** 2))

    for i in range(1, n + 1):
        for k in range(1, n - i + 2):
            rf[:, k - 1] = xf / (
                u + k + np.sqrt((u + k) ** 2 + xf ** 2 * rf[:, k] / rf[:, k - 1])
            )

    yf = rf[:, 0]
    for i in range(u, v, -1):
        yf = 1.0 / (yf + 2 * i / xf)

    # Assemble result
    y = np.empty_like(x)
    y[finite] = yf
    y[edge_zero] = 0.0
    y[edge_inf] = 1.0
    y = y.reshape(s)
    return y

v_besselratioi

V_BESSELRATIOI - Inverse Bessel function ratio.

v_besselratioi 

v_besselratioi(r, v=0, p=5) -> ndarray

Calculate the inverse Bessel function ratio.

Given r = besseli(v+1,s)/besseli(v,s), find s.

Parameters:

Name Type Description Default
r array_like

Value of the Bessel function ratio.

 required
v int

Denominator Bessel function order (default 0, must be 0 for now).

 0
p int

Digits precision <=13 (default 5).

 5

Returns:

Name Type Description
s ndarray

Value of s such that r = besseli(v+1,s)/besseli(v,s).
Source code in pyvoicebox/v_besselratioi.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def v_besselratioi(r, v=0, p=5) -> np.ndarray:
    """Calculate the inverse Bessel function ratio.

    Given r = besseli(v+1,s)/besseli(v,s), find s.

    Parameters
    ----------
    r : array_like
        Value of the Bessel function ratio.
    v : int, optional
        Denominator Bessel function order (default 0, must be 0 for now).
    p : int, optional
        Digits precision <=13 (default 5).

    Returns
    -------
    s : ndarray
        Value of s such that r = besseli(v+1,s)/besseli(v,s).
    """
    p = min(max(int(np.ceil(p)), 1), 13)
    if v != 0:
        raise ValueError('v must be zero (for now)')

    r = np.asarray(r, dtype=float)
    sr = r.shape
    r = r.ravel()
    nr = r.size

    s = np.full(nr, -1.0)
    y = 2.0 / (1.0 - r)

    m1 = (r >= 0) & (r < 1)  # valid range
    mn = np.copy(m1)  # Newton iteration needed

    # Use inverse taylor series for 0 <= r <= 0.85
    m = m1 & (r <= 0.85)
    if np.any(m):
        rm = r[m]
        mn[m] = rm >= 0.642
        xm = rm ** 2
        sm = ((rm - 5.6076) * rm + 5.0797) * rm - 4.6494
        sm = sm * y[m] * xm - 1.0
        s[m] = ((((sm * xm + 15.0) * xm + 60.0) * xm / 360.0 + 1.0) * xm - 2.0) * rm / (xm - 1.0)

    # Use continued fraction for 0.85 < r < 1
    m = m1 & ~(r <= 0.85)
    if np.any(m):
        rm = r[m].copy()
        mn[m] = rm < 0.95
        ym = y[m]
        mc = rm < 0.95
        if np.any(mc):
            rm[mc] = (-2326.0 * rm[mc] + 4317.5526) * rm[mc] - 2001.035224
        if np.any(~mc):
            rm[~mc] = 32.0 / (120.0 * rm[~mc] - 131.5 + ym[~mc])
        s[m] = (ym + 1.0 + 3.0 / (ym - 5.0 - 12.0 / (ym - 10.0 - rm))) * 0.25

    # Newton iterations
    if np.any(mn):
        rmn = r[mn]
        smn = s[mn]
        ymn = y[mn]
        ymn = ((0.00048 * ymn - 0.1589) * ymn + 0.744) * ymn - 4.2932
        smn = smn + (v_besselratio(smn, 0, p + 1) - rmn) * ymn
        mr = (rmn >= 0.75) & (rmn <= 0.875)
        if np.any(mr):
            smn[mr] = smn[mr] + (v_besselratio(smn[mr], 0, p + 1) - rmn[mr]) * ymn[mr]
        s[mn] = smn

    s = s.reshape(sr)
    return s

v_besratinv0

V_BESRATINV0 - Inverse of Modified Bessel Ratio I1(k)/I0(k).

v_besratinv0 

v_besratinv0(r) -> ndarray

Inverse function of the Modified Bessel Ratio I1(k)/I0(k).

Parameters:

Name Type Description Default
r array_like

Input argument in range [0,1].

 required

Returns:

Name Type Description
k ndarray

Output satisfying r = I1(k)/I0(k).
Source code in pyvoicebox/v_besratinv0.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def v_besratinv0(r) -> np.ndarray:
    """Inverse function of the Modified Bessel Ratio I1(k)/I0(k).

    Parameters
    ----------
    r : array_like
        Input argument in range [0,1].

    Returns
    -------
    k : ndarray
        Output satisfying r = I1(k)/I0(k).
    """
    r = np.asarray(r, dtype=float)
    k = np.zeros_like(r)

    m1 = r <= 0.85
    a = r[m1]
    if a.size > 0:
        y = 2.0 / (1.0 - a)
        x = a * a
        s = ((a - 5.6076) * a + 5.0797) * a - 4.6494
        s = s * y * x - 1.0
        s = ((((s * x + 15.0) * x + 60.0) * x / 360.0 + 1.0) * x - 2.0) * a / (x - 1.0)

        m2 = a >= 0.642
        b = a[m2]
        if b.size > 0:
            z = y[m2]
            yy = ((0.00048 * z - 0.1589) * z + 0.744) * z - 4.2932
            t = s[m2]
            t = (v_besselratio(t, 0, 9) - b) * yy + t
            m3 = b >= 0.75
            c = b[m3]
            if c.size > 0:
                u = t[m3]
                t[m3] = (v_besselratio(u, 0, 9) - c) * yy[m3] + u
            s[m2] = t
        k[m1] = s

    a = r[~m1]
    if a.size > 0:
        y = 2.0 / (1.0 - a)
        x = np.zeros_like(a)
        m2 = a > 0.95
        x[m2] = 32.0 / (120.0 * a[m2] - 131.5 + y[m2])
        x[~m2] = (-2326.0 * a[~m2] + 4317.5526) * a[~m2] - 2001.035224
        s = (y + 1.0 + 3.0 / (y - 5.0 - 12.0 / (y - 10.0 - x))) * 0.25

        m2 = a < 0.95
        b = a[m2]
        if b.size > 0:
            z = y[m2]
            yy = ((0.00048 * z - 0.1589) * z + 0.744) * z - 4.2932
            t = s[m2]
            t = (v_besselratio(t, 0, 9) - b) * yy + t
            m3 = b <= 0.875
            c = b[m3]
            if c.size > 0:
                u = t[m3]
                t[m3] = (v_besselratio(u, 0, 9) - c) * yy[m3] + u
            s[m2] = t
        k[~m1] = s

    return k