Maxima 在相对简单的积分上崩溃答案

【问题标题】：Maxima crashes on relatively simple integralMaxima 在相对简单的积分上崩溃
【发布时间】：2011-02-08 03:02:38
【问题描述】：

我正在尝试最大化我的 Mathematica 框选项公式（https://github.com/barrycarter/bcapps/blob/master/box-option-value.m）但是 Maxima 在一个相当简单的集成上崩溃了：

load(distrib); 
pdflp(x, p0, v, p1, p2, t1, t2) := pdf_normal(x,log(p0),sqrt(t1)*v); 
cdfmaxlp(x, p0, v, p1, p2, t1, t2) := 1-erf(x/(v*sqrt(t2-t1)/sqrt(2))); 

upandin(p0, v, p1, p2, t1, t2) :=  
 integrate( 
 float( 
 pdflp(x, p0, v, p1, p2, t1, t2)* 
 cdfmaxlp(log(p1)-x, p0, v, p1, p2, t1, t2) 
 ), 
 x, minf, log(p1));

使用某些值评估 upandin 崩溃：

upandin(1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425); 

rat: replaced -.00995033085316809 by -603/60601 = -.00995033085262619 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced 8116.5 by 16233/2 = 8116.5 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced 1.0 by 1/1 = 1.0 

rat: replaced 1.792882852833688 by 4484/2501 = 1.792882846861255 

rat: replaced 180.1832400641081 by 126849/704 = 180.1832386363636 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced -1.0 by -1/1 = -1.0 

rat: replaced 1.792882852833688 by 4484/2501 = 1.792882846861255 

rat: replaced 180.1832400641081 by 126849/704 = 180.1832386363636 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced 1.0 by 1/1 = 1.0 

rat: replaced -1.0 by -1/1 = -1.0 
Maxima encountered a Lisp error: 

 The value 16090668801 is not of type FIXNUM.

没有 upandin 中的 float()，Maxima 只是将积分留在原始形式。

有人可以帮忙吗？我认为将 Mathematica 转换为 Maxima 将是容易，但现在我不太确定。

Mathematica 版本运行良好：

pdflp[x_, p0_, v_, p1_, p2_, t1_, t2_] :=  
 PDF[NormalDistribution[Log[p0],Sqrt[t1]*v]][x] 

cdfmaxlp[x_, p0_, v_, p1_, p2_, t1_, t2_] := 1-Erf[x/(v*Sqrt[t2-t1]/Sqrt[2])]; 

(* NIntegrate below "equivalent" to Maximas float(); no closed form *) 

upandin[p0_, v_, p1_, p2_, t1_, t2_] :=  
 NIntegrate[pdflp[x, p0, v, p1, p2, t1, t2]* 
           cdfmaxlp[Log[p1]-x, p0, v, p1, p2, t1, t2], 
{x, -Infinity, Log[p1]}] 

upandin[1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425] 

0.0998337

编辑：是否有任何类似 Mathematica 的开源程序可以数值近似这个函数？我真的很想公开源代码到开源平台。

【问题讨论】：

在我的机器上，当我评估它说“太大而无法表示为 DOUBLE-FLOAT:”，然后给出一个 9774 位整数。我实际上并不熟悉 Maxima，但也许您可以在早期将所有内容都转换为浮点数？
奇怪的是“rat：”消息暗示浮点数正在被转换回有理数，这可能会使事情变得更糟。而且千里马没有任意精度数？！
也许你可以把问题发到ask.sagemath.org，因为sage 中的很多微积分仍然使用Maxima。

标签： wolfram-mathematica maxima

【解决方案1】：

（我可能没有资格回答这个问题，但是……）

只是一个猜测，但似乎集成想要再次使输入精确，并且可能正在执行一些涉及有理算术的困难的 bignum 计算。它使您的近似 e（欧拉数）合理化，这意味着它的行为可能与带有精确输入的积分（0）不同。

可能要检查

http://eagle.cs.kent.edu/MAXIMA/maxima_21.html

或

http://www.delorie.com/gnu/docs/maxima/maxima_62.html

用于专用数字代码，例如来自 Quadpack。

