0
# Probability and Combinatorics Functions
1

2
This document covers Math.js's probability functions, combinatorics operations, random number generation, and specialized mathematical functions for probability calculations and discrete mathematics.
3

4
## Import
5

6
```typescript
7
import {
8
  // Combinatorics
9
  combinations, factorial, gamma, lgamma, multinomial, permutations,
10
  // Special combinatorics  
11
  bellNumbers, catalan, composition, stirlingS2,
12
  // Random functions
13
  pickRandom, random, randomInt,
14
  // Distance/divergence
15
  kldivergence,
16
  // Set operations (related to combinatorics)
17
  setCartesian, setDifference, setDistinct, setIntersect, 
18
  setIsSubset, setMultiplicity, setPowerset, setSize,
19
  setSymDifference, setUnion
20
} from 'mathjs'
21
```
22

23
## Basic Combinatorics
24

25
### Combinations (Binomial Coefficients)
26

27
```typescript
28
combinations(n: MathType, k: MathType): MathType
29
```
30
{ .api }
31

32
```typescript
33
// "n choose k" - number of ways to choose k items from n items
34
combinations(5, 2) // 10 (ways to choose 2 from 5)
35
combinations(10, 3) // 120
36
combinations(4, 4) // 1 (only one way to choose all items)
37
combinations(4, 0) // 1 (only one way to choose none)
38

39
// Large numbers using BigNumber
40
combinations(bignumber('100'), bignumber('50')) // Very large result
41

42
// Properties
43
combinations(n, k) === combinations(n, subtract(n, k)) // Symmetry
44
combinations(n, 0) === 1 // Base case
45
combinations(n, n) === 1 // Base case
46

47
// Pascal's triangle relationship  
48
combinations(5, 2) === add(combinations(4, 1), combinations(4, 2)) // 10 = 4 + 6
49
```
50

51
### Permutations
52

53
```typescript
54
permutations(n: MathType, k?: MathType): MathType
55
```
56
{ .api }
57

58
```typescript
59
// Number of ways to arrange k items from n items
60
permutations(5, 2) // 20 (5 * 4)
61
permutations(5, 5) // 120 (5!)
62
permutations(10, 3) // 720 (10 * 9 * 8)
63

64
// Without second argument: n!
65
permutations(5) // 120 (same as factorial(5))
66
permutations(0) // 1
67

68
// Relationship to combinations
69
permutations(n, k) === multiply(combinations(n, k), factorial(k))
70

71
// Circular permutations
72
function circularPermutations(n) {
73
  return divide(factorial(n), n) // (n-1)!
74
}
75
circularPermutations(4) // 6 (ways to arrange 4 people around a table)
76
```
77

78
### Factorial
79

80
```typescript
81
factorial(n: MathType): MathType
82
```
83
{ .api }
84

85
```typescript
86
// Standard factorial: n! = n * (n-1) * ... * 2 * 1
87
factorial(0) // 1 (by definition)
88
factorial(1) // 1  
89
factorial(5) // 120
90
factorial(10) // 3628800
91

92
// Large factorials with BigNumber
93
factorial(bignumber('50')) // Very large result, maintains precision
94

95
// Double factorial (not built-in, custom implementation)
96
function doubleFactorial(n) {
97
  if (smaller(n, 0)) throw new Error('Negative input')
98
  if (smallerEq(n, 1)) return 1
99
  return multiply(n, doubleFactorial(subtract(n, 2)))
100
}
101
doubleFactorial(8) // 384 (8 * 6 * 4 * 2)
102
doubleFactorial(9) // 945 (9 * 7 * 5 * 3 * 1)
103

104
// Subfactorial (derangements)
105
function subfactorial(n) {
106
  if (equal(n, 0)) return 1
107
  if (equal(n, 1)) return 0
108
  return multiply(subtract(n, 1), add(subfactorial(subtract(n, 1)), subfactorial(subtract(n, 2))))
109
}
110
```
111

112
### Multinomial Coefficients
113

114
```typescript
115
multinomial(a: MathCollection): MathType
116
```
117
{ .api }
118

