JoyDivisionMan

50 Reputation

3 Badges

5 years, 4 days

MaplePrimes Activity


These are replies submitted by JoyDivisionMan

@sand15 

Thank you for your help.  I now understand how 2*Convolution of X and Y can result in the density function of Z, where Z = |X - Y|.  In fact, a few months ago, I wrestled with this same issue when considering Z = |X - Y| when X and Y were both uniformly distributed on the [0,1] interval.

Of course, X and Y both being uniformly distributed on [0,1] is more complicated than my current dilemma.

I looked at the worksheet you provided.  Some of the things you asked for I did not understand.  I pruned away some of the items that were not necessary for me.

I understand this is to simulate values of Z:

S__Z := Sample(Z, 10^6):

plots:-display(
  Histogram(S__Z, transparency=0.7, minbins=100, style=polygon)
  , plot(f__U(u), u=0..2, color=red, thickness=2, legend=typeset(:-PDF('Z')))
)

I did not understand the section below, but I think it is a numerical estimate.

F__U := proc(u)
  option remember:
  local F1;
  F1 := 2*evalf(Int(op([2, 2], f__U(t)), t=0..1)):
  if _passed::numeric then 
    if u <= 1 then
      2*evalf(Int(op([2, 2], f__U(t)), t=0..u))
    elif u <= 2 then
      F1 + 2*evalf(Int(Re(op([2, 4], f__U(t))), t=1..u, method = _d01ajc))
    end if:
  else
    procname(_passed)
  end if:
end proc:

plot([seq([u, F__U(u)], u=0..2, 0.1)], color=red, thickness=2, legend=typeset(:-CDF('Z')), gridlines=true)
Error, (in F__U) improper op or subscript selector

I (think I) understand the below to be the exact Mean of Z, which is great.

:-Mean('Z') = Mean(Z);

But when I ask for the PDF of Z, I get a bunch of garbage.  But I do get the expected result when I entered plot(f__U(u), u = 0 .. 2).

My modification of your work is included in PDFDerivation.mw.

As such, I moved ahead with convolving X and Y to see what would happen.  Unfortunately, I was unable to integrate the product of the X and Y.  This work is also attached (Convolution.mw).  I "shifted" the t variable to the right so that t is always positive.

restart

kernelopts(version)

`Maple 2023.2, X86 64 WINDOWS, Nov 24 2023, Build ID 1762575`

(1)

with(Statistics):

f := x -> piecewise(abs(x) < 1, 2*sqrt(1 - x^2)/Pi, 0);

f := proc (x) options operator, arrow; piecewise(abs(x) < 1, 2*sqrt(1-x^2)/Pi, 0) end proc

(2)

F := unapply(int(f(t), t=-1..x), x);

F := proc (x) options operator, arrow; piecewise(x <= -1, 0, x <= 1, (1/2)*(2*x*sqrt(1-x^2)+Pi+2*arcsin(x))/Pi, 1 < x, 1) end proc

(3)

d := Distribution(
       PDF = (x -> f(x)),
       CDF = (x -> F(x))
     ):

Y__1 := RandomVariable(d):
Y__2 := RandomVariable(d):

# From symmetry considerations the PDF of Z is twice the PDF of U restricted to U > 0

Z := abs(U):
U := Y__1 - Y__2:
 

 

# To help Maple find the correct PDF we may observe that, from symmetry considerations,
# the PDF of Z is twice the PDF of U = Y__1 - Y__2 once restricted to U > 0

U := Y__1 - Y__2:

f__U := unapply( (2*PDF(U, u) assuming u > 0), u);

f__U := proc (u) options operator, arrow; (4/3)*piecewise(u <= 1, ((u^2+4)*EllipticE((-2+u)/(u+2))-4*EllipticK((-2+u)/(u+2))*u)*(u+2), u <= 2, u*(u^2*EllipticE(I*sqrt(-u^2+4)/u)+4*EllipticE(I*sqrt(-u^2+4)/u)-8*EllipticK(I*sqrt(-u^2+4)/u)), 2 < u, 0)/Pi^2 end proc

(4)

S__Z := Sample(Z, 10^6):

plots:-display(
  Histogram(S__Z, transparency=0.7, minbins=100, style=polygon)
  , plot(f__U(u), u=0..2, color=red, thickness=2, legend=typeset(:-PDF('Z')))
)

 

F__U := proc(u)
  option remember:
  local F1;
  F1 := 2*evalf(Int(op([2, 2], f__U(t)), t=0..1)):
  if _passed::numeric then
    if u <= 1 then
      2*evalf(Int(op([2, 2], f__U(t)), t=0..u))
    elif u <= 2 then
      F1 + 2*evalf(Int(Re(op([2, 4], f__U(t))), t=1..u, method = _d01ajc))
    end if:
  else
    procname(_passed)
  end if:
end proc:

