如何使用 scipy.optimize.fmin_l_bfgs_b 拟合递归关系？答案

【问题标题】：How do I fit recurrence relations using scipy.optimize.fmin_l_bfgs_b?如何使用 scipy.optimize.fmin_l_bfgs_b 拟合递归关系？
【发布时间】：2020-02-10 05:13:08
【问题描述】：

我最近观看了Veritasium 的video，他在其中对logistic map 进行了有趣的讨论。这让我开始思考递推关系，以及参数递推关系如何适合数据。

我想参数化 $\theta$ in $\hat y_{n+1} = \theta \hat y_{n}(1-\hat y_{n})$ 在序列 $y_k \in [ 0, 1]$ 其中$k=n+1$ 通过使用L-BFGS-B algorithm 最小化均方误差Scipy。这个例子对我推广以将其他递归关系拟合到真实世界的数据具有指导意义。如何实现目标函数，其中预测值是递归关系的输出，以便可以将其传递给 fmin_l_bfgs_b 的 func 参数？

【问题讨论】：

我打算发布这个问题以进行交叉验证，但我觉得这与编码有关，而不是统计/机器学习。我看到我的 LaTeX 代码没有翻译...

标签： python scipy mathematical-optimization numerical-methods recurrence

【解决方案1】：

如果我正确理解您的问题，您想找到 $\theta$ 的值，以最小化 $\sum_{n=0}^{N-1} (\hat y_{n+1} - \theta \hat y_{n} (1 - \hat y_{n}))^2$，给定一个序列 ${\hat y_k}_{k=0}^N$。如果是这样，假设您的数据是 ys 并且您最初的猜测是 x0，您可以通过

def f(l): 
    t = l[0] 
    return ((ys[1:] - (t * ys[:-1] * (1 - ys[:-1])))**2).sum()

fmin_l_bfgs_b(f, x0=(x0,), approx_grad=True)

例如，如果我们创建一些 theta 约为 3 的数据：

In [43]: import numpy as np 
    ...: ys = [0.3] 
    ...: theta = 3 
    ...: for _ in range(100): 
    ...:     ys.append((np.random.uniform(-0.02, 0.02) + theta)*ys[-1] * (1 - ys[-1])) 
    ...: ys = np.array(ys) 
    ...:                                                                                        

In [44]: def f(l): 
    ...:     t = l[0] 
    ...:     return ((ys[1:] - (t * ys[:-1] * (1 - ys[:-1])))**2).sum() 
    ...: fmin_l_bfgs_b(f, x0=(0.5,), approx_grad=True)                                          
Out[44]: 
(array([2.99949454]),
 0.0006258145273212467,
 {'grad': array([-5.70908338e-07]),
  'task': b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL',
  'funcalls': 6,
  'nit': 2,
  'warnflag': 0})

当然，您也可以在这里提供渐变；我只是有点懒。

但是，如果这确实是您想要做的事情，那么您很可能需要针对最小二乘问题量身定制的东西（例如 Levenberg--Marquardt）；在 SciPy 中，这些方法在 scipy.optimize.least_squares 中可用。有了这些，您的问题归结为以下几点：

def F(t): 
    return ys[1:] - (t * ys[:-1] * (1 - ys[:-1]))

least_squares(F, x0=x0)

根据上面的数据：

In [53]: def F(t): 
    ...:     return ys[1:] - (t * ys[:-1] * (1 - ys[:-1])) 
    ...:                                                                                        

In [54]: least_squares(F, x0=0.5)                                                               
Out[54]: 
 active_mask: array([0.])
        cost: 0.00031290726365087974
         fun: ...
        grad: array([-2.43698605e-09])
         jac: ...
     message: '`gtol` termination condition is satisfied.'
        nfev: 4
        njev: 4
  optimality: 2.4369860459044074e-09
      status: 1
     success: True
           x: array([2.9994946])

【讨论】：