Is that Sage vs Maxima inconsistency known ?

[Emmanuel Charpentier]

Paedagogic exercise :
Goal : simulate from a lognormal with known mean and standard deviation (for meta-analytic purposes...).
==> Intermediate goal 1 : estimate lognormal parameters mu and sigma from published mean and variance (method of moments).
==> Intermediate goal 2 : express mean and variance in terms of mu and sigma.
==> Intermediate goal 3 : define the lognormal density (starting with the normal).
Maxima does this relatively easily. Using Sage's maxima :

charpent@asus16-ec:~$ sage -maxima
;;; Loading #P"/usr/local/sage-6.3/local/lib/ecl/cmp.fas"
;;; OPTIMIZE levels: Safety=2, Space=0, Speed=3, Debug=0
;;;
;;; End of Pass 1.
;;; OPTIMIZE levels: Safety=2, Space=0, Speed=3, Debug=0
;;;
;;; End of Pass 1.
;;; OPTIMIZE levels: Safety=2, Space=0, Speed=3, Debug=0
;;;
;;; End of Pass 1.
;;; OPTIMIZE levels: Safety=2, Space=0, Speed=3, Debug=0
;;;
;;; End of Pass 1.
;;; Loading #P"/usr/local/sage-6.3/local/lib/ecl/sb-bsd-sockets.fas"
;;; Loading #P"/usr/local/sage-6.3/local/lib/ecl/sockets.fas"
;;; Loading #P"/usr/local/sage-6.3/local/lib/ecl/defsystem.fas"
Maxima 5.34.1 http://maxima.sourceforge.net
using Lisp ECL 13.5.1
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) declare(x,real,y,real,mu,real,sigma,real);
(%o1)                                done
(%i2) assume(sigma>0);
(%o2)                             [sigma > 0]
(%i3) define(phi(x,mu,sigma),%e^-((x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*%pi)));
                                                       2
                                               (x - mu)
                                             - ---------
                                                      2
                                               2 sigma
                                           %e
(%o3)            phi(x, mu, sigma) := -----------------------
                                      sqrt(2) sqrt(%pi) sigma
(%i4) define(dln(x,mu,sigma),subst([y=log(x)],phi(y,mu,sigma))*diff(log(x),x));
                                                         2
                                            (log(x) - mu)
                                          - --------------
                                                      2
                                               2 sigma
                                        %e
(%o4)           dln(x, mu, sigma) := -------------------------
                                     sqrt(2) sqrt(%pi) sigma x
(%i5) M1:integrate(x*dln(x,mu,sigma),x,0,inf);
                                       2
                                  sigma  + 2 mu
                                  -------------
                                        2
(%o5)                           %e
(%i6) V1:factor(integrate((M1-x)^2*dln(x,mu,sigma),x,0,inf));
                                2             2
                           sigma         sigma  + 2 mu
(%o6)                   (%e       - 1) %e
(%i7)

But the same computation in Sage seems to need (needlessly) further hypotheses :

