willmwade's profile - activity

The only other item would be a mess I keep running into at defining the variable of integration. integral(f(x),x,0,1) verses numerical_integral(f(x),0,1). Numerical integral errors when the variable of integration is declared. If it would not do that, it would help bring the two into uniformity.

For example in the answer to this one. I did not know that f=lambda x: x^2 would allow me to do f(2) result 4. Nor has any of the python docs on it been the most help. Mind I do more in Java :( and PHP :) than the little I have done in Python, but still some Sage documentation pointing at others?

The main area I think that could use some documentation is the use of the Python inline lambda function. Most if not all of the integration issues I have had, have been solved with using this. However there is little in how to use it with Sage specifically.

This worked great!

This works for me:

sage: var('d')
d
sage: f = lambda x,y: numerical_integral(1/d*2*x*y,.01,Infinity)[0]
sage: f(3,1)
406.69135669845832

You could try with numerical_integral. You didn't provide definition of B(), so here is example with simpler function:

sage: var('d h')
  (d, h)
sage:  integral(integral(exp(-(h-1)^2-(d-2)^2), d, 0,Infinity), h, 0, Infinity)
  1/4*pi + 1/4*(pi + pi*erf(1))*erf(2) + 1/4*pi*erf(1)
sage: n(_)
  2.88773776713433
sage: numerical_integral(lambda h: numerical_integral(lambda d:exp(-(h-1)^2-(d-2)^2), 0,Infinity)[0], 0, Infinity)
  (2.8877377671374997, 2.8275945874944256e-06)

BTW: I'd like to know how to do the above without 'lambda'.

BTW2: Are you sure your integral converges at d=0?

f(x,y)=numerical_integral(1/d*2*x*y,.01,Infinity)[0]
Error

So what I want is to the integration wait until after the variables have been substituted so that it is able to numerically integrate. (Yes I need to numerically integrate. This is a simplified form that reproduces the same result.)

f(3,1)=numerical_integral(1/d*2*3*1,.01,Infinity)[0]
406.69135669845798

I thought there might be a way using a lambda defined function, but I was unable to find one.

I found out how. Deleting ~/.sage did it. rm -Rf ~/.sage

Hmm... I still get the same error with your options. Perhaps I broke something while trying to fix the problem. Any way to reset everything? Note I tried reset() and it still produced the error.

integral(integral(e^(-0.00260657639223762*(h - 11.1600000000000)^2 - 1.34372480515*(d - 2.85000000000000)^2)/pi, d, 0,Infinity), h, 0, Infinity).n()

I have been working through an error I am receiving with the above command. The error I get when I run it is:

TypeError: Error executing code in Maxima
CODE:
    sage297 : integrate(sage293,sage294,sage295,sage296)$
Maxima ERROR:

rat: replaced -0.0026065763922 by -200/76729 = -.00260657639223762

rat: replaced -11.16 by -279/25 = -11.16

rat: replaced -1.343724 by -3233/2406 = -1.34372402327515

rat: replaced -2.85 by -57/20 = -2.85
Maxima encountered a Lisp error:

 5883790627 is not of type FIXNUM.

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.

However if I shorten the decimals to something like:

integral(integral(e^(-0.0026*(h - 11.1600000000000)^2 - 1.343724*(d - 2.85)^2)/pi, d, 0,Infinity), h, 0, Infinity).n()

Then the command runs fine without any issues. Now perhaps this is a bug, but really I do not need the level of precision that is in these commands. Of course I do not want numbers here but variables and it is evaluating the variables to this level of precision that is causing the issue. For example I say:

sage: f(d,h) = (1/(2*pi*0.61*13.85))*exp(-1/2*((d-2.85)^2/0.61^2+(h-13.85)^2/13.85^2))
sage: f
(h, d) |--> 0.05918210333195241758892111025625850742735*e^(-0.002606576392237615503916381029337017294634*(h - 11.16000000000000000000000000000000000000)^2 - 1.343724805159903251814028486965869389949*(d - 2.850000000000000000000000000000000000000)^2)/pi

Any ideas? P.S. I have searched around, but I have not found any that apply to just basic numbers in Sage.

So I have equations:

B(x,d,h,G)=(a01+b01*x+c01*x^2+d01*x^3) + (a02+b02*x+c02*x^2+d02*x^3)*h + (a11+b11*x+c11*x^2+d11*x^3)*(G/d) +  (a12+b12*x+c12*x^2+d12*x^3)*(G/d)*h +  (a21+b21*x+c21*x^2+d21*x^3)*(G/d)^2 +  (a22+b22*x+c22*x^2+d22*x^3)*((G/d)^2)*h