plot([seq([u, F__U(u)], u=0..2, 0.1)], color=red, thickness=2, legend=typeset(:-CDF('Z')), gridlines=true)

Error, (in F__U) improper op or subscript selector

 

# Because Y__1 and Y__2 both have the same symmetric PDF:
#
#      U = Y__1 - Y__2 = Y__1 + (-Y__2) = Y__1 + Y__2
#
# Thus the PDF of U can be obtained through a convolution product

f__Y := unapply(PDF(Y__1, y), y):

C := unapply( int(f__Y(y)*f__Y(u-y), y = -infinity..+infinity), u ):

plots:-display(
  Histogram(S__Z, transparency=0.7, minbins=100, style=polygon)
  , plot(2*C(u), u=0..2, color=red, thickness=2, legend=typeset(:-PDF('Z')))
)

 

Depending on what statistic you want it is sometimes more powerfull to use U instead of Z

:-Mean('Z') = Mean(Z);

Mean(Z) = (256/45)/Pi^2

(5)

:-Variance('Z') = Variance(Z);

Variance(Z) = (1/4050)*(2025*Pi^4-131072)/Pi^4

(6)

2*PDF(U,u);

4*piecewise(u <= -2, 0, u <= -1, -((u^2+4)*EllipticE((2+u)/(-2+u))+4*EllipticK((2+u)/(-2+u))*u)*(-2+u), u <= 0, (1/2)*(-2*u^3-8*u)*EllipticE(I*u/(sqrt(-u+2)*sqrt(2+u)), I*sqrt(-u^2+4)/u)+8*EllipticF(I*u/(sqrt(-u+2)*sqrt(2+u)), I*sqrt(-u^2+4)/u)*u+infinity, u <= 1, ((u^2+4)*EllipticE((-2+u)/(2+u))-4*EllipticK((-2+u)/(2+u))*u)*(2+u), u <= 2, u*(u^2*EllipticE(I*sqrt(-u^2+4)/u)+4*EllipticE(I*sqrt(-u^2+4)/u)-8*EllipticK(I*sqrt(-u^2+4)/u)), 2 < u, 0)/(3*Pi^2)

(7)

plot(f__U(u), u = 0 .. 2)

 

NULL


 

Download PDFDerivation.mw
 

restart

with(Statistics)

X := 2*sqrt(1-(1-tau)^2)/Pi

2*(1-(1-tau)^2)^(1/2)/Pi

(1)

Y := 2*sqrt(1-(1-t+tau)^2)/Pi

2*(1-(1-t+tau)^2)^(1/2)/Pi

(2)

CF := int(X*Y, tau = 0 .. t)

int(4*(1-(1-tau)^2)^(1/2)*(1-(1-t+tau)^2)^(1/2)/Pi^2, tau = 0 .. t)

(3)
 

 

Download Convolution.mw

@sand15 

Thank you for your reply.  

Y1 and Y2 DO have the same distribution.  I was incorrect when I stated my simulated result of |Y1 - Y2| takes on a logrithmic form.  I had the same result as your simulated PDF.

I need to find the exact PDF of |Y1 - Y2|.

I know how to sum PDFs to find the new PDF.  I know how to multiply PDFs to find the new PDF.  But the absolute difference between two PDFs into a new PDF is something new to me.

@mmcdara 

Thank you -- this is helpful.

@acer 

Thank you.  I do not need the symbolic solution -- I was hoping it would be "tidier." 

The evalf solution was just a check to see if the result was as expected when compared to the simulated result.

I am good.  Thanks again.

@mmcdara 

Thank you for your help with this issue.  I have had success approximating PDFs with the Pearson Distribution.  This is what you suggested:  PDF_analytic_approximation.mw

I have recently attempted to approximate another random variable that I wasn't able to simulate via Maple.  I had to simulate the data via Java.  I then imported the simulated data into Maple with the intent of approximating the PDF via the Pearson Distribution.  I was successful in bringing the data into Maple via the Excel import and create a histogram.  I ran into problems when I attempted to approximate the Pearson disbtibution in Maple.  I think the "data type" is the problem.  The data was imported as a matrix.  I then tried to convert the data to Vector and Array form, but without success.

The data and my attempt at approximating the Pearson Distribution as as follows:  MyData.xlsxMyApproximationAttempt.mw.

@mmcdara thank you.  I had no idea these integration tools are available.

This brings me to a follow-up question.  This density function:  fw:= 2*sqrt(1-t2)/Pi...can it be classified as a Random Variable?

Thank you, acer.  I did not realize that the PDF function could be followed by the "assuming" keyword.  

My guess would be that my omitting this ended up confusing Maple.

All of this is very helpful.  Thank you for your responses.

Thank you.  This was the problem.  The current directory was set to a folder that no longer exists.

Thank you very much -- this is exactly what I needed.

Page 1 of 1