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 sage-support+unsubscr...@googlegroups.com. To post to this group, send email to sage-support@googlegroups.com. Visit this group at http://groups.google.com/group/sage-support. For more options, visit https://groups.google.com/d/optout.
