标签: numerical-integration

我想计算曲线下面积来进行积分而不定义像in这样的函数integrate().

我的数据看起来像这样:

Date          Strike     Volatility
2003-01-01    20         0.2
2003-01-01    30         0.3
2003-01-01    40         0.4
etc.

我试图plot(strike, volatility)看看波动性微笑.有没有办法整合这个绘制的"曲线"？

计划目的:整合.我正在实现高维度(高达100)的自适应正交(又称数值积分)算法.该想法是通过使用与该点处的误差估计成比例的采样密度来评估点来随机地将体积分成更小的部分.在早期我"老化"一个统一的样本,然后根据估计误差的高斯分布随机选择点.以类似于模拟退火的方式,I"降低温度"并随着时间的推移降低高斯的标准偏差,因此低误差点最初有很大的选择机会,但后来选择稳定下降可能性.这使程序能够偶然发现由于错误功能的缺陷而可能错过的尖峰.(我的算法在精神上类似于马尔可夫链蒙特卡罗积分.)

功能特征.要整合的功能是由于自然灾害导致的多个建筑物的保险政策损失估计.政策功能并不顺利:有免赔额,最高限额,分层数(例如,零赔偿支出高达100万美元,100%支付从1-2百万美元,然后零支付超过200万美元)和其他奇怪的政策条款.这引入了非线性行为和在许多平面中没有衍生物的函数.在政策功能之上是损坏功能,它根据建筑类型和飓风强度而变化,绝对不是钟形.

问题上下文:错误功能.困难在于选择一个好的错误功能.对于每个点,我记录对此有用的度量:函数的大小,由于先前的度量(一阶导数的代理)而改变了多少,该点占据的区域的体积(更大的体积可以更好地隐藏错误),以及与区域形状相关的几何因子.我的误差函数将是这些度量的线性组合,其中每个度量被赋予不同的权重.(如果我得到的结果很差,我会考虑非线性函数.)为了帮助我完成这项工作,我决定对每个权重的各种可能值进行优化,因此微软解决方案基金会.

要优化的内容:错误等级.我的测量标准化,从零到一.随着集成的进行以反映函数值,变化等的近期平均值,这些错误值会逐步修改.因此,我不是要创建一个产生实际错误值的函数,而是生成一个对其进行排序的数字.真正的错误,即如果所有采样点都按此估计误差值排序,则它们应该接收与按真实误差排序时所接收的等级相似的等级.

并非所有点都相同.如果#1真正错误的点区域排名#1000(反之亦然),我非常在意,但如果#500点排名#1000,则非常小心.我衡量成功的标准是在算法执行的中途将许多区域的以下总和最小化:

ABS(Log2(trueErrorRank) - Log2(estimatedErrorRank))

对于Log2,我使用的函数返回小于或等于数字的最大2的幂.从这个定义来看,得到有用的结果.交换#1和#2需要花费一分,但交换#2和#3不需要任何费用.这具有将点分为两个范围的幂的效果.在范围内交换的点不会添加到该函数中.

我如何评估.我已经构建了一个名为Rank的类来执行此操作:

1. 按真实误差对所有地区排名一次.

2. 对于每组独立的参数化权重,它计算该区域的试验(估计)误差.

3. 按试验错误对区域进行排序.

4. 计算每个地区的试验等级.

5. 添加两个等级的日志的绝对差异,并将其称为参数化的值,因此要最小化的值.

C#代码.完成所有这些后,我只需要一种方法来设置Microsoft Solver Foundation以找到最佳参数.语法让我难过.这是我到目前为止的C#代码.在其中,您将看到我已经确定的三个问题的评论.也许你可以发现更多!任何想法如何使这项工作？

