Convolution of piecewise polynomial functions

[Co-Processor]

Given two piecewise polynomial functions, how could one compute the (symbolic representation of the) function generated by the convolution of these functions in REDUCE? E.g. I want to get the (piecewise polynomial) probability density function of X1+X2 given the probability density functions of X1,X2, which are piecewise polynomial.

[Co-Processor]

I guess it is not possible to use any of the existing possibilities to model piecewise functions (let {... => e1 when ...} or procedure ...; if ... then e1 ...) because the limits of the scope of validity of the single expressions cannot be retrieved, but they are necessary to compute the integration limits. So the way to go seems to be to model a piecewise function f as a list of pairs like {{l1,e1},{l2,e2}, ...} with the intention of let {f(~x) => e1 when x <= l1, f(~x) => e2 when l1 < x <= l2, ...}. And then one can use REDUCE to compute the integrals int(sub({x=t},ei)*sub({x=t-x},ej),t,...) over appropriate intervals. Or is there any smarter possibility?

[Arthur Norman]

On Thu, 24 Sep 2020, Co-Processor wrote:

I guess it is not possible to use any of the existing possibilities to
model piecewise functions (let {... => e1 when ...} or procedure ...; if ... then e1 ...) because the limits of the scope of validity of the
single expressions cannot be retrieved, but they are necessary to compute
the integration limits. So the way to go seems to be to model a piecewise
function f as a list of pairs like {{l1,e1},{l2,e2}, ...} with the
intention of let {f(~x) => e1 when x <= l1, f(~x) => e2 when l1 < x <= l2, ...}. And then one can use REDUCE to compute the integrals
int(sub({x=t,ei})*sub({x=t-x},ej),t,...) over appropriate intervals. Or
is there any smarter possibility?

I think your suggestion using a list is liable to be as easy as
anything... most certainly in the short term.

One can insepect the slightly-more-internal form that Reduce parses things
into using "on defn;" and so (e.g).

3: on defn;
4: let (a(~x) => e1 when x < 0, a(~x) => e2 );
(let '((replaceby (a (!~ x)) (when e1 (lessp x 0))) (replaceby (a (!~ x))
e2)))

so in certain contexts (a(~x) => e1 when x < 0, a(~x) => e2) cane be
parsed, eg I can go

16: q (a(~x) => e1 when x < 0, a(~x) => e2 );
(aeval
(list
'q
(list
'replaceby
(list 'a (list '!~ 'x))
(list 'when 'e1 (list 'lessp 'x 0)))
(list 'replaceby (list 'a (list '!~ 'x)) 'e2)))

but it looks as if I can not write the list of cases stand-alone.
If you (or somebody else) became keen to implement a package for "doing
things" with piecewise-defined functions this would give a starting point
for the input fragments, and the stuff in doc/primers might help further.

Deciding on exact syntax such that the new syntax was not ambiguous
against any existing would also be a good start, and tnen enhancing the
parser. Which for unambiguous cases is typically not too hard in
packages/rlisp/*.red.

I might like to be able to write
... f(~x) => x+1 when x<0, x+1 when x<10, x^2 ...
and if the word "when" or use of commas leads to parsing pain to pick on
some other notation.

Code such as the taylor series packages illustrate the introduction of new
"domains" and "piecewise" could be similar in various respects to
"truncated power series" in overall status within the system.

You prototyping stuff using just lists could be a good start to show the
value of the concept and if what you create would be useful to others we
could incororate it as an extra package...

Arthur

---

Convolution of piecewise polynomial functions

---

Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/reduce-algebra/discussion/899364/

To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

[Co-Processor]

Ok, I managed to implement a first version of a convolution procedure pw_conv for piecewise functions represented as list as proposed above. The standard uniform distribution is represented as {{0,0},{1,1},{INFINITY,0}} (a trailing {INFINITY,0} could be omitted). The successive convolutions with itself give the expected results {{0,0},{1,x},{2,2-x}}, ... (see attached code + screenshot).
Now I am stuck on how to convert this representation to a procedure for numerical evaluation and plotting. Is it possible to generate a procedure from this representation of piecewise functions?

[Co-Processor]

Here is the screenshot:

[Arthur Norman]

The code here does not do what you want (!) because I did not think at all
hard about where I thought I wanted to be in terms of the ranges. And
Reduce is unhappy about attempts to compare "4 < infinity" so I made the
thing detect your (perhaps unnecessary!) final list item explicitly.

But this may give you a START.

Arthur

a := {{0, 0}, {1, x}, {2, -x-2}, {infinity, 0}};
a := {{0,0},{1,x},{2, - x - 2},{infinity,0}}

procedure piecewiseeval(e, x);
if e = {} or first first e = infinity then 0
else if x < first first e then second first e
else piecewiseeval(rest e, x);
piecewiseeval

for i := -2:4 do
write i, " ", piecewiseeval(a, i);
-2 0

-1 0

0 x

1 - x - 2

2 0

3 0

4 0

end;

On Fri, 25 Sep 2020, Co-Processor wrote:

Ok, I managed to implement a first version of a convolution procedure pw_conv for piecewise functions represented as proposed above. The standard uniform distribution is represented as {{0,0},{1,1},{INFINITY,0}} (a trailing {INFINITY,0} could be omitted). The successive convolutions with itself give the expected results {{0,0},{1,x},{2,2-x}}, ... (see attached code + screenshot).
Now I am stuck on how to convert this representation to a procedure for numerical evaluation and plotting. Is it possible to generate a procedure from this representation of piecewise functions?

[Eberhard Schruefer]

Using Arthur's procedure, wouldn't just

procedure piecewisevalue(e,v);
   sub(x=v, piecewiseeval(e, v));

procedure piecewiseeval(e, x);
   if e = {} or first first e = infinity then 0
    else if x < first first e then second first e
    else piecewiseeval(rest e, x);

do what you want? E.g.

48: pwc := {{0,0},{1,x**2/2},{2,( - 2*x**2 + 6*x - 3)/2},{3,(x**2 - 6*x + 9)/2},{infinity,0}};

pwc := {{0,0},

                2
        {1,0.5*x },

               2
        {2, - x  + 3*x - 1.5},

                2
        {3,0.5*x  - 3*x + 4.5},

        {infinity,0}}

49: 
49: on rounded;

50: piecewisevalue(pwc,1.99);

0.5099

51: piecewisevalue(pwc,2.01);

0.49005

Eberhard

[Co-Processor]

Basically that would do it. For plotting one would have to generate points and pass them to the plot command, if I got that correctly. But I was looking for a way of (programmatically) generating a function/procedure mapping the x-value to the y-value. In python one could do something like this (using strings here for simplicity instead of symbolic expressions like provided by sympy):

def genPieceFunc(pdf:list):
  def pieceFunc(x):
    for upper,expr in pdf:
      if x < upper: return eval(expr, {'x':x})
    return 0
  return pieceFunc

pdf = [[0,0], [1,'x'], [2,'2-x']]
pdfFunc = genPieceFunc(pdf)

for x in [0,0.5,1,1.5,2,2.5]:
  print(x,pdfFunc(x))

In REDUCE one would have to go to the Lisp level to do something like this?