sage: var("x, y, mu, sigma")
(x, y, mu, sigma)
sage: assume(x,"real",y,"real",mu,"real",sigma,"real",sigma>0)
sage: phi(x,mu,sigma)=e^-((x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*pi)).simplify()
dln(x,mu,sigma)=phi(y,mu,sigma).subs(y=log(x))*diff(log(x),x).simplify()
sage: dln
(x, mu, sigma) |--> 1/2*sqrt(2)*e^(-1/2*(mu - log(x))^2/sigma^2)/(sqrt(pi)*sigma*x)
sage: M2=integrate(x*dln(x,mu,sigma),x,0,oo).simplify()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-11-2c2ea6d23ef8> in <module>()
----> 1 M2=integrate(x*dln(x,mu,sigma),x,Integer(0),oo).simplify()

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/misc/functional.pyc in integral(x, *args, **kwds)
    800     """
    801     if hasattr(x, 'integral'):
--> 802         return x.integral(*args, **kwds)
    803     else:
    804         from sage.symbolic.ring import SR

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.integral (build/cythonized/sage/symbolic/expression.cpp:50867)()

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.pyc in integrate(expression, v, a, b, algorithm, hold)
    710         return indefinite_integral(expression, v, hold=hold)
    711     else:
--> 712         return definite_integral(expression, v, a, b, hold=hold)
    713
    714 integral = integrate

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:9283)()

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/function.so in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:5924)()

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.pyc in _eval_(self, f, x, a, b)
    173         for integrator in self.integrators:
    174             try:
--> 175                 return integrator(*args)
    176             except NotImplementedError:
    177                 pass

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/symbolic/integration/external.pyc in maxima_integrator(expression, v, a, b)
     19         result = maxima.sr_integral(expression,v)
     20     else:
---> 21         result = maxima.sr_integral(expression, v, a, b)
     22     return result._sage_()
     23

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/interfaces/maxima_lib.pyc in sr_integral(self, *args)
    783                 raise ValueError("Integral is divergent.")
    784             elif "Is" in s: # Maxima asked for a condition
--> 785                 self._missing_assumption(s)
    786             else:
    787                 raise

/usr/local/sage-6.3/local/lib/python2.7/site-packages/sage/interfaces/maxima_lib.pyc in _missing_assumption(self, errstr)
    992              + errstr[jj+1:k] +">0)', see `assume?` for more details)\n" + errstr
    993         outstr = outstr.replace('_SAGE_VAR_','')
--> 994         raise ValueError(outstr)
    995
    996 def is_MaximaLibElement(x):

ValueError: Computation failed since Maxima requested additional constraints; using the 'assume' command before evaluation *may* help (example of legal syntax is 'assume(sigma^2+mu>0)', see `assume?` for more details)
Is sigma^2+mu positive, negative or zero?

Trying to explicitely delegate the problem to Sage's Maxima fails for seemingly the same reason. Similarly, setting Maxima's domain to "real" (via maxima.domain("real")) does not allow to solve the problem.

One notes that, conversively, setting domain to complex in Sage's Maxima does not create a similar problem. It simply returns unsimplified expressions.

Is that one known ? Should it be reported to either a new ticket or an existing one ? Skimming trac for "domain" didn't turn up anything significant...

HTH,

--
Emmanuel Charpentier

[Emmanuel Charpentier]

Six month and a few versions of Sage and Maxima later, I've checked (in a different way, see below) that the same problem still exists. Nobody has a clue about this problem ?

sage: reset()
sage: var("x,mu,sigma")
(x, mu, sigma)
sage: dlnorm(x,mu,sigma)=e^-((log(x)-mu)^2/(2*sigma^2))/(sigma*x*sqrt(2*pi))
sage: integrate(dlnorm(x,mu,sigma),x,0,oo).canonicalize_radical()
1
sage: m=integrate(x*dlnorm(x,mu,sigma),x,0,oo).canonicalize_radical()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/all_cmdline.pyc in <module>()
----> 1 m=integrate(x*dlnorm(x,mu,sigma),x,Integer(0),oo).canonicalize_radical()

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/misc/functional.pyc in integral(x, *args, **kwds)
    661     """
    662     if hasattr(x, 'integral'):
--> 663         return x.integral(*args, **kwds)
    664     else:
    665         from sage.symbolic.ring import SR

/usr/local/sage-6.7/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.integral (build/cythonized/sage/symbolic/expression.cpp:52709)()
  10695                     R = ring.SR
  10696             return R(integral(f, v, a, b, **kwds))
> 10697         return integral(self, *args, **kwds)
  10698
  10699     integrate = integral

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.pyc in integrate(expression, v, a, b, algorithm, hold)
    759         return indefinite_integral(expression, v, hold=hold)
    760     else:
--> 761         return definite_integral(expression, v, a, b, hold=hold)
    762
    763 integral = integrate

/usr/local/sage-6.7/src/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:10862)()
    992             res = self._evalf_try_(*args)
    993             if res is None:
--> 994                 res = super(BuiltinFunction, self).__call__(
    995                         *args, coerce=coerce, hold=hold)
    996

/usr/local/sage-6.7/src/sage/symbolic/function.pyx in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6798)()
    500             for i from 0 <= i < len(args):
    501                 vec.push_back((<Expression>args[i])._gobj)