（仍然想知道我为什么要回答这个问题。Stack Overflow 上一定有 Maxima 的专业知识。）

丹尼尔·利希特布劳 Wolfram 研究

【讨论】：

你还签了公司的名字？嘘！正确的答案应该是“看，这就是为什么你需要购买 Mathematica 而不是在这些蹩脚的‘开源’东西上闲逛”。
@Timo 我与 Maxima 没有任何争执，无论是程序还是实现者。此外，我认为领导这项工作的 Robert Dodier 有时会问或回答发布在 MathGroup 上的问题。总是友好和尊重，所以我有动力不去惹他。只希望他回答你的问题。综上所述，您说得对，我在公司的名字上签名可能很愚蠢。应该补充说我不是发言人......
@Barry 抱歉，我指的是你，而不是 Timo（仍在考虑更早的回复，不同的问题）。
您回答问题的努力应该得到更多尊重。继续前进！

【解决方案2】：

使用 quad_qagi 对无限区间上的积分进行数值近似。 ?? quad_ 显示有关 Quadpack 函数的信息。

load (distrib);
pdflp (x, p0, v, p1, p2, t1, t2) := pdf_normal (x, log(p0), sqrt(t1)*v); 
cdfmaxlp (x, p0, v, p1, p2, t1, t2) := 1 - erf(x/(v * sqrt(t2 - t1)/sqrt(2))); 

upandin (p0, v, p1, p2, t1, t2) := block ([integrand],
   integrand : pdflp (x, p0, v, p1, p2, t1, t2) * cdfmaxlp (log(p1) - x, p0, v, p1, p2, t1, t2),
   quad_qagi (integrand, x, minf, log(p1))); 

upandin (1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425);
 => [.09983372557898755, 2.839204848435967E-10, 225, 0]

抱歉回复晚了。将其留在这里以防有人通过搜索找到它。

【讨论】：

【解决方案3】：

我知道 Maxima 非常努力地避免浮动，我认为这就是它在这里试图做的，但我还不足以作为 Maxima 专家来解释如何防止它。几乎任何数字都可以处理这个问题，尽管您可能必须手动打破间隔或转换被积函数。请注意，您说它相当简单，但非常陡峭：对于这些参数，被积函数在 0.1 时为 ~6*10^(-34)，在 -0.1 时为 ~3*10^(-206)。这足以让许多简单的集成算法适合。

无论如何，您可以在 Sage 中使用来自 scipy 和 gsl 的工具在幕后轻松完成：

import scipy.stats

def pdflp(x,p0,v,t1):
    return scipy.stats.norm(log(p0), sqrt(t1)*v).pdf(x)

def cdfmaxlp(x,v,t1,t2):
    return (1-erf(x/(v*sqrt(t2-t1)/sqrt(2.))))

def upandin(p0, v, p1, p2, t1, t2):
    integrand = lambda x: (pdflp(x,p0,v,t1) * cdfmaxlp(log(p1)-x,v,t1,t2))
    return numerical_integral(integrand, -Infinity, log(p1))

sage: upandin(1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425)
(0.099833725578983457, 7.5174412058308382e-07)

如果您需要任意精度，也可以使用 mpmath 的四边形。 [我对这里的“正确”值进行了猜测，但由于我们一开始没有那么精确，所以有点傻。]

【讨论】：

【解决方案4】：

Maxima 的“积分”函数进行符号积分，而不是数字积分。当它从积分返回名词形式时，这意味着它不能执行（符号）积分。将表达式的参数从精确更改为浮动（使用 'float'）不会改变这一点。

我认为您正在寻找的是一个 numeric 集成例程——Maxima 提供了多种集成例程，从非常基本的 romberg 到各种 Quadpack 方法（尝试 ?? quad 获取文档)。

-s

PS 至于“这种糟糕的‘开源’东西”——是什么造成的？您可能想从 Wikipedia 文章中查看 Macsyma/Maxima 的历史以了解一些观点。

【讨论】：

我相当肯定“蹩脚的‘开源’东西”只是对公司代表的期望的嘲弄，相比之下，是对 Daniel 方法的赞扬。