Sympy函数

Question

到目前为止，我还没有使用过渐近，但我想尝试一下，因为它似乎是定义您想要优化的函数的一个很好的方法。

，这就是我正在尝试用sympy编写的函数。我知道我可以简单地使用一个普通的函数来做这件事，但是这将是很酷的，例如，在渐近中编写它，然后它将非常容易，例如，我认为python可以找到导数。我要把这笔钱算出来

再说一次，使用普通函数很简单，但如果有人能帮助我在渐变中编写这个函数，我会很感激！总而言之，我想要编写这个函数，并具有灵活性，让和运行在一个向量上，如我提供的示例。另外，我这样做基本上是为了优化这个函数，使用类似牛顿的拉普森方法。这只是一个了解算法的玩具示例。注意: theta是一个未知数，我将对其进行优化。谢谢你的提示：)！

是

Answer 1

你当然可以使用sympy来找出导数：

In [132]: x = IndexedBase('x', real=True)                                                                                                      

In [133]: n, i = symbols('n, i', integer=True, positive=True)                                                                                  

In [134]: theta = Symbol('theta', real=True)                                                                                                   

In [135]: objective = -n*log(pi) - Sum(log(1 + (x[i] - theta)**2), (i, 0, n-1))                                                                

In [136]: objective                                                                                                                            
Out[136]: 
            n - 1                      
             ___                       
             ╲                         
              ╲      ⎛           2    ⎞
-n⋅log(π) -   ╱   log⎝(-θ + x[i])  + 1⎠
             ╱                         
             ‾‾‾                       
            i = 0                      

In [137]: objective.diff(theta)                                                                                                                
Out[137]: 
 n - 1                 
  ____                 
  ╲                    
   ╲     2⋅θ - 2⋅x[i]  
    ╲  ────────────────
-   ╱             2    
   ╱   (-θ + x[i])  + 1
  ╱                    
  ‾‾‾‾                 
 i = 0

您还可以对这些表达式进行lambdify以加快计算速度(请注意，我将总和从0到n-1以匹配Python索引)：

In [138]: lambdify((theta, x, n), objective)(1, [1,2,3], 3)                                                                                    
Out[138]: -5.736774750542246

您还可以尝试对小输入样本进行解析求解，例如：

In [142]: vec = [1, 2, 3]                                                                                                                      

In [143]: objective.diff(theta).subs(n, 3).doit().subs({x[i]:vec[i] for i in range(3)})                                                        
Out[143]: 
    2⋅θ - 6        2⋅θ - 4        2⋅θ - 2   
- ──────────── - ──────────── - ────────────
         2              2              2    
  (3 - θ)  + 1   (2 - θ)  + 1   (1 - θ)  + 1

In [144]: solve(_, theta)                                                                                                                      
Out[144]: 
⎡           _____________          _____________          _____________          _____________⎤
⎢          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ ⎥
⎢2, 2 -   ╱  - ─ - ───── , 2 +   ╱  - ─ - ───── , 2 -   ╱  - ─ + ───── , 2 +   ╱  - ─ + ───── ⎥
⎣       ╲╱     3     3         ╲╱     3     3         ╲╱     3     3         ╲╱     3     3   ⎦

对于较大的输入，解析解将不起作用(方程是多项式的，因此Abel-Ruffini对此进行了限制)，但如果您愿意，可以使用渐近进行数值求解：

In [150]: vec = [1, 2, 3, 4, 5]                                                                                                                

In [151]: objective.diff(theta).subs(n, 5).doit().subs({x[i]:vec[i] for i in range(5)})                                                        
Out[151]: 
    2⋅θ - 10       2⋅θ - 8        2⋅θ - 6        2⋅θ - 4        2⋅θ - 2   
- ──────────── - ──────────── - ──────────── - ──────────── - ────────────
         2              2              2              2              2    
  (5 - θ)  + 1   (4 - θ)  + 1   (3 - θ)  + 1   (2 - θ)  + 1   (1 - θ)  + 1

In [152]: nsolve(_, theta, 1)                                                                                                                  
Out[152]: 3.00000000000000

对于较大的输入数据，使用lambdified的目标函数和派生函数以及类似scipy's minimize的函数可以获得更快的速度：

In [158]: f = lambdify((theta, x, n), -objective)                                                                                              

In [159]: fp = lambdify((theta, x, n), -objective.diff(theta))                                                                                 

In [160]: from scipy.optimize import minimize                                                                                                  

In [161]: minimize(f, 1, jac=fp, args=([1, 2, 3], 3))                                                                                          
Out[161]: 
      fun: 4.820484018668093
 hess_inv: array([[0.50002255]])
      jac: array([9.77560778e-08])
  message: 'Optimization terminated successfully.'
     nfev: 4
      nit: 3
     njev: 4
   status: 0
  success: True
        x: array([2.00000005])

请注意，结果并不像使用nsolve计算的结果那样准确。

此外，这种情况下的lambdified函数没有手写的numpy代码那么有效，因为它没有向量化总和：

In [169]: print(inspect.getsource(f))                                                                                                          
def _lambdifygenerated(theta, Dummy_167, n):
    return (n*log(pi) + (builtins.sum(log((-theta + Dummy_167[i])**2 + 1) for i in range(0, n - 1+1))))