public void Optimize()
{
    // Get the parameters from the GUI and figures out the low and high values for each weight.
    ParseParameters();

    // Computes the true rank for each region according to true error.
    var myRanker = new Rank(ErrorData, false);

    // Obtain Microsoft Solver Foundation's core solver object. …

我有一个符号数组,可以表示为:

from sympy import lambdify, Matrix

g_sympy = Matrix([[   x,  2*x,  3*x,  4*x,  5*x,  6*x,  7*x,  8*x,   9*x,  10*x],
                  [x**2, x**3, x**4, x**5, x**6, x**7, x**8, x**9, x**10, x**11]])

g = lambdify( (x), g_sympy )

所以每个x我得到一个不同的矩阵:

g(1.) # matrix([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.],
      #         [  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.]])
g(2.) # matrix([[  2.00e+00,   4.00e+00,   6.00e+00,   8.00e+00,   1.00e+01, 1.20e+01,   1.40e+01,   1.60e+01,   1.80e+01,   2.00e+01],
      #         [  4.00e+00,   8.00e+00, …

好吧我知道之前有一个有限的例子已经被要求缩放[-1, 1]间隔[a, b] 高斯 - 勒让德正交在numpy中的不同间隔但是没有人发布如何概括它[-a, Infinity](如下所述,但不是(还)快).这也显示了如何使用多个实现调用复杂函数(无论如何在定量选项定价中).有基准quad代码,后面跟着leggauss代码示例链接到如何实现自适应算法.我已经解决了大多数链接的adaptive algorithm困难 - 它目前打印分割积分的总和以显示它正常工作.在这里,你会发现功能的范围从转换[-1, 1]到[0, 1]以[a, Infinity](感谢@AlexisClarembeau).要使用自适应算法,我不得不创建另一个函数来转换[-1, 1]到[a, b]被反馈到[a, Infinity]功能.

import numpy as np
from scipy.stats import norm, lognorm
from scipy.integrate import quad

a = 0
degrees = 50

flag=-1.0000
F = 1.2075
K = 0.1251
vol = 0.43
T2 = 0.0411
T1 = 0.0047

def integrand(x, flag, F, K, …

我试图绘制傅里叶积分,但我在积分时得到错误

X <- seq(-10, 10, by = 0.05)
f_fourier <- function(X) {
    Y <- sapply(X, function(x) {
        integrand <- function(l) {
            y <- (2 / pi) * cos(l * x) / (l^2 + 1)
        }
        integrate(integrand, lower = 0, upper = Inf)$value
    })
}
plot(X,f_fourier(X))

错误:

maximum number of subdivisions reached

我发现"cos(l*x)"会导致这个错误,但Wolfram给了我正常的结果.你能提出什么建议吗？

我正在尝试开发物理模拟,我想实现一个四阶辛积分方法.问题是我必须弄错数学,因为在使用辛积分器时我的模拟根本不起作用(与模拟工作相当好的四阶Runge-Kutta积分器相比).我一直在谷歌搜索这个,我能找到的只是关于这个主题的科学文章.我试图改编文章中使用的方法,但我没有运气.我想知道是否有人有使用辛积分器的模拟的源代码,最好是模拟引力场,但任何辛积分器都可以.源代码的语言并不重要,但我会欣赏使用C风格语法的语言.谢谢!

我试图在MATLAB 2017a中集成一个常量函数,但我被卡住了.首先,当我使用以下脚本进行集成时,我得到了正确的输出.所以脚本适用于x0依赖的脚本t.

function E=sol(n,k)
x0 = @(t)  t^(2);
j = 0;
E = zeros(n,1);
 while j < n+1 ;
    K = matlabFunction(subs(po(j,k))) ; 
    eval(sprintf('x%d = integral(K,0,1);',j+1)) ;
    E(j+1,1) = subs(sprintf('x%d',j+1))
    j = j+1;
 end
end

功能po(j,k)如下,

function A_j = po(j,k)          % Adomian polynomials
 if j >0
  x =  sym('x',[1 j]);
    syms p;                            % Assinging a symbolic variable for p
    syms x0;
    S = x0+ sum(p.^(1:j) .* x) ;     % Sum of p*x up to …

如果有人可以帮助解决以下问题,我将不胜感激.我有以下ODE:

dr/dt = 4*exp(0.8*t) - 0.5*r   ,r(0)=2, t[0,1]       (1)

我用两种不同的方式解决了(1).通过Runge-Kutta方法(第4顺序)和ode45Matlab中的方法.我将这两个结果与分析解决方案进行了比较,分析解决方案由下式给出:

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)

当我根据确切的解决方案绘制每个方法的绝对误差时,我得到以下结果:

对于RK方法,我的代码是:

h=1/50;                                            
x = 0:h:1;                                        
y = zeros(1,length(x)); 
y(1) = 2;    
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;                   
for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end

并为ode45:

tspan = 0:1/50:1;
x0 = 2;
f = @(t,r) 4.*exp(0.8*t) - 0.5*r;
[tid, y_ode45] …

我想以数字方式整合以下内容:

哪里

和a,b并且?是常量,为简单起见,都可以设置为1.

Matlab使用dblquad,Mathematica使用NIntegrate都不能处理分母创建的奇点.由于它是双积分,我无法指定Mathematica中奇点的位置.

我确信它不是无限的,因为这个积分是基于扰动理论而没有

之前已被发现(只是不是我,所以我不知道它是如何完成的).

有任何想法吗？

我编写了一个函数来计算函数的拉普拉斯变换scipy.integrate.quad。它不是一个非常复杂的函数，目前在 Erlang 分布的概率密度函数上表现不佳。

我在下面列出了我的所有工作。我首先计算拉普拉斯变换，然后计算逆变换，以便将其与 Erlang 的原始 pdf 进行比较。我用mpmath这个。这mpmath.invertlaplace不是问题，因为它成功地将封闭式拉普拉斯变换非常完美地转换回原始 pdf。

请帮助我理解我的数值拉普拉斯变换的问题是什么。我收到以下错误，但无法解决它。

IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. a=0,b=np.inf,limit=limit)[0]

阴谋

代码

import numpy as np
import matplotlib.pyplot as plt
import math as m
import mpmath as mp
import scipy.stats as st
from scipy.integrate import quad


def get_laplace(func,limit=10000):
    
    '''
    Returns laplace transfrom function
    '''
    
    def laplace(s):
        '''Numerical laplace transform'''
        # Seperate into …