0
# Mathematical Functions
1

2
PyMC3 provides a comprehensive set of mathematical functions built on Theano for tensor operations, link functions, and specialized mathematical computations commonly used in Bayesian modeling. These functions support automatic differentiation and efficient computation.
3

4
## Capabilities
5

6
### Link Functions
7

8
Functions for transforming between constrained and unconstrained parameter spaces.
9

10
```python { .api }
11
def logit(p):
12
    """
13
    Logit (log-odds) transformation for probabilities.
14
    
15
    Transforms probability values to the real line using logit(p) = log(p/(1-p)).
16
    
17
    Parameters:
18
    - p: tensor, probability values in (0, 1)
19
    
20
    Returns:
21
    - tensor: logit-transformed values on real line
22
    """
23

24
def invlogit(x, eps=None):
25
    """
26
    Inverse logit (sigmoid) transformation.
27
    
28
    Transforms real values to probabilities using invlogit(x) = 1/(1+exp(-x)).
29
    
30
    Parameters:
31
    - x: tensor, real-valued input
32
    - eps: float, epsilon for numerical stability
33
    
34
    Returns:
35
    - tensor: probability values in (0, 1)
36
    """
37

38
def probit(p):
39
    """
40
    Probit transformation for probabilities.
41
    
42
    Transforms probability values using inverse normal CDF.
43
    
44
    Parameters:
45
    - p: tensor, probability values in (0, 1)
46
    
47
    Returns:
48
    - tensor: probit-transformed values
49
    """
50

51
def invprobit(x):
52
    """
53
    Inverse probit transformation.
54
    
55
    Transforms real values to probabilities using normal CDF.
56
    
57
    Parameters:
58
    - x: tensor, real-valued input
59
    
60
    Returns:
61
    - tensor: probability values in (0, 1)
62
    """
63
```
64

65
### Log-Space Operations
66

67
Numerically stable operations in log space for probability computations.
68

69
```python { .api }
70
def logsumexp(x, axis=None, keepdims=True):
71
    """
72
    Numerically stable log-sum-exp operation.
73
    
74
    Computes log(sum(exp(x))) in a numerically stable way by factoring
75
    out the maximum value before exponentiation.
76
    
77
    Parameters:
78
    - x: tensor, input values
79
    - axis: int or tuple, axis along which to sum
80
    - keepdims: bool, keep reduced dimensions
81
    
82
    Returns:
83
    - tensor: log-sum-exp result
84
    """
85

86
def logaddexp(a, b):
87
    """
88
    Numerically stable log(exp(a) + exp(b)).
89
    
90
    Parameters:
91
    - a: tensor, first input
92
    - b: tensor, second input
93
    
94
    Returns:
95
    - tensor: log(exp(a) + exp(b))
96
    """
97

98
def logdiffexp(a, b):
99
    """
100
    Numerically stable log(exp(a) - exp(b)) for a > b.
101
    
102
    Parameters:
103
    - a: tensor, larger input (a > b)
104
    - b: tensor, smaller input
105
    
106
    Returns:
107
    - tensor: log(exp(a) - exp(b))
108
    """
109

110
def log1pexp(x):
111
    """
112
    Numerically stable log(1 + exp(x)).
113
    
114
    Uses different computational strategies depending on input range
115
    to maintain numerical stability.
116
    
117
    Parameters:
118
    - x: tensor, input values
119
    
120
    Returns:
121
    - tensor: log(1 + exp(x))
122
    """
123

124
def log1mexp(x):
125
    """
126
    Numerically stable log(1 - exp(x)) for x < 0.
127
    
128
    Parameters:
129
    - x: tensor, input values (should be < 0)
130
    
131
    Returns:
132
    - tensor: log(1 - exp(x))
133
    """
134
```
135

136
### Special Functions
137

138
Specialized mathematical functions for Bayesian modeling.
139

