当前位置：文江博客话题详情

wolfram-mathematica maxima

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

发布于 2024-10-16 21:38:42 字数 2331 浏览 2 评论 0 原文

我正在尝试最大化我的 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 w/某些值崩溃：

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 的开源程序可以数值上近似这个函数？我真的很想释放开放源代码到开源平台。

原文

I'm trying to Maxima-fy my Mathematica box options formula
(https://github.com/barrycarter/bcapps/blob/master/box-option-value.m)
but Maxima crashes on a fairly simple integration:

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));

Evaluating upandin w/ certain values crashes:

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.

Without the float() in upandin, Maxima just leaves the integral in
original form.

Can someone help? I thought converting Mathematica to Maxima would be
easy, but now I'm not as sure.

The Mathematica version works fine:

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

EDIT: Is there any open source Mathematica-like program that WILL
numerically approximate this function? I'd really like to release open
source code to an open source platform.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

久而酒知 2024-10-23 21:38:42

（我可能没有权利回答这个问题，但是......）

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

可能需要检查

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

或

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

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

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

Daniel Lichtblau
沃尔夫勒姆研究公司

回复收藏 0 原文

吹梦到西洲 2024-10-23 21:38:42

使用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]

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

Use quad_qagi to numerically approximate an integral over an infinite interval. ?? quad_ shows info about Quadpack functions.

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]

Sorry for the late reply. Leaving this here in case someone finds it by searching.

回复收藏 0 原文

猫九 2024-10-23 21:38:42

我知道 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 的四元组。 [我在这里猜测了“正确”值，但由于我们一开始就没有那么精确，所以这有点愚蠢。]

I know Maxima tries very hard to avoid floats, and I think that's what it's trying to do here, but I'm not enough of a Maxima guru to explain how to prevent it. Pretty much anything numerical can handle this, although you might have to break the interval or transform the integrand manually. Note that you say it's fairly simple, but it's awfully steep: for these parameters, the integrand is ~6*10^(-34) at 0.1 and ~3*10^(-206) at -0.1. That's enough of a range to give lots of naive integration algorithms fits.

Anyway, you can do it easily enough in Sage using tools from scipy and gsl behind the scenes:

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)

or use mpmath's quad if you need arbitrary precision. [I had a guess at the "correct" value here, but since we don't have that much precision to start with, it's kind of silly.]

回复收藏 0 原文

梦途 2024-10-23 21:38:42

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

我认为您正在寻找的是数字集成例程 - Maxima 提供了各种集成例程，从非常基本的 romberg 到各种 Quadpack 方法（尝试使用 ??quad 获取文档）。

-s

PS 至于“这个蹩脚的‘开源’东西”——是什么造成的？您可能想查看 Wikipedia 文章中 Macsyma/Maxima 的历史以获得一些观点。

Maxima's 'integrate' function does symbolic, not numeric, integration. When it returns a noun form from an integral, that means it can't perform the (symbolic) integration. Changing the parameters of the expression from exact to floating (using 'float') won't change that.

I think what you're looking for is a numeric integration routine -- Maxima offers a variety of those, from the very basic romberg to a variety of Quadpack methods (try ?? quad for documentation).

-s

PS As for "this crappy 'open source' stuff" -- what brought that on? You might want to look at the history of Macsyma/Maxima in the Wikipedia article for some perspective.

回复收藏 0 原文

~没有更多了~