SciPy - 特殊函数

SciPy - 特殊函数 (scipy.special)

scipy.special 包是数学特殊函数的综合集合，这些函数常在物理、工程、统计和其他科学学科中遇到。其中许多函数被实现为 NumPy 的通用函数（ufuncs），这意味着它们可以在 NumPy 数组上按元素进行操作，并支持 广播（broadcasting）。

让我们探索一些常用的特殊函数：

• 立方根函数（cbrt）
• 指数函数（例如，exp10、expm1、expn）
• 对数函数（例如，logsumexp、loggamma）
• 误差函数（例如，erf、erfc）
• 伽马函数（gamma）及相关函数（例如，gammaln、digamma）
• 贝塞尔函数（例如，jv、jn、kv、kn）
• 组合学：排列（perm）和组合（comb）
• 兰伯特W函数（lambertw）
• 以及许多其他函数，包括艾里函数（Airy functions）、斯特鲁夫函数（Struve functions）、超几何函数（hypergeometric functions）等。

我们将简要介绍几个示例。

立方根函数（cbrt）

scipy.special.cbrt(x) 函数计算 x 的按元素的实数立方根。

from scipy.special import cbrt
import numpy as np

values = np.array([8.0, 27.0, -64.0, 0.125])
cube_roots = cbrt(values)

print(f"Original values: {values}")
print(f"Cube roots: {cube_roots}")

输出：

Original values: [  8.    27.   -64.     0.125]
Cube roots: [ 2.  3. -4.  0.5]

指数函数（exp10）

scipy.special.exp10(x) 函数计算 $10^x$ 的按元素结果。

from scipy.special import exp10
import numpy as np

exponents = np.array([0, 1, 2, -1, 3.5])
results_exp10 = exp10(exponents)

print(f"Exponents: {exponents}")
print(f"10^x results: {results_exp10}")

输出：

Exponents: [ 0.   1.   2.  -1.   3.5]
10^x results: [1.00000000e+00 1.00000000e+01 1.00000000e+02 1.00000000e-01
 3.16227766e+03]

相对误差指数函数（exprel）

scipy.special.exprel(x) 函数计算 $(\exp(x) - 1)/x$。这种形式对于 $x$ 的小值特别有用，因为直接计算 $\exp(x) - 1$ 可能会发生灾难性抵消（precision loss，精度损失）。exprel 提供了一个数值上更稳定的结果。

from scipy.special import exprel
import numpy as np

x_values = np.array([-0.2, -0.01, 1e-9, 0.01, 0.2])
results_exprel = exprel(x_values)

print(f"x values: {x_values}")
print(f"(e^x - 1)/x results: {results_exprel}")

# For comparison, direct computation for a very small x:
x_small = 1e-9
direct_comp = (np.exp(x_small) - 1) / x_small
print(f"\nFor x_small = {x_small}:")
print(f"  exprel(x_small) = {exprel(x_small)}")
print(f"  Direct computation = {direct_comp}")

输出：

x values: [-2.0e-01 -1.0e-02  1.0e-09  1.0e-02  2.0e-01]
(e^x - 1)/x results: [0.90634623 0.99501663 1.         1.00501671 1.10701379]

For x_small = 1e-09:
  exprel(x_small) = 1.0
  Direct computation = 1.0000000005

对数和指数函数（logsumexp）

scipy.special.logsumexp(a, axis=None, b=None, keepdims=False, return_sign=False) 函数计算 $\log(\sum \exp(a_i))$。这在机器学习中非常有用（例如，计算对数似然），因为它可以帮助避免在计算指数和时出现下溢（underflow）或上溢（overflow）问题，尤其对于非常大或非常小的数。

from scipy.special import logsumexp
import numpy as np

large_numbers = np.array([1000, 1001, 999])
result_lse = logsumexp(large_numbers)

print(f"Numbers: {large_numbers}")
print(f"log(sum(exp(numbers))): {result_lse}")