140
```python { .api }
141
def logbern(log_p):
142
    """
143
    Log-Bernoulli sampling from log probabilities.
144
    
145
    Samples binary values from Bernoulli distribution given log probabilities,
146
    useful for categorical and discrete choice models.
147
    
148
    Parameters:
149
    - log_p: tensor, log probabilities
150
    
151
    Returns:
152
    - tensor: binary samples
153
    """
154

155
def expand_packed_triangular(n, packed, lower=True, diagonal_only=False):
156
    """
157
    Expand packed triangular matrix to full matrix.
158
    
159
    Converts packed lower or upper triangular representation to full
160
    symmetric matrix, commonly used for covariance matrices.
161
    
162
    Parameters:
163
    - n: int, matrix dimension
164
    - packed: tensor, packed triangular values
165
    - lower: bool, expand as lower triangular
166
    - diagonal_only: bool, expand only diagonal elements
167
    
168
    Returns:
169
    - tensor: full symmetric matrix
170
    """
171

172
def flatten_list(tensors):
173
    """
174
    Flatten list of tensors into single vector.
175
    
176
    Concatenates multiple tensors after flattening each to 1D,
177
    useful for parameter vectorization.
178
    
179
    Parameters:
180
    - tensors: list, tensors to flatten and concatenate
181
    
182
    Returns:
183
    - tensor: flattened concatenated vector
184
    """
185
```
186

187
### Matrix Operations
188

189
Specialized matrix operations for multivariate distributions and linear algebra.
190

191
```python { .api }
192
def kronecker(*Ks):
193
    """
194
    Kronecker product of multiple matrices.
195
    
196
    Computes kronecker product K1 ⊗ K2 ⊗ ... ⊗ Kn for efficient
197
    computation of multivariate and matrix-normal distributions.
198
    
199
    Parameters:
200
    - Ks: matrices, input matrices for kronecker product
201
    
202
    Returns:
203
    - tensor: kronecker product matrix
204
    """
205

206
def kron_diag(*diags):
207
    """
208
    Kronecker product of diagonal matrices.
209
    
210
    Efficient computation of kronecker product when all inputs
211
    are diagonal matrices.
212
    
213
    Parameters:
214
    - diags: tensors, diagonal elements of matrices
215
    
216
    Returns:
217
    - tensor: diagonal of resulting kronecker product
218
    """
219

220
def batched_diag(C):
221
    """
222
    Extract diagonal elements from batched matrices.
223
    
224
    Extracts diagonal elements from a batch of matrices efficiently.
225
    
226
    Parameters:
227
    - C: tensor, batch of matrices (..., n, n)
228
    
229
    Returns:
230
    - tensor: batch of diagonal elements (..., n)
231
    """
232

233
def block_diagonal(matrices, sparse=False, format='csr'):
234
    """
235
    Create block diagonal matrix from list of matrices.
236
    
237
    Constructs block diagonal matrix with input matrices on diagonal
238
    and zeros elsewhere.
239
    
240
    Parameters:
241
    - matrices: list, matrices to place on diagonal
242
    - sparse: bool, return sparse matrix
243
    - format: str, sparse matrix format if sparse=True
244
    
245
    Returns:
246
    - tensor: block diagonal matrix
247
    """
248

249
def cartesian(*arrays):
250
    """
251
    Cartesian product of input arrays.
252
    
253
    Generates cartesian product of input arrays, useful for
254
    creating parameter grids and design matrices.
255
    
256
    Parameters:
257
    - arrays: arrays to compute cartesian product of
258
    
259
    Returns:
260
    - tensor: cartesian product array
261
    """
262

263
def flat_outer(a, b):
264
    """
265
    Flattened outer product of two vectors.
266
    
267
    Computes outer product and flattens result to 1D vector.
268
    
269
    Parameters:
270
    - a: tensor, first input vector
271
    - b: tensor, second input vector
272
    
273
    Returns:
274
    - tensor: flattened outer product
275
    """
276
```
277

278
### Utility Functions
279

280
General utility functions for tensor operations and numerical computations.
281

282
```python { .api }
283
def tround(*args, **kwargs):
284
    """
285
    Theano-compatible rounding function.
286
    
287
    Rounds tensor values to nearest integer with proper gradient handling.
288
    
289
    Parameters:
290
    - args: positional arguments for rounding
291
    - kwargs: keyword arguments for rounding
292
    
293
    Returns:
294
    - tensor: rounded values
295
    """
296
```
297

298
## Usage Examples
299

