[sage-devel] Change to integration of piecewise function (consensus needed)
I've submitted a patch that changes the default behaviour of integration on piecewise functions. The patch allows indefinite integration of piecewise functions, and makes it the default. Given that this changes current behaviour, consensus is needed before the patch is accepted, which is why I bring it up here. Current behaviour: sage: pw = Piecewise([[(0,1), x*2], [(1,2), x + 3]]) sage: pw.integral() 11/2 New behaviour: sage: pw = Piecewise([[(0,1), x*2], [(1,2), x + 3]]) sage: pw.integral() Piecewise defined function with 2 parts, [[(0, 1), x^2], [(1, 2), (x^2 + 6*x)/2 - 5/2]] sage: pw.integral(definite=True) 11/2 Justification: - Consistency within Sage. For most functions, f.integrate() will return the indefinite integral by default, so it is inconsistent that piecewise functions currently return the definite integral by default. - Consistency with other math software. Maple and Mathematica (with the online integrator tool) both return indefinite integrals of piecewise functions by default. - Indefinite integrals of piecewise functions are useful (for example, in dealing with continuous probability distributions), and this seems to me to be the most intuitive way to add them to Sage. I won't go into the details of the implementation because this has already been discussed in the trac thread and sage-devel thread. Trac ticket #4721 http://trac.sagemath.org/sage_trac/ticket/4721 Previous sage-devel thread: http://groups.google.com/group/sage-devel/browse_thread/thread/f6cb26796d39c67/c19f18e80a65164e --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to sage-devel-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: Integral of piecewise functions
David, thanks.
What do you mean by an+ and an- in the notation?
From what I can tell, that does look like the same algorithm I plan to
use. I have some working code now, I just want to review it and make
sure I didn't miss any edge cases before re-posting the patch. The
code posted as a patch now uses the same algorithm, but doesn't
properly handle unbounded functions.
-- Paul
On Tue, Dec 9, 2008 at 8:32 AM, David Joyner [EMAIL PROTECTED] wrote:
Paul:
It was not at all clear to me when you first started working on the
antiderivatives of piecewise defined functions that it was possible or
made sense.
Now I have changed my mind to an extent. If you think of the
antiderivative as the
(multi-valued) inverse of the differentiation map, as I usually do, then I'm
not sure what to do. But if you just want a function whose derivavtive is the
original and which satisfies the FTC, then I think I have an algorithm which
might work.
Let f(x) be the piecewise defined function which is fn(x) for anxa(n+1),
for n in ZZ (and let us ignore endpots for now). Here a_n is an
infinite sequence
a_{-infty}=-infty and a_{infty}=+infty and each f_n(x) is Riemann integrable
with antiderivative Fn(x)+Cn. The question is can we piece the Cns's
together so that
F(x) = Fn(x)+Cn for anxa(n+1) satisfies the FTC? These constants must
satisfy
Cn = F_{n-1}(an-)-Fn(an+)+C_{n-1},
If you fix C0 = C to be an arbitrary constant (you could take this to
be 0 for example),
this is a 1-step recursionrelation which can be solved.
I think this should do the job.
Thanks for raising this issue!
On Sat, Dec 6, 2008 at 5:06 PM, Paul Butler [EMAIL PROTECTED] wrote:
When a1 = -infinity, I would make F1 = integrate(f, x, a2, x) instead of
integrate(f, x, a1, x). Then I would not calculate the definite integral of
the first interval, which would align my constants so that F(a2) = 0. When I
get a chance, I'll add this to my code.
Functions like floor with an infinite number of pieces are beyond the scope
of the Piecewise class, because piecewise functions in Sage can only (at the
moment) have finitely many pieces. I can't give an algorithm for integration
of all piecewise functions with multiple pieces off the top of my head, but
I'll give it some thought.
-- Paul
On Sat, Dec 6, 2008 at 4:13 PM, David Joyner [EMAIL PROTECTED] wrote:
Okay, this helps me understand what you mean.
Still, the case a1=-ifinty is precisely the special case which I don't
understand.
For example, take a function such as f(x) = max(1,floor(x)), x real.
How do you define an antiderivative F(x) so that
F(b)-F(a) = area under the y=f(x) for axb?
(And mayeb you can do it for that special function,
and let us ignore points of discontinuity for the sake of
discussion.) In other words, I am asking for the algorithmic procedure
you would use to create an area function of a piecewise-defined
function on the reals.
With this definition, F(b) - F(a) can be used to find the Riemann sum
between a and b. Also, F'(x) = f(x) seems to hold, except at points
where
f(x) goes from defined to undefined or vice-versa.
The antiderivative is only well-defined up to an additive constant.
IMHO, the piecewise defined function of antiderivavtives
int f1(x) dx +C1 , a1x=a2,
int f2(x) dx +C2, a2x=a3,
...
int fn(x) dx +Cn, anx=a{n+1}
does not make sense.
I agree that it doesn't make sense where C1 .. Cn are arbitrary
constants.
-- Paul
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-devel] Integral of piecewise functions
Currently, taking the integral of a piecewise function in Sage gives you the definite integral. I've proposed on trac that the integral of piecewise functions be indefinite by default. This would be consistent with how integration works on other functions in Sage, as well as piecewise functions in Maple and Mathematica. The main concern is whether the integral of a piecewise function is even well-defined. It seems to me that at least for continuous piecewise functions, the indefinite integral is well-defined. The anti-derivative is well defined, and by the fundamental theorem of calculus, the indefinite integral of a continuous function is the anti-derivative. As for discontinuous piecewise functions, I'm finding it difficult to convince myself either way. The trac ticket is 4721 ( http://trac.sagemath.org/sage_trac/ticket/4721 ) -- Paul --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Re: Integral of piecewise functions
On Sat, Dec 6, 2008 at 12:37 PM, David Joyner [EMAIL PROTECTED] wrote:
On Sat, Dec 6, 2008 at 11:39 AM, Paul Butler [EMAIL PROTECTED] wrote:
Currently, taking the integral of a piecewise function in Sage gives
you the definite integral. I've proposed on trac that the integral of
piecewise functions be indefinite by default. This would be consistent
with how integration works on other functions in Sage, as well as
piecewise functions in Maple and Mathematica.
The main concern is whether the integral of a piecewise function is
even well-defined. It seems to me that at least for continuous
piecewise functions, the indefinite integral is well-defined. The
anti-derivative is well defined, and by the fundamental theorem of
calculus, the indefinite integral of a continuous function is the
anti-derivative. As for discontinuous piecewise functions, I'm finding
it difficult to convince myself either way.
First, the indefinite integral is the anti-derivative, by definition.
What the FTC says is that although the indefinite integral evalated
at b is not
well-defined, and the same antiderivative evalated at a is also not
well-defined,
their difference *is* well-defined. Moreover, this difference agrees with
the
definite integral defined by the Riemann sum between a and b.
Second, I think (but I am not sure), what you want when you say
indefinite intergal is not the indefinite integral but is a function who
derivavtive is the original function, defined as follows:
Actually, a function F which is simply the antiderivative of f is no use to
me, I need a function with the property that F(b) - F(a) is the Riemann sum
between a and b. (They are only guaranteed by the FTC to be the same
function if f is continuous.) Maybe there is a better term for this?
if the orginial function f(x) is
f1(x), a1x=a2,
f2(x), a2x=a3,
...
fn(x), anx=a{n+1}
(and 0 outside (a1,a{n+1}) then I guess you want to define the integral,
call it F, by
int_{a1}^x f1(t) dt, a1x=a2,
int_{a2}^x f2(t) dt, a2x=a3,
...
int_{an}^x fn(t) dt, anx=a{n+1}
Is this correct? This is not the antiderivative
but it does have the property that F'(x)=f(x).
I think you and I are defining antiderivative differently. I'm using the
definition that F is an antiderivative of f if F'(x) = f(x) for all x in the
domain of f(x). (Also stated here:
http://planetmath.org/encyclopedia/Antiderivative.html .)
Either way, the property F'(x) = f(x) is not necessarily true for piecewise
antiderivatives defined that way. Consider this function.
f(x) = x, 0 = x = 1
f(x) = 1, 1 x
If we use the definition you gave to find F = integral(f), F'(1) is
undefined so it is not true that F'(x) == f(x) for all x.
Instead, we use the definition that F=
integrate(f1, t, a1, x), a1 x = a2
integrate(f2, t, a2, x) + integrate(f1, t, a1, a2), a2 x = a3
integrate(f3, t, a3, x) + integrate(f2, t, a2, a3) + integrate(f1, t, a1,
a2), a3 x = a4
...
integrate(fn, t, an, x) + integrate(f[n-1], t, a[n-1], an) + ... +
integrate(f1, t, a1, a2), an x
(We also need a special case for when a1 = -infinity, which I didn't show.)
With this definition, F(b) - F(a) can be used to find the Riemann sum
between a and b. Also, F'(x) = f(x) seems to hold, except at points where
f(x) goes from defined to undefined or vice-versa.
The antiderivative is only well-defined up to an additive constant.
IMHO, the piecewise defined function of antiderivavtives
int f1(x) dx +C1 , a1x=a2,
int f2(x) dx +C2, a2x=a3,
...
int fn(x) dx +Cn, anx=a{n+1}
does not make sense.
I agree that it doesn't make sense where C1 .. Cn are arbitrary constants.
-- Paul
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-devel] Re: Integral of piecewise functions
When a1 = -infinity, I would make F1 = integrate(f, x, a2, x) instead of
integrate(f, x, a1, x). Then I would not calculate the definite integral of
the first interval, which would align my constants so that F(a2) = 0. When I
get a chance, I'll add this to my code.
Functions like floor with an infinite number of pieces are beyond the scope
of the Piecewise class, because piecewise functions in Sage can only (at the
moment) have finitely many pieces. I can't give an algorithm for integration
of all piecewise functions with multiple pieces off the top of my head, but
I'll give it some thought.
-- Paul
On Sat, Dec 6, 2008 at 4:13 PM, David Joyner [EMAIL PROTECTED] wrote:
Okay, this helps me understand what you mean.
Still, the case a1=-ifinty is precisely the special case which I don't
understand.
For example, take a function such as f(x) = max(1,floor(x)), x real.
How do you define an antiderivative F(x) so that
F(b)-F(a) = area under the y=f(x) for axb?
(And mayeb you can do it for that special function,
and let us ignore points of discontinuity for the sake of
discussion.) In other words, I am asking for the algorithmic procedure
you would use to create an area function of a piecewise-defined
function on the reals.
With this definition, F(b) - F(a) can be used to find the Riemann sum
between a and b. Also, F'(x) = f(x) seems to hold, except at points where
f(x) goes from defined to undefined or vice-versa.
The antiderivative is only well-defined up to an additive constant.
IMHO, the piecewise defined function of antiderivavtives
int f1(x) dx +C1 , a1x=a2,
int f2(x) dx +C2, a2x=a3,
...
int fn(x) dx +Cn, anx=a{n+1}
does not make sense.
I agree that it doesn't make sense where C1 .. Cn are arbitrary
constants.
-- Paul
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-devel] Re: Expected value of probability space
Hi David,
Thanks for explaining that, I see how that causes problems when S is not a
set of numbers. Even so, would it make sense for the random variable ps to
be the identity function X(x) = x on the probability space ps? Currently the
random variable ps is the function X(x) = P(x). Is this a useful random
variable that I'm just not aware of?
One way to accomplish this is to override expectation in the
DiscreteProbabilitySpace class with the one you sent earlier. We would also
need to do this for variance, covarience, etc., so I admit this is probably
not the ideal solution.
-- Paul
On Mon, Dec 1, 2008 at 3:48 AM, William Stein [EMAIL PROTECTED] wrote:
[Another response in this thread from David Kohel (who maybe should be
posting on list)]
Hi William and Paul,
Actually, I correct myself -- the average should be over the values
of the function, weighted by the probabilities. The domain of the
function (the keys) can be in any set (e.g. A,B,C), so the
current behavior is correct.
In the defintition of the probability space itself:
sage: ps = DiscreteProbabilitySpace([1,2,3],{1:1/3,2:1/3,3:1/3})
Consider replacing the domain of ps with such a set. Then it
should be clear that you can't average over A, B, and C.
S = [A,B,C]
P = {}
for i in range(3):
P[S[i]] = 1/3
ps = DiscreteProbabilitySpace(S,P)
ps.expectation() # 0.33
This is the random variable with
f(A) = 1,
f(B) = 2,
f(C) = 3,
for which the expectation is 2:
f = {}
for i in range(3):
f[S[i]] = i+1
rv = DiscreteRandomVariable(ps,f)
rv.expectation() # 2.00
On the other hand, I'd be happy with some syntax for creating a
random variable with the uniform distribution P(x) = 1/n on some
set or tuple. Currently this gives an error:
sage: rv = DiscreteRandomVariable(S,f)
but it could easily create the uniform distribution on S behind
the scenes. The reason I didn't do this is the following:
sage: ps1 = DiscreteProbabilitySpace(S,P)
sage: ps2 = DiscreteProbabilitySpace(S,P)
sage: ps1 == ps2
False
Probably this should be considered a bug. But to avoid creating
many copies of the same space, one should consider caching the
probability spaces. However, since dictionaries are mutable, it
can't be cached.
The correct solution might be to define additional types for the
uniform distribution and other standard distributions (which are
the most likely candidates to be created over and over) and have
only one instance of each.
My type DiscreteProbabilitySpace is also too naive -- it really
only implements a FiniteProbabilitySpace.
It is suitable for the examples I've taught in classical cryptography,
but for a statistics class one certainly needs infinite discrete
probability spaces and more classes for distributions on continuous
intervals.
Some thought also needs to be given to their ease of use. The
usual approach to this is to first write some documentation which
describes the desired behavior and then implement the classes and
functions.
--David
- Forwarded message from David R. Kohel [EMAIL PROTECTED]
-
Date: Mon, 1 Dec 2008 09:04:12 +0100
From: David R. Kohel [EMAIL PROTECTED]
To: William Stein [EMAIL PROTECTED]
Cc: Paul Butler [EMAIL PROTECTED]
Subject: Re: Fwd: [sage-devel] Expected value of probability space
User-Agent: Mutt/1.5.6i
Dear William, Paul,
Indeed, the function definition should be:
def expectation(self):
r
The expectation of the discrete random variable, namely
$\sum_{x \in S} p(x) X[x]$,
where $X$ = self and $S$ is the probability space of $X$.
E = 0
Omega = self.probability_space()
for x in self._function.keys():
E += Omega(x) * x
return E
rather than:
def expectation(self):
r
The expectation of the discrete random variable, namely
$\sum_{x \in S} p(x) X[x]$,
where $X$ = self and $S$ is the probability space of $X$.
E = 0
Omega = self.probability_space()
for x in self._function.keys():
E += Omega(x) * self(x)
return E
Cheers,
David
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
[sage-devel] Re: Expected value of probability space
Hi David, When I was referring to the probability space ps and the random variable ps, I was referring to the fact that ps is by inheritance a probability space and a random variable (is_DiscreteProbabilitySpace(ps) == is_DiscreteRandomVariable(ps) == True). I do see that the values are necessarily probabilities, I'm not suggesting that that be changed. But right now if you ask for the expected value of ps, it uses the probability function (from the constructor) as both the probability function and the random variable. In other words, currently ps.expectation() = sum(P(x)^2 over all x). I only brought it up because after looking at the code I was unsure if it was intentional or just a side-effect of the class organization. I'm not aware of any formal definition of the expectation of a probability space, so I wanted to ask to clarify. You mentioned that continuous random variables are missing, which is something I wouldn't mind working on. I may follow up later with an email to get some input on what needs to be done, if you don't mind. -- Paul On Mon, Dec 1, 2008 at 4:50 PM, David Kohel [EMAIL PROTECTED] wrote: Hi Paul, Thanks for explaining that, I see how that causes problems when S is not a set of numbers. Even so, would it make sense for the random variable ps to be the identity function X(x) = x on the probability space ps? Currently the random variable ps is the function X(x) = P(x). Is this a useful random variable that I'm just not aware of? Well, ps is a probability space (you did create it as DiscreteProbabilitySpace, rather than DiscreteRandomVariable, which I had missed in my first reply), hence its values are necessarily probabilities. There is no other valid choice. To create a DiscreteRandomVariable, you currently must first create a probability space, then the function on that space. One could have shorter constructors which assume a finite uniform probability space, if not given. The order of the arguments is probability space, then function or values, but if a probability space is no longer required, this order should change. --David --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~--~~~~--~~--~--~---
[sage-devel] Expected value of probability space
I've been experimenting with probability and found that in Sage, a
probability space is also a random variable by inheritance. This may be
useful. Without it, creating a random variable requires two classes: a
probability space and a random variable on that probability space.
Unfortunately, the random variable doesn't work like I expected. For
example:
sage: ps = DiscreteProbabilitySpace([1,2,3],{1:1/3,2:1/3,3:1/3})
sage: ps.expectation()
0.333
(I expected 2.00)
I've prepared a patch that gives me the value I'd expect, but I'd like to
make sure this is the proper behavior.
-- Paul
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~--~~~~--~~--~--~---
