Advanced Search

Question: Why does "discont" give this result

Question:Why does "discont" give this result

mmcdara 7359 Maple 2015
piecewise discont inttrans heaviside invfourier
October 27 2022
0

I'm working on finding the analytic expression of the PDF of a sum of abstract Uniform Random Variables URV).
Here "abstract" means that the supports are not numeric but litteral.

Maple is capable to find such a PDF for numeric supports but unable to determine the PDF of U1+U2 where 

U1 := RandomVariable(Uniform(a1, b1)):
U2 := RandomVariable(Uniform(a2, b2)):

A way to deal with abstract URV is to complute explicitely the convolution product of th PDFs.
Ir seems that this fails (MAPLE 2015.2) for these PDF are piecewise functions.
A workaround is to convert them first into Heaviside(s).

One done the explicit expression of the convolution product can be obtained for a sum of 2 abstract URVs, but not for a sum of a larger number of abstract URVs.

Thus the second workaround which consists in using direct and inverse Fourier transform.

The question is :
obviously, the PDF f(t ; a1...aN, b1...bN) of  U1 + ... + UN is a continuous function of t: why does discont(f(t ; ...), t) returns the non empty set of the values where the Heaviside functions are undefined ?

Here is a very simple result

u1 := RandomVariable(Uniform(-1, 1)): 
u2 := RandomVariable(Uniform(-1, 1)):
p := PDF(u1+u2, t):
discont(p, t);

print("-------------------------------------");

f1 := convert(PDF(u1, t), Heaviside):
f2 := convert(PDF(u2, t), Heaviside):
g1 := fourier(f1, t, xi):
g2 := fourier(f2, t, xi):
g  := g1*g2:
f  := invfourier(g, xi, t):
discont(f, t);

                               {}
            "-------------------------------------"
                           {-2, 0, 2}

The problem is (IMO) that discont(f, t) should return { }, but that some function (does it exists) should say that d is undefined at points t=-2, t=0, t=2.
The output of discont(f, t) doesn't seem consistent with the definition of the continuity

limit(f, t=-2, left);
limit(f, t=-2, right);
eval(f, t=-2)
                               0
                               0
                           undefined
> 

restart:
> 

with(inttrans);
with(Statistics):

[addtable, fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, invhilbert, invlaplace, invmellin, laplace, mellin, savetable]

 (1)
> 

N := 3:
for n from 1 to N do
  U||n := RandomVariable(Uniform(a__||n, b__||n)):
end do;

_R

  

_R0

  

_R1

 (2)
> 

# Maple fails to compute the PDF of a sum of abstract (meaning with symbolic support) uniform RVS
# PDF(U1+U2, t)
> 

# For N <=3 the computation of the convolution product is possible.
#
# For N > 4 (Maple 2015) it seems this is no longer the case. The trick used here is based on
# the fact that Fourier(Conv(f, g)) : Fourier(f)*Fourier(g)
# Thus PDF(U1+U2) = conv(PDF(U1), PDF(U2)) = invFourier(Fourier(PDF(U1)).Fourier(PDF(U1)))


for n from 1 to N do
  f||n := convert(PDF(U||n, t), Heaviside):
end do:

for n from 1 to N do
  g||n := fourier(f||n, t, xi):
end do:
g := mul(g||n, n=1..N):

hyp := seq(b__||n > a__||n, n=1..N);
f   := invfourier(g, xi, t) assuming hyp;

a__1 < b__1, a__2 < b__2, a__3 < b__3

  

((1/2)*(t-a__1-a__2-a__3)^2*Heaviside(-t+a__1+a__2+a__3)-(1/2)*(t-a__1-a__2-b__3)^2*Heaviside(-t+a__1+a__2+b__3)-(1/2)*(t-a__1-a__3-b__2)^2*Heaviside(-t+a__1+a__3+b__2)+(1/2)*(t-a__1-b__2-b__3)^2*Heaviside(-t+a__1+b__2+b__3)-(1/2)*(t-a__2-a__3-b__1)^2*Heaviside(-t+a__2+a__3+b__1)+(1/2)*(t-a__2-b__1-b__3)^2*Heaviside(-t+a__2+b__1+b__3)+(1/2)*(t-b__1-a__3-b__2)^2*Heaviside(-t+b__1+a__3+b__2)-(1/2)*(t-b__1-b__2-b__3)^2*Heaviside(-t+b__1+b__2+b__3))/((-b__1+a__1)*(-b__2+a__2)*(-b__3+a__3))

 (3)
> 

# f is obviously a continuous function of t, but I get this strange result

discont(f, t);

{a__1+a__2+a__3, a__1+a__2+b__3, a__1+a__3+b__2, a__1+b__2+b__3, a__2+a__3+b__1, a__2+b__1+b__3, b__1+a__3+b__2, b__1+b__2+b__3}

 (4)
> 

# note that this "error" also appears if the a__n's and b__n's are numeric

r := rand(0. .. 1.):
P := convert(indets(g, name) minus{xi}, list):
E := NULL:
for n from 1 to N do
  a := r():
  b := a + r():
  E := E, a__||n = a, b__||n = b
end do:
E := [E];

G := eval(g, E):
f := invfourier(G, xi, t);
discont(f, t);

[a__1 = .3055679837, b__1 = .7906643786, a__2 = .8311025583, b__2 = .9857257095, a__3 = .4223879539, b__3 = .7034757826]

  

0.1333206533e-8*(-0.1000000000e11*t+0.1136670542e11)*Heaviside(-t+1.136670542)+0.1333206533e-8*(0.1000000000e11*t-0.1291293693e11)*Heaviside(-t+1.291293693)+0.1333206533e-8*(0.1000000000e11*t-0.1621766937e11)*Heaviside(-t+1.621766937)+0.1333206533e-8*(-0.1000000000e11*t+0.1776390088e11)*Heaviside(-t+1.776390088)

  

{1.136670542, 1.291293693, 1.621766937, 1.776390088}

 (5)
> 

u1 := RandomVariable(Uniform(-1, 1)):
u2 := RandomVariable(Uniform(-1, 1)):
p := PDF(u1+u2, t):
discont(p, t);

print("-------------------------------------");

f1 := convert(PDF(u1, t), Heaviside):
f2 := convert(PDF(u2, t), Heaviside):
g1 := fourier(f1, t, xi):
g2 := fourier(f2, t, xi):
g  := g1*g2:
f  := invfourier(g, xi, t):
discont(f, t);
 

{}

  

"-------------------------------------"

  

{-2, 0, 2}

 (6)
> 

limit(f, t=-2, left);
limit(f, t=-2, right);
eval(f, t=-2);

print("-------------------------------------");

limit(f, t=0, left);
limit(f, t=0, right);
eval(f, t=0)

0

  

0

  

undefined

  

"-------------------------------------"

  

1/2

  

1/2

  

undefined

 (7)
> 

 

Download Sum_of_Uniform_RVs.mw
Please Wait...