119
```typescript
120
// Multinomial coefficient: n! / (k1! * k2! * ... * km!)
121
// where n = k1 + k2 + ... + km
122

123
multinomial([2, 3, 4]) // 1260 = 9! / (2! * 3! * 4!)
124
multinomial([1, 1, 1]) // 6 = 3! / (1! * 1! * 1!)
125

126
// Number of ways to partition n objects into groups of sizes k1, k2, ..., km
127
// Example: Ways to divide 9 people into groups of 2, 3, and 4
128
multinomial([2, 3, 4]) // 1260 ways
129

130
// Relationship to binomial coefficients (special case)
131
multinomial([k, subtract(n, k)]) === combinations(n, k)
132
```
133

134
## Advanced Combinatorics
135

136
### Bell Numbers
137

138
```typescript
139
bellNumbers(n: MathType): MathType
140
```
141
{ .api }
142

143
```typescript
144
// Number of ways to partition a set of n elements
145
bellNumbers(0) // 1 (empty set has one partition)
146
bellNumbers(1) // 1 (one element: {1})
147
bellNumbers(2) // 2 (two elements: {{1},{2}} or {{1,2}})
148
bellNumbers(3) // 5 ({{1},{2},{3}}, {{1,2},{3}}, {{1,3},{2}}, {{2,3},{1}}, {{1,2,3}})
149
bellNumbers(4) // 15
150
bellNumbers(5) // 52
151

152
// Bell's triangle (similar to Pascal's triangle)
153
function bellTriangle(rows) {
154
  const triangle = [[1]]
155
  
156
  for (let i = 1; i < rows; i++) {
157
    const row = [triangle[i-1][triangle[i-1].length - 1]]
158
    
159
    for (let j = 1; j <= i; j++) {
160
      row.push(add(row[j-1], triangle[i-1][j-1]))
161
    }
162
    
163
    triangle.push(row)
164
  }
165
  
166
  return triangle
167
}
168
```
169

170
### Catalan Numbers
171

172
```typescript
173
catalan(n: MathType): MathType  
174
```
175
{ .api }
176

177
```typescript
178
// nth Catalan number: C_n = (2n)! / ((n+1)! * n!)
179
catalan(0) // 1
180
catalan(1) // 1
181
catalan(2) // 2  
182
catalan(3) // 5
183
catalan(4) // 14
184
catalan(5) // 42
185

186
// Applications:
187
// - Number of binary trees with n internal nodes
188
// - Number of ways to triangulate a polygon with n+2 sides  
189
// - Number of paths from (0,0) to (n,n) that don't cross the diagonal
190
// - Number of ways to parenthesize n+1 factors
191

192
// Recurrence relation: C_0 = 1, C_{n+1} = sum(C_i * C_{n-i}) for i=0 to n
193
function catalanRecurrence(n) {
194
  if (equal(n, 0)) return 1
195
  let sum = 0
196
  for (let i = 0; i < n; i++) {
197
    sum = add(sum, multiply(catalan(i), catalan(subtract(n, add(i, 1)))))
198
  }
199
  return sum
200
}
201
```
202

203
### Stirling Numbers of the Second Kind
204

205
```typescript
206
stirlingS2(n: MathType, k: MathType): MathType
207
```
208
{ .api }
209

210
```typescript
211
// Number of ways to partition n objects into k non-empty subsets
212
stirlingS2(4, 2) // 7 (ways to partition 4 objects into 2 groups)
213
stirlingS2(5, 3) // 25
214
stirlingS2(n, 1) // 1 (only one way: all objects in one subset)
215
stirlingS2(n, n) // 1 (only one way: each object in its own subset)
216

217
// Relationship to Bell numbers
218
function bellFromStirling(n) {
219
  let sum = 0
220
  for (let k = 0; k <= n; k++) {
221
    sum = add(sum, stirlingS2(n, k))
222
  }
223
  return sum
224
}
225
bellFromStirling(4) === bellNumbers(4) // true
226

227
// Recurrence: S(n,k) = k*S(n-1,k) + S(n-1,k-1)
228
```
229

230
### Compositions
231

232
```typescript
233
composition(n: MathType, k: MathType): MathType
234
```
235
{ .api }
236

