Ask Your Question
0

To obtain multiple integral from distributions functions

asked 2023-05-22 17:12:48 +0200

Cyrille gravatar image

In the dominance stochastic problems for the density $f(x)$, we need to evaluate

$\int_{-\infty}^x f(s) ds$, $\int_{-\infty}^t\int_{-\infty}^x f(s) dsdt$, $\int_{-\infty}^u\int_{-\infty}^t\int_{-\infty}^x f(s) dsdxdu$ and so on even if, the most part of the time one stops for serious reasons. Sagemath has a powerfull set of pre programed distribution. For instance for the uniform $U[a,b]$ For instance, I can define

def uniform(x,a,b) :
    T = RealDistribution('uniform', [a, b])
    return T.distribution_function(x)

and ask

uniform(1,0,2)

to obtain $.5$. But if I ask

integral(uniform(x,0,2), (x, 0,1))

I will have an error of the type unable to simplify to float approximation. After a lot of reading about this error, it seems that I should either add somewhere hold=True or some thing else or perhaps I should define all the distribution functions from scratch ?

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
0

answered 2023-05-25 13:34:29 +0200

slelievre gravatar image

For a numerical approximation, use numerical_integral.

It also gives an estimate for the error.

You can extract the value and discard the estimate.

sage: f = lambda x: RealDistribution('uniform', [0, 2]).distribution_function(x)

sage: numerical_integral(f, 0, 1)
(0.5, 5.551115123125783e-15)

sage: numerical_integral(f, 0, 1)[0]
0.5
edit flag offensive delete link more
0

answered 2023-05-25 08:38:36 +0200

Emmanuel Charpentier gravatar image

First (but incidental): by naming your function 'uniform', you scratch the predefined uniform function, which may be of some usefulness in your problem :

Signature:      uniform(a, b)
Docstring:     
   Get a random number in the range [a, b).

   Equivalent to a + (b-a) * random().

   EXAMPLES:

      sage: s = uniform(0, 1); s  # random
      0.111439293741037
      sage: 0.0 <= s <= 1.0
      True

      sage: s = uniform(e, pi); s  # random
      0.5143475134191677*pi + 0.48565248658083227*e
      sage: bool(e <= s <= pi)
      True

But your problem isn't here. consider :

sage: def unif(x,a,b) :
....:     T = RealDistribution('uniform', [a, b])
....:     return T.distribution_function(x)
....:

That's what you did. but note that :

sage: type(unif)
<class 'function'>

Which tells you that your function is a Python (numeric) function (wrapping in Python a gsl library C function), not a symbolic expression nor a symbolic function (a. k. a. callable symbolic expression).

This function can handle numerical arguments fine :

sage: unif(0.5,0,2)
0.5

but fails to handle symbolic arguments :

sage: unif(x,0,2)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
File /usr/local/sage-10/src/sage/symbolic/expression.pyx:1962, in sage.symbolic.expression.Expression.__float__()
   1961 try:
-> 1962     ret = float(self._eval_self(float))
   1963 except TypeError:

File /usr/local/sage-10/src/sage/symbolic/expression.pyx:1644, in sage.symbolic.expression.Expression._eval_self()
   1643 else:
-> 1644     raise TypeError("cannot evaluate symbolic expression to a numeric value")
   1645 

TypeError: cannot evaluate symbolic expression to a numeric value

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
File /usr/local/sage-10/src/sage/symbolic/expression.pyx:1965, in sage.symbolic.expression.Expression.__float__()
   1964 try:
-> 1965     c = (self._eval_self(complex))
   1966     if imag(c) == 0:

File /usr/local/sage-10/src/sage/symbolic/expression.pyx:1644, in sage.symbolic.expression.Expression._eval_self()
   1643 else:
-> 1644     raise TypeError("cannot evaluate symbolic expression to a numeric value")
   1645 

TypeError: cannot evaluate symbolic expression to a numeric value

Here is the source of your problem, plainly stated.

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In [50], line 1
----> 1 unif(x,Integer(0),Integer(2))

Cell In [46], line 3, in unif(x, a, b)
      1 def unif(x,a,b) :
      2     T = RealDistribution('uniform', [a, b])
----> 3     return T.distribution_function(x)

File /usr/local/sage-10/src/sage/probability/probability_distribution.pyx:851, in sage.probability.probability_distribution.RealDistribution.distribution_function()
    849 """
    850 if self.distribution_type == uniform:
--> 851     return sage.rings.real_double.RDF(gsl_ran_flat_pdf(x, self.parameters[0], self.parameters[1]))
    852 elif self.distribution_type == gaussian:
    853     return sage.rings.real_double.RDF(gsl_ran_gaussian_pdf(x, self.parameters[0]))

File /usr/local/sage-10/src/sage/symbolic/expression.pyx:1971, in sage.symbolic.expression.Expression.__float__()
   1969             raise
   1970     except TypeError:
-> 1971         raise TypeError("unable to simplify to float approximation")
   1972 return ret
   1973 

TypeError: unable to simplify to float approximation

Therefore, you cannot pass `unif` to `integral` whose docstring states :

   Return an indefinite or definite integral of an object "x".

   First call "x.integral()" and if that fails make an object and
   integrate it using Maxima, maple, etc, as specified by algorithm.

   For symbolic expression calls "sage.calculus.calculus.integral()" -
   see this function for available options.

Checking that unif has neither integral nor maxima, gap, etc... methods is left to the reader as a sanitary exercise.

BTW : brute force nested integrations is probably not the "right" way to tackle your problem. You should look up Convolution and Characteristic function. And think...

HTH,

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2023-05-22 17:12:48 +0200

Seen: 97 times

Last updated: May 25 '23