300
### Link Function Transformations
301

302
```python
303
import pymc3 as pm
304
import numpy as np
305

306
with pm.Model() as model:
307
    # Unconstrained parameter
308
    alpha_raw = pm.Normal('alpha_raw', 0, 1)
309
    
310
    # Transform to probability using inverse logit
311
    alpha = pm.math.invlogit(alpha_raw)
312
    
313
    # Use in Bernoulli likelihood
314
    y = pm.Bernoulli('y', p=alpha, observed=data)
315
    
316
    trace = pm.sample(1000)
317

318
# Manual transformation
319
prob_values = np.array([0.1, 0.5, 0.9])
320
logit_values = pm.math.logit(prob_values)
321
back_to_prob = pm.math.invlogit(logit_values)
322
```
323

324
### Numerically Stable Log Operations
325

326
```python
327
import pymc3 as pm
328
import numpy as np
329

330
# Log-sum-exp for mixture models
331
with pm.Model() as model:
332
    # Mixture weights (in log space)
333
    log_weights = pm.Dirichlet('log_weights', a=[1, 1, 1]).log()
334
    
335
    # Component means
336
    mu = pm.Normal('mu', 0, 1, shape=3)
337
    
338
    # Log likelihoods for each component
339
    log_likes = pm.Normal.dist(mu, 1).logp(data[:, None])
340
    
341
    # Mixture log likelihood using logsumexp
342
    mixture_logp = pm.math.logsumexp(log_weights + log_likes, axis=1)
343
    
344
    # Custom likelihood
345
    pm.Potential('mixture', mixture_logp.sum())
346
    
347
    trace = pm.sample(1000)
348
```
349

350
### Matrix Operations for Multivariate Models
351

352
```python
353
import pymc3 as pm
354
import numpy as np
355

356
# Kronecker-structured covariance
357
with pm.Model() as model:
358
    # Row and column covariance matrices
359
    Sigma_row = pm.InverseWishart('Sigma_row', nu=5, V=np.eye(3))
360
    Sigma_col = pm.InverseWishart('Sigma_col', nu=4, V=np.eye(2))
361
    
362
    # Kronecker product covariance
363
    Sigma = pm.math.kronecker(Sigma_row, Sigma_col)
364
    
365
    # Matrix normal distribution
366
    mu = pm.Normal('mu', 0, 1, shape=(3, 2))
367
    X = pm.MvNormal('X', 
368
                    mu=mu.flatten(), 
369
                    cov=Sigma, 
370
                    observed=data.flatten())
371
    
372
    trace = pm.sample(1000)
373

374
# Block diagonal covariance
375
with pm.Model() as model:
376
    # Individual covariance blocks
377
    Sigma1 = pm.InverseWishart('Sigma1', nu=3, V=np.eye(2))
378
    Sigma2 = pm.InverseWishart('Sigma2', nu=4, V=np.eye(3))
379
    
380
    # Block diagonal structure
381
    Sigma = pm.math.block_diagonal([Sigma1, Sigma2])
382
    
383
    # Multivariate normal with block structure
384
    mu = pm.Normal('mu', 0, 1, shape=5)
385
    y = pm.MvNormal('y', mu=mu, cov=Sigma, observed=data)
386
    
387
    trace = pm.sample(1000)
388
```
389

390
### Packed Triangular Matrices
391

392
```python
393
import pymc3 as pm
394
import numpy as np
395

396
# Cholesky parameterization
397
with pm.Model() as model:
398
    # Packed lower triangular elements
399
    n_dim = 3
400
    n_packed = n_dim * (n_dim + 1) // 2
401
    
402
    packed_L = pm.Normal('packed_L', 0, 1, shape=n_packed)
403
    
404
    # Expand to full Cholesky factor
405
    L = pm.math.expand_packed_triangular(n_dim, packed_L, lower=True)
406
    
407
    # Covariance matrix
408
    Sigma = pm.math.dot(L, L.T)
409
    
410
    # Multivariate normal
411
    mu = pm.Normal('mu', 0, 1, shape=n_dim)
412
    y = pm.MvNormal('y', mu=mu, chol=L, observed=data)
413
    
414
    trace = pm.sample(1000)
415
```