--> 502             res = g_function_evalv(self._serial, vec, hold)
    503         elif self._nargs == 1:
    504             res = g_function_eval1(self._serial,

/usr/local/sage-6.7/src/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction._evalf_or_eval_ (build/cythonized/sage/symbolic/function.cpp:11519)()
   1063         res = self._evalf_try_(*args)
   1064         if res is None:
-> 1065             return self._eval0_(*args)
   1066         else:
   1067             return res

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.pyc in _eval_(self, f, x, a, b)
    174         for integrator in self.integrators:
    175             try:
--> 176                 return integrator(*args)
    177             except NotImplementedError:
    178                 pass

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/symbolic/integration/external.pyc in maxima_integrator(expression, v, a, b)

     19         result = maxima.sr_integral(expression,v)
     20     else:
---> 21         result = maxima.sr_integral(expression, v, a, b)
     22     return result._sage_()
     23

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/interfaces/maxima_lib.pyc in sr_integral(self, *args)
    782                 raise ValueError("Integral is divergent.")
    783             elif "Is" in s: # Maxima asked for a condition
--> 784                 self._missing_assumption(s)
    785             else:
    786                 raise

/usr/local/sage-6.7/local/lib/python2.7/site-packages/sage/interfaces/maxima_lib.pyc in _missing_assumption(self, errstr)
    991              + errstr[jj+1:k] +">0)', see `assume?` for more details)\n" + errstr
    992         outstr = outstr.replace('_SAGE_VAR_','')
--> 993         raise ValueError(outstr)
    994
    995 def is_MaximaLibElement(x):

ValueError: Computation failed since Maxima requested additional constraints; using the '

assume' command before evaluation *may* help (example of legal syntax is 'assume(sigma^2+mu>0)', see `assume?` for more details)

Is sigma^2+mu positive, negative or zero?

sage: maxima("declare(x,real,mu,real,sigma,real);")
done
sage: maxima("assume(sigma>0);")
[sigma>0]
sage: maxima("define(dlnorm(x,mu,sigma),%e^-((log(x)-mu)^2/(2*sigma^2))/(sigma*x*sqrt(2*%pi)));")
dlnorm(x,mu,sigma):=%e^-((log(x)-mu)^2/(2*sigma^2))/(sqrt(2*%pi)*sigma*x)
sage: maxima("radcan(integrate(dlnorm(x,mu,sigma),x,0,inf));")
1
sage: maxima("m:radcan(integrate(x*dlnorm(x,mu,sigma),x,0,inf));")
%e^((sigma^2+2*mu)/2)
sage: maxima("v:radcan(integrate((x-m)^2*dlnorm(x,mu,sigma),x,0,inf));")
%e^(2*sigma^2+2*mu)-%e^(sigma^2+2*mu)
sage:

HTH,
--
Emmanuel Charpentier

[Nils Bruin]

On Monday, May 4, 2015 at 1:42:21 PM UTC-7, Emmanuel Charpentier wrote:

Six month and a few versions of Sage and Maxima later, I've checked (in a different way, see below) that the same problem still exists. Nobody has a clue about this problem ?

Well, at least if you pass to maxima exactly the result that sage ends up passing to the maxima integrator, we get a useless result from maxima as well:

 (%i8) F: 1/2*sqrt(2)*e^(-1/2*(mu - log(x))^2/sigma^2)/(sqrt(%pi)*sigma*x);
(%o8) 1/(sqrt(2)*sqrt(%pi)*e^((mu-log(x))^2/(2*sigma^2))*sigma*x)
(%i9) integrate(x*F,x,0,inf);
(%o9) ('integrate(1/e^((mu-log(x))^2/(2*sigma^2)),x,0,inf)) /(sqrt(2)*sqrt(%pi)*sigma)

it seems the maxima that sage uses gets prepped with additional methods (which then fail on some assumption) but it does not seem to be any of

display2d : false;
domain : complex;
keepfloat : true;
load(to_poly_solve);
load(simplify_sum);
load(abs_integrate);
load(diag);

which I thought was the preamble we provide. Perhaps something else gets added in some other dark corner?