237
```typescript
238
// Number of ways to write n as a sum of k positive integers (order matters)
239
composition(4, 2) // 3 (4 = 1+3, 2+2, 3+1)
240
composition(5, 3) // 6 (5 = 1+1+3, 1+2+2, 1+3+1, 2+1+2, 2+2+1, 3+1+1)
241

242
// Formula: C(n,k) = C(n-1, k-1) where C is combinations
243
// This is because we place k-1 dividers among n-1 positions
244
composition(n, k) === combinations(subtract(n, 1), subtract(k, 1))
245

246
// Total compositions of n (into any number of parts)
247
function totalCompositions(n) {
248
  return pow(2, subtract(n, 1))
249
}
250
totalCompositions(4) // 8 (2^3)
251
```
252

253
## Gamma Function and Related Functions
254

255
### Gamma Function
256

257
```typescript
258
gamma(n: MathType): MathType
259
```
260
{ .api }
261

262
```typescript
263
// Gamma function: Γ(n) = (n-1)! for positive integers
264
gamma(1) // 1 = 0!
265
gamma(2) // 1 = 1!
266
gamma(3) // 2 = 2!
267
gamma(5) // 24 = 4!
268

269
// For non-integers
270
gamma(0.5) // √π ≈ 1.772
271
gamma(1.5) // 0.5 * √π ≈ 0.886
272
gamma(2.5) // 1.5 * 0.5 * √π ≈ 1.329
273

274
// Complex numbers supported
275
gamma(complex(1, 1)) // Complex result
276

277
// Relationship to factorial
278
gamma(add(n, 1)) === factorial(n) // For non-negative integers
279

280
// Reflection formula: Γ(z) * Γ(1-z) = π / sin(πz)
281
multiply(gamma(z), gamma(subtract(1, z))) === divide(pi, sin(multiply(pi, z)))
282
```
283

284
### Log Gamma Function
285

286
```typescript
287
lgamma(n: MathType): MathType
288
```
289
{ .api }
290

291
```typescript
292
// Natural logarithm of the gamma function (more numerically stable)
293
lgamma(5) // log(24) ≈ 3.178 = log(gamma(5))
294
lgamma(100) // Much more stable than log(gamma(100))
295

296
// Useful for large arguments where gamma(n) would overflow
297
lgamma(bignumber('1000')) // Computable, whereas gamma(1000) is huge
298

299
// Relationship
300
lgamma(n) === log(gamma(n)) // But lgamma is more numerically stable
301

302
// Stirling's approximation for large n
303
function stirlingApprox(n) {
304
  return add(
305
    multiply(subtract(n, 0.5), log(n)),
306
    multiply(-1, n),
307
    multiply(0.5, log(multiply(2, pi)))
308
  )
309
}
310
// lgamma(n) ≈ stirlingApprox(n) for large n
311
```
312

313
## Random Number Generation
314

315
### Basic Random Numbers
316

317
```typescript
318
random(size?: number | number[], min?: number, max?: number): number | Matrix
319
```
320
{ .api }
321

322
```typescript
323
// Single random number [0, 1)
324
random() // e.g., 0.7234...
325

326
// Random number in range [min, max)
327
random(5, 10) // e.g., 7.8234...
328

329
// Array of random numbers
330
random(5) // [0.1234, 0.5678, ...] (5 numbers)
331
random([2, 3]) // 2x3 matrix of random numbers
332

333
// Seeded random (for reproducibility)
334
import { create, all, config } from 'mathjs'
335
const math = create(all, {
336
  randomSeed: 'deterministic-seed'
337
})
338
math.random() // Reproducible sequence
339
```
340

341
### Random Integers
342

343
```typescript
344
randomInt(size?: number | number[], min?: number, max?: number): number | Matrix  
345
```
346
{ .api }
347

