python / scipy 中的牛顿法

Question

我正在尝试在 python 上使用牛顿算法获取函数的根。即使更改精度级别，也会出现运行时错误。你能帮我了解一下我该如何改进吗？

最好的，国标

下面是我的'simple'代码和求根部分：

from scipy.stats import norm
import numpy as np 
import matplotlib.pyplot as plt  
import scipy.optimize as opt


def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

x0=112
r=0
T=1
d=0
sigma=0.2
q=0
K1=110
K2=130




res2= opt.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)


the error i get is: res2= opt.zeros.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)
/Users/anaconda/lib/python3.5/site-packages/scipy/optimize/zeros.py:173: RuntimeWarning: Tolerance of 0.011300000000005639 reached
  warnings.warn(msg, RuntimeWarning)

Answer 1

很难在如此稀疏的上下文中提供建议。但有些评论：

甲：

您没有按照描述使用 newton here:

The Newton-Raphson method is used if the derivative fprime of func is provided, otherwise the secant method is used.

乙：

你的错误来自 here 我会说：

由于这些 tol 值是硬编码的，我想这不应该发生！

# Secant method
p0 = x0
if x0 >= 0:
    p1 = x0*(1 + 1e-4) + 1e-4
else:
    p1 = x0*(1 + 1e-4) - 1e-4
q0 = func(*((p0,) + args))
q1 = func(*((p1,) + args))
for iter in range(maxiter):
    if q1 == q0:
        if p1 != p0:
            msg = "Tolerance of %s reached" % (p1 - p0)
            warnings.warn(msg, RuntimeWarning)
        return (p1 + p0)/2.0

意思是：你的代码可能有问题！

C:

让我们试试更慢但更安全的方法bisection-method:

# brackets not tuned! it's just some trial!
res2 = opt.bisect(vega_callspread, 0, 200, args=(r,T,d,sigma,q,K1,K2))

输出：

X:\so_newton.py:9: RuntimeWarning: divide by zero encountered in log
  d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
X:\so_newton.py:10: RuntimeWarning: divide by zero encountered in log
  d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

这是个坏兆头。

D（跟随 C）：

您的函数读取：

def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

你打电话给：

x0=112
r=0
T=1
d=0
sigma=0.2
q=0
K1=110
K2=130

args=(r,T,d,sigma,q,K1,K2,)

没有S！

所以要么你正在危险地重命名变量，要么你这边出现了一些错误！

S 的值是 log(S/K1) 和 co.

中的问题

这不是关于 S/K1 而是关于这个：

import numpy as np
np.log(0)
# __main__:1: RuntimeWarning: divide by zero encountered in log
# -inf

我从 .

得到的

您现在可以尝试解释这对您的函数做了什么，因为此 log-eval 将获得 -np.inf 的值。

E scipy 如何处理这个问题（遵循 D）：

我懒得阅读 args-handling 的 docs/sources，但让我们检查一下（向 func 添加打印；像之前一样使用二分法）：

def vega_callspread(r,S,T,d,sigma,q,K1,K2):

    print('S: ', S)

    d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
    d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))

    u1 = S*np.exp(-q*T)*norm.pdf(d1,0,1)
    u2 = S*np.exp(-q*T)*norm.pdf(d2,0,1)

    return u1-u2;

输出：

S:  0
X:\so_newton.py:11: RuntimeWarning: divide by zero encountered in log
  d1 = (np.log(S/K1)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
X:\so_newton.py:12: RuntimeWarning: divide by zero encountered in log
  d2 = (np.log(S/K2)+(r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
S:  0

所以 arg-handling 似乎是通过 arg-name 进行的（而不仅仅是调用中的参数位置；我在这里遗漏了正确的编程语言术语；再次懒惰！）

这个 (S=0) 可能非常糟糕（而不是你想要的），这是你的错误！

编辑

根据您的评论，您似乎尝试优化 S。这让我很清楚，你正在使用一些优化算法优化 x 而你的函数中没有 x！

我不是在这里分析你的任务，但你可能想让你的函数使用一些 x（由 x0 初始化），因为这是 [=29= 的一般想法].我们可以保留名称 S，但它需要是函数的第一个参数。这在文档中都有解释。

所以：

def vega_callspread(S, r,T,d,sigma,q,K1,K2):  # S now first argument !!!
    ...

res2= opt.newton(vega_callspread, x0, args=(r,T,d,sigma,q,K1,K2,),tol=10**(-1),maxiter=1000000)  # S removed from args; S is init by x0 -> read docs!

输出：

117.214682594

python / scipy 中的牛顿法

Newton method in python / scipy

python

mathematical-optimization

scipy

newtons-method

甲：

乙：

C:

D（跟随 C）：

E scipy 如何处理这个问题（遵循 D）：