normal(x,av,sd)=(1/(sd*sqrt(2*pi)))*exp(-(x-av)^2/(2*sd^2))

brk(x,d,h,D,H,Dstd,Hstd,G)=B(x,d,h,G)*normal(d,D,Dstd)*normal(h,H,Hstd)

(a##-d## are all decimals)

And I want to double integrate over d and h. So just for example I can integrate over just d:

f(x) = break1(x,d,11.16,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity)
f(500).numerical_approx()                                              
1.04299933381999

and it evaluates just fine. Likewise if I put in an value for d and integrate over h, it also produces a value.

So I want the numerical approximation of this double integral, but when I try for it using:

f(x) = break1(x,d,h,2.85,11.16,0.61,13.85,.5).integrate(d,0,infinity).integrate(h,0,infinity)
f(500).numerical_approx()

/home/wil/sage/local/lib/python2.6/site-packages/sage/symbolic/integration/integral.pyc in _evalf_(self, f, x, a, b, parent)
    199         # The gsl routine returns a tuple, which also contains the error.
    200         # We only return the result.
--> 201         return numerical_integral(f, a, b)[0]
    202 
    203     def _tderivative_(self, f, x, a, b, diff_param=None):

/home/wil/sage/local/lib/python2.6/site-packages/sage/gsl/integration.so in sage.gsl.integration.numerical_integral (sage/gsl/integration.c:1551)()

ValueError: Integrand has wrong number of parameters

At any rate, I think what is happening is that it is trying to evaluate the inside integral numerically first perhaps, which it is not able to do as it has a variable?

I tried using the lambda in the first answer, and I was able to evaluate, but I had to set it to max_points=10 to get an answer that was even close to correct, plus I could find no way to plot that one.

Thanks for any help! Wil

Thanks; I don't quite understand what's broken, but I think I have a workaround.

sage: normal(x,av,sd)=((1/(sd*sqrt(2*pi)))*exp(-(x-av)^2/(2*sd^2)))
sage: f(x,y)= x*y^3*normal(y,1,2)
sage: g = f.integrate(y,-2,2)
sage: g
(x, y) |--> 1/4*sqrt(2)*x*integrate(y^3*e^(-1/8*(y - 1)^2), y, -2, 2)/sqrt(pi)

(note that sage seems to think g is a function of x and y, even though it shouldn't be)

sage: h(x) = f.integrate(y,-2,2)
sage: h
x |--> 1/4*sqrt(2)*x*integrate(y^3*e^(-1/8*(y - 1)^2), y, -2, 2)/sqrt(pi)

(this is looking better, but still doesn't work . . . the problem seems to be the integrate command in the definition of h)

sage: plot(h,0,5)
...
ValueError: free variable: y

(the first way I could think of to get Sage to evaluate the integral was to use the numerical_approx method; maybe there's a better way, but this works :)

sage: k = lambda x: h(x).numerical_approx()
sage: plot(k,0,5)

(graph is shown!)

Note that, unfortunately, the following does not work:

k(x) = h(x).numerical_approx()
...
TypeError: cannot evaluate symbolic expression numerically

Unless someone knows better, I'm inclined to think this should be filed as a Trac ticket.

EDIT: after writing all this, I think the problem is the "symbolic expression" integrate(...). The following works as expected:

sage: k(x) = 1/4*sqrt(2)*x*(integrate(y^3*e^(-1/8*(y - 1)^2), y, -2,2).numerical_approx())/sqrt(pi)
sage: plot(k,0,5)

See below.

Oh and before this I did run var('x d') so both of those are declared.

I only have trouble with it when I want to integrate with the same variable as is in normal() and having that variable outside of normal(). Example: x/(d+1)*normal(d,2.85,0.61) the integral(normal(d,2.85,0.61),d,0,Infinity)= just under 1

normal(x,av,sd)=(1/(sd*sqrt(2*pi)))*exp(-(x-av)^2/(2*sd^2)) If I integrate this alone at say the standard normal (x,0,1) or any other mean and standard deviation I get a number just like you would expect. It is just the normal distribution.

I figured out what the [0] was! It was limiting the return to only the value and leaving out the error calculation.

This should be a rather simple question, and I think I am just not understanding the error.

I have a function say f(x,y)= x*y^3*normal(y,1,2) where normal() is the normal pdf equation (x,mean,standard deviation).

So I want to integrate across the pdf and then plot x. What I thought I would do was:

g = f.integral(y,0,Infinity)
plot(g,-10,10)

But that give the error: ValueError: free variable: x

What am I missing?

The full equation takes a really long time. I am guessing that it is just trying to be too precise. I don't need much. Anyway to decrease the precision?