348
```typescript
349
// Single random integer [min, max)
350
randomInt(1, 7) // Dice roll: 1, 2, 3, 4, 5, or 6
351
randomInt(0, 2) // Coin flip: 0 or 1
352

353
// Array of random integers
354
randomInt(10, 1, 11) // 10 numbers from 1 to 10
355
randomInt([3, 4], 0, 100) // 3x4 matrix of integers [0, 100)
356

357
// Uniform distribution over integers
358
function uniformInt(n, trials = 1000) {
359
  const counts = new Array(n).fill(0)
360
  for (let i = 0; i < trials; i++) {
361
    counts[randomInt(0, n)]++
362
  }
363
  return counts
364
}
365
```
366

367
### Pick Random Elements
368

369
```typescript
370
pickRandom(array: MathCollection, number?: number, weights?: MathCollection): MathType | MathCollection
371
```
372
{ .api }
373

374
```typescript
375
const items = ['apple', 'banana', 'cherry', 'date', 'elderberry']
376

377
// Pick single random element
378
pickRandom(items) // e.g., 'cherry'
379

380
// Pick multiple elements (with replacement)
381
pickRandom(items, 3) // e.g., ['banana', 'apple', 'banana']
382

383
// Weighted random selection
384
const weights = [0.1, 0.3, 0.4, 0.15, 0.05] // Must sum to 1
385
pickRandom(items, 1, weights) // More likely to pick 'cherry' (weight 0.4)
386

387
// Multiple weighted picks
388
pickRandom(items, 10, weights) // 10 picks with given probabilities
389

390
// Lottery/sampling applications
391
function lottery(tickets, winners) {
392
  return pickRandom(tickets, winners) // Pick winners from ticket holders
393
}
394

395
// Bootstrap sampling
396
function bootstrap(data, samples = 1000) {
397
  return range(0, samples).map(() => pickRandom(data, data.length))
398
}
399
```
400

401
## Probability Distributions (Custom Implementations)
402

403
### Discrete Distributions
404

405
```typescript
406
// Binomial probability mass function
407
function binomialPMF(k, n, p) {
408
  return multiply(combinations(n, k), multiply(pow(p, k), pow(subtract(1, p), subtract(n, k))))
409
}
410

411
// Poisson probability mass function
412
function poissonPMF(k, lambda) {
413
  return multiply(divide(pow(lambda, k), factorial(k)), exp(multiply(-1, lambda)))
414
}
415

416
// Geometric probability mass function
417
function geometricPMF(k, p) {
418
  return multiply(pow(subtract(1, p), subtract(k, 1)), p)
419
}
420

421
// Examples
422
binomialPMF(3, 10, 0.5) // P(X=3) in Binomial(10, 0.5)
423
poissonPMF(2, 3) // P(X=2) in Poisson(3)
424
geometricPMF(5, 0.2) // P(X=5) in Geometric(0.2)
425
```
426

427
### Continuous Distributions (approximations)
428

429
```typescript
430
// Standard normal PDF (approximation)
431
function normalPDF(x, mu = 0, sigma = 1) {
432
  const coefficient = divide(1, multiply(sigma, sqrt(multiply(2, pi))))
433
  const exponent = exp(divide(multiply(-0.5, pow(divide(subtract(x, mu), sigma), 2)), 1))
434
  return multiply(coefficient, exponent)
435
}
436

437
// Standard normal CDF (approximation using error function)
438
function normalCDF(x, mu = 0, sigma = 1) {
439
  const z = divide(subtract(x, mu), sigma)
440
  return multiply(0.5, add(1, erf(divide(z, sqrt(2)))))
441
}
442
```
443

444
## Information Theory and Divergence
445

446
### Kullback-Leibler Divergence
447

448
```typescript
449
kldivergence(q: MathCollection, p: MathCollection): MathType
450
```
451
{ .api }
452

453
```typescript
454
// KL divergence: D_KL(P||Q) = sum(p_i * log(p_i / q_i))
455
// Measures how one probability distribution differs from another
456

457
const p = [0.5, 0.3, 0.2] // True distribution  
458
const q = [0.4, 0.4, 0.2] // Approximating distribution
459

460
kldivergence(p, q) // KL divergence from P to Q
461

462
// Properties:
463
// - D_KL(P||Q) >= 0 (Gibbs' inequality)
464
// - D_KL(P||Q) = 0 iff P = Q
465
// - Generally D_KL(P||Q) ≠ D_KL(Q||P) (not symmetric)
466

467
// Cross entropy
468
function crossEntropy(p, q) {
469
  return add(entropy(p), kldivergence(p, q))
470
}
471

472
// Entropy (custom implementation)
473
function entropy(p) {
474
  return multiply(-1, sum(p.map(x => multiply(x, log(x)))))
475
}
476
```
477

