JoyDivisionMan

45 Reputation

3 Badges

4 years, 288 days

MaplePrimes Activity


These are questions asked by JoyDivisionMan

I constructed a density function, and I am certain it shows me what I want.  The problem I am having is parsing the Elliptic functions. Is there a way to "get rid" of the ones I don't want or need.

I generated a plot of the function -- the plot tells me what I expected based on simulation. I need to know if there is a way to express the density function (y) as a function of t and without the elliptic functions..  Even a numerical solution would be fine.

I assume the denomenator term is correct. I also assume that I don't need complex values. My input file is below.

Basics.mw

I have a piecewise density function (f).  I am trying to find the median value.

I have tried the Median function and the Percentile function, but neither work for me.  I am not sure why.

I have also tried to integrate the density function on the (0,x) interval such that the area under the curve is 1/2 and then solve for x.  This works for simple problems, but not the one attached.

There is something very simple that I'm not doing, but I am not sure what.

restart

kernelopts(version)

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

(1)

with(Statistics)

 

(2)

f := 2*t*piecewise(t <= 0, 0, t <= 1, t^2+Pi-4*t, t <= sqrt(2), -(sqrt(t^2-1)*t^2+2*sqrt(t^2-1)*arcsin((t^2-2)/t^2)-4*t^2+2*sqrt(t^2-1)+4)/sqrt(t^2-1), sqrt(2) < t, 0)

2*t*piecewise(t <= 0, 0, t <= 1, t^2+Pi-4*t, t <= 2^(1/2), -((t^2-1)^(1/2)*t^2+2*(t^2-1)^(1/2)*arcsin((t^2-2)/t^2)-4*t^2+2*(t^2-1)^(1/2)+4)/(t^2-1)^(1/2), 2^(1/2) < t, 0)

(3)

plot(f, t = 0 .. sqrt(2))

 

Median(f)

2*t*piecewise(t <= 0, 0, t <= 1, t^2+Pi-4*t, t <= 2^(1/2), -((t^2-1)^(1/2)*t^2+2*(t^2-1)^(1/2)*arcsin((t^2-2)/t^2)-4*t^2+2*(t^2-1)^(1/2)+4)/(t^2-1)^(1/2), 2^(1/2) < t, 0)

(4)

Percentile(f, 50)

2*t*piecewise(t <= 0, 0, t <= 1, t^2+Pi-4*t, t <= 2^(1/2), -((t^2-1)^(1/2)*t^2+2*(t^2-1)^(1/2)*arcsin((t^2-2)/t^2)-4*t^2+2*(t^2-1)^(1/2)+4)/(t^2-1)^(1/2), 2^(1/2) < t, 0)

(5)

myMedian := solve(int(f, t = 0 .. x) = 1/2, x)

Warning, solutions may have been lost

 

(6)

NULL

Download Median.mw

I am trying to find the standard deviation for a piecewise density function (f).  The interval is (0, sqrt(3))  I am convinced the density function is what I want.  I calculate the expected value of the density function (evE) and the answer is correct -- about 0.66145

When I try to calculate the standard deviation (stdE), I get an answer that is "off" by a large degree.  Via simulating values, I should get a value of about 0.24936.  Each time I "re-run" the calculation, I get varying results, all of which are "off" by a large degree.

I am only guessing, but my integration function might be missing some sort of assumption and/or option.

My work is attached.  Does anyone know what I am doing wrong?

restart

with(Statistics)

