Максима вылетает на относительно простом интеграле

Question

Я пытаюсь Maxima-fy мою формулу опций Mathematica box (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. 

Без функции float () в upandin 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, которая будет численно аппроксимировать эту функцию?на платформу с открытым исходным кодом.

Answer 1

(у меня, вероятно, нет никакого дела, отвечая на это, но ...)

Просто предположение, но кажется, что интеграция хочет снова сделать ввод точным, и, возможно, выполняет некоторые сложные вычисления Бигнума с рациональнымарифметика.Он рационализирует ваше приблизительное e (число Эйлера) так, чтобы оно могло вести себя иначе, чем интегрирование (0 с точным вводом.

Возможно, вы захотите проверить

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

или

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

для выделенного числового кода, например, из Quadpack.

(Все еще задаюсь вопросом, почему я даже пытаюсь ответить на этот вопрос. Должен быть опыт Maxima где-то в стекеПереполнение.)

Даниэль Лихтблау Вольфрам Исследования

Answer 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]

Извините за поздний ответ. Оставьте это здесь на случай, если кто-то найдет его в поиске.

Answer 3

Я знаю, что Максима очень старается избегать поплавков, и я думаю, что это то, что она пытается здесь сделать, но мне не хватает гуру Максима, чтобы объяснить, как это предотвратить. Практически все числовые могут справиться с этим, хотя вам, возможно, придется разбить интервал или преобразовать подынтегральное выражение вручную. Обратите внимание, что вы говорите, что это довольно просто, но очень круто: для этих параметров подынтегральное выражение составляет ~ 6 * 10 ^ (- 34) при 0,1 и ~ 3 * 10 ^ (- 206) при -0,1. Этого достаточно, чтобы дать множество наивных алгоритмов интеграции.

В любом случае, вы можете сделать это достаточно легко в 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, если вам нужна произвольная точность. [Я догадывался о «правильном» значении здесь, но, поскольку у нас не так много точности для начала, это немного глупо.]

Answer 4

Функция Maxima 'интегрировать' выполняет символическую, а не числовую интеграцию. Когда он возвращает существительную форму из интеграла, это означает, что он не может выполнить (символическую) интеграцию. Изменение параметров выражения с точного на плавающее (используя 'float') не изменит этого.

Я думаю, что вы ищете числовую процедуру интеграции - Maxima предлагает различные из них, от базового romberg до множества методов Quadpack (попробуйте quad для документации) .

    -s

PS Что касается "этого дерьмового 'открытого' 'материала" - что вызвало это? Возможно, вы захотите взглянуть на историю Macsyma / Maxima в статье в Википедии для некоторой перспективы.