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
Le samedi 25 octobre 2014 11:28:47 UTC+2, Emmanuel Charpentier a écrit :
>
> 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
>
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