f := -t*piecewise(t <= 0, 0, t <= 1, t*(t^3+6*Pi*t-8*t^2-4*Pi), t <= sqrt(2), -(2*(4*t^2*sqrt(t^2-1)*arcsin((t^2-2)/t^2)+4*t^2*arctan(sqrt(t^2-1))*sqrt(t^2-1)+(t^4+(2*Pi+3)*t^2-4*t*Pi+3*Pi-1/2)*sqrt(t^2-1)-8*t^4+4*t^2+4))/sqrt(t^2-1), t <= sqrt(3), (8*arctan((t^3+t^2-3*t-1)/(sqrt(t^2-2)*(t^3-t^2-t-1)))*t*sqrt(t^2-2)-4*arccot((t^2-t-1)/sqrt(t^2-2))*t^2*sqrt(t^2-2)+4*t*sqrt(t^2-2)*(t+2)*arctan((t^2+t-1)/sqrt(t^2-2))+(8*t^2-4)*sqrt(t^2-2)*arcsin((t^2-3)/(t^2-1))+16*arcsin(sqrt(t^2-2)/sqrt(t^2-1))*sqrt(t^2-2)-16*arcsin(1/sqrt(t^2-1))*sqrt(t^2-2)+((t^2+5)*sqrt(t^2-2)-8*t^2+16)*(t^2+1))/sqrt(t^2-2), sqrt(3) < t, 0)

-t*piecewise(t <= 0, 0, t <= 1, t*(t^3+6*Pi*t-8*t^2-4*Pi), t <= 2^(1/2), -2*(4*t^2*(t^2-1)^(1/2)*arcsin((t^2-2)/t^2)+4*t^2*arctan((t^2-1)^(1/2))*(t^2-1)^(1/2)+(t^4+(2*Pi+3)*t^2-4*t*Pi+3*Pi-1/2)*(t^2-1)^(1/2)-8*t^4+4*t^2+4)/(t^2-1)^(1/2), t <= 3^(1/2), (8*arctan((t^3+t^2-3*t-1)/((t^2-2)^(1/2)*(t^3-t^2-t-1)))*t*(t^2-2)^(1/2)-4*arccot((t^2-t-1)/(t^2-2)^(1/2))*t^2*(t^2-2)^(1/2)+4*t*(t^2-2)^(1/2)*(t+2)*arctan((t^2+t-1)/(t^2-2)^(1/2))+(8*t^2-4)*(t^2-2)^(1/2)*arcsin((t^2-3)/(t^2-1))+16*arcsin((t^2-2)^(1/2)/(t^2-1)^(1/2))*(t^2-2)^(1/2)-16*arcsin(1/(t^2-1)^(1/2))*(t^2-2)^(1/2)+((t^2+5)*(t^2-2)^(1/2)-8*t^2+16)*(t^2+1))/(t^2-2)^(1/2), 3^(1/2) < t, 0)

(1)

plot(f, t = 0 .. sqrt(3))

 

evE := simplify(int(t*f, t = 0 .. sqrt(3)))

(1/105)*(294*(2^(1/2)-3/2)*(-2+3^(1/2))*ln(1+2^(1/2))+168*(2^(1/2)-3/2)*(-2+3^(1/2))*ln(1+3^(1/2))-252*(2^(1/2)-3/2)*(-2+3^(1/2))*arccoth(2^(1/2))+((-14*Pi-84*ln(2)-19)*3^(1/2)+28*Pi+168*ln(2)+50)*2^(1/2)+(21*Pi+126*ln(2)+20)*3^(1/2)-42*Pi-252*ln(2)-58)/((2*2^(1/2)-3)*(-2+3^(1/2)))

(2)

evalf(evE)

.6617071838

(3)

stdE := simplify(sqrt(int(f*(t-evE)^2, t = 0 .. sqrt(3))))