478
## Set Operations (Combinatorial Applications)
479

480
### Basic Set Operations
481

482
```typescript
483
setUnion(a1: MathCollection, a2: MathCollection): MathCollection
484
setIntersect(a1: MathCollection, a2: MathCollection): MathCollection  
485
setDifference(a1: MathCollection, a2: MathCollection): MathCollection
486
setSymDifference(a1: MathCollection, a2: MathCollection): MathCollection
487
```
488
{ .api }
489

490
```typescript
491
const A = [1, 2, 3, 4]
492
const B = [3, 4, 5, 6]
493

494
setUnion(A, B) // [1, 2, 3, 4, 5, 6]
495
setIntersect(A, B) // [3, 4]
496
setDifference(A, B) // [1, 2] (A - B)
497
setSymDifference(A, B) // [1, 2, 5, 6] (A ∪ B) - (A ∩ B)
498

499
// Principle of inclusion-exclusion
500
setSize(setUnion(A, B)) === add(
501
  setSize(A), 
502
  setSize(B), 
503
  multiply(-1, setSize(setIntersect(A, B)))
504
)
505
```
506

507
### Advanced Set Operations
508

509
```typescript
510
setCartesian(a1: MathCollection, a2: MathCollection): MathCollection
511
setPowerset(a: MathCollection): MathCollection
512
setDistinct(a: MathCollection): MathCollection
513
```
514
{ .api }
515

516
```typescript
517
const A = [1, 2]
518
const B = ['a', 'b', 'c']
519

520
// Cartesian product
521
setCartesian(A, B) // [[1,'a'], [1,'b'], [1,'c'], [2,'a'], [2,'b'], [2,'c']]
522
setSize(setCartesian(A, B)) === multiply(setSize(A), setSize(B)) // 2 * 3 = 6
523

524
// Power set (all subsets)
525
setPowerset([1, 2, 3]) // [[], [1], [2], [3], [1,2], [1,3], [2,3], [1,2,3]]
526
setSize(setPowerset(A)) === pow(2, setSize(A)) // 2^n subsets
527

528
// Remove duplicates
529
setDistinct([1, 2, 2, 3, 3, 3]) // [1, 2, 3]
530
```
531

532
### Set Properties and Tests
533

534
```typescript
535
setIsSubset(a1: MathCollection, a2: MathCollection): boolean
536
setMultiplicity(e: MathType, a: MathCollection): number
537
setSize(a: MathCollection): number
538
```
539
{ .api }
540

541
```typescript
542
const A = [1, 2, 3]
543
const B = [1, 2, 3, 4, 5]
544

545
setIsSubset(A, B) // true (A ⊆ B)
546
setIsSubset(B, A) // false
547

548
setMultiplicity(2, [1, 2, 2, 2, 3]) // 3 (element 2 appears 3 times)
549
setSize([1, 2, 3, 4, 5]) // 5 (cardinality)
550

551
// Set equality
552
function setEqual(A, B) {
553
  return setSize(setSymDifference(A, B)) === 0
554
}
555
```
556

557
## Probability Applications
558

559
### Simulation and Monte Carlo
560

561
```typescript
562
// Simulate coin flips
563
function coinFlips(n, p = 0.5) {
564
  return range(0, n).map(() => random() < p ? 1 : 0)
565
}
566

567
// Estimate π using Monte Carlo
568
function estimatePi(samples = 1000000) {
569
  let insideCircle = 0
570
  
571
  for (let i = 0; i < samples; i++) {
572
    const x = random(-1, 1)
573
    const y = random(-1, 1) 
574
    if (add(pow(x, 2), pow(y, 2)) < 1) {
575
      insideCircle++
576
    }
577
  }
578
  
579
  return multiply(4, divide(insideCircle, samples))
580
}
581

582
// Birthday paradox simulation
583
function birthdayParadox(people = 23, trials = 10000) {
584
  let matches = 0
585
  
586
  for (let trial = 0; trial < trials; trial++) {
587
    const birthdays = randomInt(people, 1, 366)
588
    const unique = setDistinct(birthdays)
589
    if (setSize(unique) < people) {
590
      matches++
591
    }
592
  }
593
  
594
  return divide(matches, trials)
595
}
596
```
597