# Direct computation might overflow:
# print(np.log(np.sum(np.exp(large_numbers)))) # This would likely cause an OverflowError

输出：

Numbers: [1000 1001  999]
log(sum(exp(numbers))): 1001.4076059644443

兰伯特W函数（lambertw）

scipy.special.lambertw(z, k=0, tol=1e-8) 函数计算兰伯特W函数，$W(z)$，它是 $f(w) = w \exp(w)$ 的反函数。也就是说，$z = W(z) \exp(W(z))$。兰伯特W函数是多值函数；参数 k 用于选择分支（k=0 是主分支，principal branch）。

from scipy.special import lambertw
import numpy as np

z_value = 1.0
w_principal = lambertw(z_value) # Principal branch k=0
w_k_minus_1 = lambertw(z_value, k=-1) # Another branch

print(f"For z = {z_value}:")
print(f"  Principal branch W_0(z): {w_principal}")
print(f"  Check W_0(z) * exp(W_0(z)): {w_principal * np.exp(w_principal)}")
# print(f"  Branch W_-1(z): {w_k_minus_1}") # This might be complex for z=1

输出（$W_0(1)$ 的值也称为 Omega 常数）：

For z = 1.0:
  Principal branch W_0(z): (0.5671432904097838+0j)
  Check W_0(z) * exp(W_0(z)): (1.0+0j)

排列（perm）和组合（comb）

这些函数计算组合学量：

组合： scipy.special.comb(N, k, exact=False, repetition=False) 计算从 N 个项目中选取 k 个项目的方式数量，不考虑选择的顺序。这就是“N选k”或 $\binom{N}{k}$。

• exact=True 返回精确的整数结果（可以处理大数）。
• repetition=True 允许重复选取项目。

from scipy.special import comb

# How many ways to choose 3 items from 10 (no repetition)
num_combinations = comb(10, 3, exact=True)
print(f"Combinations C(10, 3) without repetition: {num_combinations}")

# How many ways to choose 3 items from 10 (with repetition allowed)
num_combinations_rep = comb(10, 3, exact=True, repetition=True)
print(f"Combinations C(10, 3) with repetition: {num_combinations_rep}")

输出：

Combinations C(10, 3) without repetition: 120
Combinations C(10, 3) with repetition: 220

排列： scipy.special.perm(N, k, exact=False) 计算从 N 个项目中选取 k 个项目并进行排列的方式数量，此时选择的顺序很重要。这就是 $P(N, k) = N! / (N-k)!$。

from scipy.special import perm

# How many ways to arrange 3 items chosen from 10
num_permutations = perm(10, 3, exact=True)
print(f"Permutations P(10, 3): {num_permutations}")

输出：

Permutations P(10, 3): 720

伽马函数（gamma）

scipy.special.gamma(z) 函数计算伽马函数 $\Gamma(z)$，它是阶乘函数向实数和复数的扩展。对于正整数 $n$，$\Gamma(n) = (n-1)!$。

from scipy.special import gamma
import numpy as np

values_gamma = np.array([1, 2, 3, 4, 5, 0.5])
gamma_results = gamma(values_gamma)

print(f"Input values z: {values_gamma}")
print(f"Gamma(z) results: {gamma_results}")
# Note: gamma(n) = (n-1)!  So gamma(5) = 4! = 24
# gamma(0.5) = sqrt(pi)
print(f"gamma(0.5) = {gamma(0.5):.6f}, sqrt(pi) = {np.sqrt(np.pi):.6f}")

输出：

Input values z: [1.  2.  3.  4.  5.  0.5]
Gamma(z) results: [1.00000000e+00 1.00000000e+00 2.00000000e+00 6.00000000e+00
 2.40000000e+01 1.77245385e+00]
gamma(0.5) = 1.772454, sqrt(pi) = 1.772454

scipy.special 模块非常庞大。当你遇到 NumPy 或标准 Python math 库中没有的数学函数时，很可能可以在这里找到它。请始终查阅 官方 SciPy 文档 以获取完整函数列表和详细说明。