(1/210)*(-4808607312140100*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*ln(1+2^(1/2))^2+(-5495551213874400*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*ln(1+3^(1/2))+8243326820811600*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*arccoth(2^(1/2))+((457962601156200*Pi+2747775606937200*ln(2)+631290043332900)*3^(1/2)-793214493249360*Pi-4759286959496160*ln(2)-1093426429714584)*2^(1/2)+(-647656921614780*Pi-3885941529688680*ln(2)-892778941072368)*3^(1/2)+1121774694224184*Pi+6730648165345104*ln(2)+1546338486359352)*ln(1+2^(1/2))-1570157489678400*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*ln(1+3^(1/2))^2+(4710472469035200*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*arccoth(2^(1/2))+((261692914946400*Pi+1570157489678400*ln(2)+360737167618800)*3^(1/2)-453265424713920*Pi-2719592548283520*ln(2)-624815102694048)*2^(1/2)+(-370089669494160*Pi-2220538016964960*ln(2)-510159394898496)*3^(1/2)+641014110985248*Pi+3846084665911488*ln(2)+883621992205344)*ln(1+3^(1/2))-3532854351776400*(2^(1/2)-3880899/2744210)*(3^(1/2)-70226/40545)*arccoth(2^(1/2))^2+(((-392539372419600*Pi-2355236234517600*ln(2)-541105751428200)*3^(1/2)+679898137070880*Pi+4079388822425280*ln(2)+937222654041072)*2^(1/2)+(555134504241240*Pi+3330807025447440*ln(2)+765239092347744)*3^(1/2)-961521166477872*Pi-5769126998867232*ln(2)-1325432988308016)*arccoth(2^(1/2))+((-392539372419600*ln(2)^2+(-130846457473200*Pi-180368583809400)*ln(2)-10903871456100*Pi^2-30061430634900*Pi+1205966071986042)*3^(1/2)+679898137070880*ln(2)^2+(226632712356960*Pi+312407551347024)*ln(2)+18886059363080*Pi^2+52067925224504*Pi-2088794509063308)*2^(1/2)+(555134504241240*ln(2)^2+(185044834747080*Pi+255079697449248)*ln(2)+15420402895590*Pi^2+42513282908208*Pi-1705493574764511)*3^(1/2)-961521166477872*ln(2)^2+(-320507055492624*Pi-441810996102672)*ln(2)-26708921291052*Pi^2-73635166017112*Pi+2954001523727854)^(1/2)*(5*2^(1/2)-7)*(-5+3*3^(1/2))/((70*2^(1/2)-99)*(15*3^(1/2)-26)*(-2+3^(1/2))*(4*3^(1/2)-7)*(12*2^(1/2)-17)*(2*2^(1/2)-3))

(4)

evalf(stdE)

7031348.295

(5)
 

NULL

Download StdE.mw

I have constructed two Random Variables (X1 and X2).  They work the way I want.  When I take the abs(X1-X2) and ask for the PDF, I get output with Elliptic functions that seems to be in complex space. I read about the Elliptic functions, but got nowhere -- above my head.

I only want the "real" output on the 0 <= t <= 2 domain.  Does anyone know how to parse this output so that I can see a PDF as a function of t?

NULL

with(Statistics)

fx1 := piecewise(t < -1, 0, And(t >= -1, t <= 1), 2*sqrt(-t^2+1)/Pi, t > 1, 0)

fx2 := piecewise(t < -1, 0, And(t >= -1, t <= 1), 2*sqrt(-t^2+1)/Pi, t > 1, 0)

X1 := RandomVariable(Distribution(PDF = unapply(fx1, t)))

X2 := RandomVariable(Distribution(PDF = unapply(fx2, t)))

Dist := abs(X1-X2)

simplify(PDF(Dist, t))

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

(1)
 

NULL

Download AbsDiff.mw

I have one random variable as follows:  r := sqrt(RandomVariable(Uniform(0,1)))

I have another as follows:  x := cos(RandomVariable(Uniform(0, 2*Pi)));

The r variable returns a PDF of 2*t, where t<1; 0 otherwise.  The variable x returns a PDF of 1/(Pi*sqrt(1-t^2)), where |t|<1; 0 otherwise.

For both random variables, I have validated through simulated values.  I am convinced the above PDF functions are correct.  When I multiply r*x in Maple and then ask for a PDF of the resultant random variable, I fail.  Maple seems to not like this.

My question is then as follows...can anyone help me get the PDF of r*x?

1 2 Page 1 of 2