598
### Probability Mass Function Calculations
599

600
```typescript
601
// Expected value for discrete distribution
602
function expectedValue(values, probabilities) {
603
  return sum(values.map((val, i) => multiply(val, probabilities[i])))
604
}
605

606
// Variance for discrete distribution  
607
function varianceDiscrete(values, probabilities) {
608
  const mu = expectedValue(values, probabilities)
609
  return sum(values.map((val, i) => 
610
    multiply(probabilities[i], pow(subtract(val, mu), 2))
611
  ))
612
}
613

614
// Probability generating function
615
function probabilityGeneratingFunction(probabilities, z) {
616
  return sum(probabilities.map((p, k) => multiply(p, pow(z, k))))
617
}
618

619
// Moment generating function
620
function momentGeneratingFunction(values, probabilities, t) {
621
  return sum(values.map((val, i) => 
622
    multiply(probabilities[i], exp(multiply(t, val)))
623
  ))
624
}
625
```
626

627
### Combinatorial Probability
628

629
```typescript
630
// Classic problems using combinatorics
631

632
// Hypergeometric distribution
633
function hypergeometric(k, N, K, n) {
634
  // k successes in n draws from population of N with K successes
635
  return divide(
636
    multiply(combinations(K, k), combinations(subtract(N, K), subtract(n, k))),
637
    combinations(N, n)
638
  )
639
}
640

641
// Derangements probability
642
function derangementProb(n) {
643
  return divide(subfactorial(n), factorial(n))
644
}
645

646
// Matching problem (random permutation has fixed points)
647
function matchingProb(n, k) {
648
  // Probability that exactly k elements are in their original positions
649
  return divide(
650
    multiply(combinations(n, k), subfactorial(subtract(n, k))),
651
    factorial(n)
652
  )
653
}
654
```
655

656
## Performance and Numerical Considerations
657

658
### Large Number Handling
659

660
```typescript
661
// Use BigNumber for large factorials and combinations
662
const largeCombination = combinations(bignumber('1000'), bignumber('500'))
663

664
// Use lgamma for very large gamma function values
665
const largeGamma = exp(lgamma(bignumber('1000'))) // More stable than gamma(1000)
666

667
// Stirling's approximation for very large factorials
668
function stirlingFactorial(n) {
669
  return multiply(
670
    sqrt(multiply(2, pi, n)),
671
    pow(divide(n, e), n)
672
  )
673
}
674
```
675

676
### Numerical Stability
677

678
```typescript
679
// Use log-space calculations for small probabilities
680
function logBinomialPMF(k, n, p) {
681
  return add(
682
    lgamma(add(n, 1)),
683
    multiply(-1, lgamma(add(k, 1))),
684
    multiply(-1, lgamma(subtract(n, k, -1))),
685
    multiply(k, log(p)),
686
    multiply(subtract(n, k), log(subtract(1, p)))
687
  )
688
}
689

690
// Convert back to probability space
691
const logProb = logBinomialPMF(3, 10, 0.5)
692
const prob = exp(logProb)
693
```
694

695
## Chain Operations and Functional Programming
696

697
```typescript
698
// Chain operations with probability functions
699
const diceRolls = chain(range(1, 7))
700
  .map(face => divide(1, 6)) // Equal probability for each face
701
  .done()
702

703
// Functional approach to combinatorial calculations  
704
const pascalRow = n => 
705
  range(0, add(n, 1))
706
    .map(k => combinations(n, k))
707
    
708
const fibonacciFromBinomial = n =>
709
  range(0, add(n, 1))
710
    .map(k => combinations(subtract(n, k), k))
711
    .reduce((a, b) => add(a, b), 0)
712
```