如何用符号回归近似拟合帕斯卡三角的某一行？

你遇到的问题核心确实和模型缺少幂运算能力有关——毕竟帕斯卡三角的每一个元素本质是组合数 ( C(n,k) = \frac{n!}{k!(n-k)!} )，不管是直接表达这个公式，还是用斯特林近似去逼近阶乘，都离不开幂运算。下面给你几个可行的解决方案，一步步帮你提升拟合效果：

1. 先给模型加上幂运算符 #

PySR底层用Julia实现，所以幂运算符是^（不是Python的**），你只需要把它添加到binary_operators列表里就行。修改后的模型初始化部分会是这样：

model = PySRRegressor(
    maxsize=15,
    niterations=5000,
    binary_operators=["+", "*", "^"],  # 新增幂运算符
    unary_operators=[
        "log",
        "exp",
        "inv",
        "square",
        "sqrt",
        "sign",
    ],
    elementwise_loss="loss(prediction, target) = (prediction - target)^2",
)

添加这个运算符后，模型就能搜索包含幂运算的表达式，这是逼近阶乘或组合数的基础。

2. 给模型提供更有针对性的特征 #

你当前只用了k（也就是0到20的整数）作为输入特征，但组合数是依赖n（这里是20）和k的关系的。不如构造更丰富的特征，帮模型更快捕捉到组合数的结构：

n = 20
y = print_pascals_triangle_row(n)
# 构造特征：k值、n值、n-k值
X = np.array([[k, n, n - k] for k in range(len(y))])

这样模型更容易发现n、k、n-k之间的乘积/阶乘关系，进而拟合出组合数公式。

3. 注册自定义阶乘运算符（进阶技巧） #

既然目标是组合数，直接给模型提供阶乘运算符能大幅降低搜索难度。PySR支持注册自定义运算符，我们可以手动添加阶乘：

from pysr import register_operator
import math

# 定义Python侧的阶乘实现
def factorial(x):
    return [math.factorial(int(round(i))) for i in x]

# 注册到PySR（需要对应Julia和SymPy的定义）
register_operator(
    "factorial",
    jl_definition="factorial(x) = factorial(Int(round(x)))",
    sympy_definition="factorial",
    arity=1,
)

然后把"factorial"加到unary_operators里，模型就能直接调用阶乘函数，很大概率会直接拟合出factorial(n) / (factorial(k) * factorial(n - k))这个精确的组合数表达式。

4. 调整搜索参数提升拟合质量 #

你当前的maxsize=15（最大表达式长度）和niterations=5000（迭代次数）可能不够搜索到复杂的组合数表达式。可以适当调大这两个参数：

• maxsize调到25左右，允许模型生成更长的表达式
• niterations调到10000甚至更高，给模型足够的搜索时间

另外打开verbosity=1可以看到搜索过程，方便你观察模型的进展。

修改后的完整示例代码 #

import math
import numpy as np
from pysr import PySRRegressor, register_operator

# 注册自定义阶乘运算符
def factorial(x):
    return [math.factorial(int(round(i))) for i in x]

register_operator(
    "factorial",
    jl_definition="factorial(x) = factorial(Int(round(x)))",
    sympy_definition="factorial",
    arity=1,
)

# 生成帕斯卡三角第n行
def print_pascals_triangle_row(n):
    row = []
    for k in range(n + 1):
       coefficient = math.comb(n, k)
       row.append(coefficient)
    return row

n = 20
y = print_pascals_triangle_row(n)
# 构造多特征输入
X = np.array([[k, n, n - k] for k in range(len(y))])

# 初始化优化后的模型
model = PySRRegressor(
    maxsize=25,
    niterations=10000,
    binary_operators=["+", "*", "^"],
    unary_operators=[
        "log",
        "exp",
        "inv",
        "square",
        "sqrt",
        "sign",
        "factorial"
    ],
    elementwise_loss="loss(prediction, target) = (prediction - target)^2",
    verbosity=1
)

# 拟合模型
model.fit(X, y)

# 查看最优表达式
model.show()

训练完成后，你可以用model.show()查看找到的最优表达式，大概率会得到和组合数公式高度吻合的结果，或者精度很高的近似式。

备注：内容来源于stack exchange，提问作者Simd