sand15

1576 Reputation

15 Badges

11 years, 121 days

MaplePrimes Activity


These are answers submitted by sand15

You wrote d = ... instead of d := ...

This should answer your question.
Example of use:

Histogram(Sample(d, 10^4))



If you want to compute formally some statistics of d, for instance its PDF or CDF that is another story.

If you directly ask Maple for the PDF/CDF of sqrt(DX2+DY2) you won't get the expected result (I can't with my 2015 version).

To get this PDF (CDF) you have to proceed step by step in order to help Maple:

  • Build the PDF (or CDF) of DX2
     
  • Build the CDF (or PDF) of DY2
     
  • Compute the CDF of DX2+DY2, either from PDF(DX2) and CDF(DY2), or from CDF(DX2) and PDF(DY2) (much simpler than computing directly the PDF of DX2+DY2 from rhe convolution product of  PDF(DX2) by CDF(DY2)).
    Let K(t) the CDF of DX2+DY2.
     
  • Define a random variable Z2 this way
    Z2 := RandomVariable(Distribution(CDF = (t -> K(t))));
  • Finally ask Maple what is the PDF of sqrt(Z2) (Maple 2015 does the job).
    Maybe your Maple version won't do it, if it is the case let me know.


You will get then a (very) lengthy expression for PDF(Z2) that you can compare to a direct sample of sqrt(DX2+DY2).
 

Here is the full code (final result at the end of the worksheet) TriangleEuclidean_sand15.mw

and a plot comparing PDF(Z2) to the histogram of a sample of sqrt(DX2+DY2):


The expression of PDF( sqrt(DX2+DY2) ) is that lengthy that it is almost unreadable (I tried a few ad hoc simplifications, mostly by hand, but the expression of PDF(Z2) remains complex, see the attached file).

To help you make an idea of how complex it is, here is its optimized Python code:

t3 = t ** 2
t4 = 4 * t3
t5 = t4 - 3
t6 = math.sqrt(t5)
t9 = t4 - 1
t10 = math.sqrt(t9)
t11 = 0.1e1 / t10
t13 = t3 < 1
t14 = t3 == 0.1e1 / 0.4e1
t15 = math.sqrt(3)
t17 = t3 + 0.1e1 / 0.3e1
t19 = 2 * t3
t20 = t19 - 1
t21 = 1 / t3
t24 = math.asin(t21 * t20 / 2)
t30 = t3 ** 2
t31 = math.sqrt(t3)
t32 = t31 * t30
t34 = t31 * t3
t36 = t30 * t3
t37 = t31 * t36
t39 = t3 * math.pi
t44 = t30 ** 2
t58 = 0.32e2 / 0.243e3 * (-27 * t24 * t17 * t3 * t10 * t15 + t10 * (t32 * (-0.144e3 / 0.5e1 * t15 + 0.152e3 / 0.5e1) + 72 * t34 - 0.192e3 / 0.35e2 * t37 + t15 * (-0.27e2 / 0.2e1 * t39 + 0.45e2 / 0.32e2 * math.pi + 0.9e1 / 0.10e2) + t44 + 3 * t36 - 0.159e3 / 0.2e1 * t30 - 0.34e2 / 0.5e1 * t3 - 0.29471e5 / 0.8960e4) + (0.402e3 / 0.5e1 * t30 + 0.384e3 / 0.5e1 * t36 - 0.231e3 / 0.10e2 * t3 - 0.9e1 / 0.20e2) * t15) * t11 if t14 else 0
t59 = t58 if t13 else 0
t60 = t6 * t59
t62 = 12 * t3 - 9
t63 = math.sqrt(t62)
t64 = t63 * t10
t67 = t3 - 0.2e1 / 0.3e1
t68 = 0 if t13 else 1
t70 = 0 if t3 < 0.3e1 / 0.4e1 else 1
t71 = t68 - t70
t72 = t71 * t67
t73 = t6 / 2
t74 = -t73 < 0
t75 = 0 if t74 else 1
t77 = t73 < 0
t78 = 0 if t77 else 1
t80 = t3 + 0.4e1 / 0.3e1
t82 = -t80 * t68 + t75 * t72 - t78 * t72
t83 = -0.3e1 / 0.4e1 + t3
t84 = t83 * t82
t85 = t10 * t3
t86 = t19 - 3
t89 = math.asin(t21 * t86 / 2)
t90 = t89 * t85
t94 = 0 if t3 < 0.1e1 / 0.4e1 else 1
t95 = t68 - t94
t96 = t15 * t95
t97 = t6 * t96
t98 = t24 * t17
t99 = t98 * t85
t105 = t15 * math.pi
t106 = t105 * t3 * t67
t108 = 9 * t36
t111 = 0.48e2 / 0.5e1 * t32 * t15 - 0.9e1 / 0.4e1 * t106 + t44 - t108 - 0.81e2 / 0.8e1 * t30 + 0.27e2 / 0.16e2 * t3 - 0.243e3 / 0.1280e4
t118 = t6 * t111 - 0.128e3 / 0.35e2 * t83 * (t36 - 0.26e2 / 0.3e1 * t30 + 0.433e3 / 0.128e3 * t3 + 0.333e3 / 0.256e3)
t119 = t118 * t71
t120 = t75 * t10
t122 = t78 * t10
t124 = t15 * t68
t125 = t6 * t124
t126 = t3 * t80
t127 = t3 - 2
t129 = math.asin(t21 * t127)
t131 = t129 * t10 * t126
t134 = t3 + 3
t135 = t3 - 1
t136 = math.sqrt(t135)
t137 = t136 * t135
t138 = t137 * t134
t141 = t135 ** 2
t142 = t136 * t141
t157 = 0.80e2 / 0.3e1 * t32 + 12 * t34 - 0.128e3 / 0.35e2 * t37 + t15 * (-0.3e1 / 0.2e1 * t39 - 0.9e1 / 0.2e1 * t30 * math.pi) - 0.53e2 / 0.2e1 * t30 + 0.4e1 / 0.3e1 * t44 - 8 * t36 + 0.17e2 / 0.448e3 - 0.34e2 / 0.15e2 * t3
t159 = 0.48e2 / 0.5e1 * t15
t160 = -0.152e3 / 0.15e2 + t159
t162 = t68 * t160 - t159 - 0.248e3 / 0.15e2
t165 = -24 * t68 + 12
t167 = t68 + 1
t169 = t3 - 0.10e2 / 0.3e1
t170 = t136 * t169
t178 = math.pi * t15 * t3
t179 = 0.9e1 / 0.2e1 * t178
t180 = -0.2e1 / 0.3e1 * t44 + 0.439e3 / 0.8e1 * t30 + 0.7493e4 / 0.240e3 * t3 - 2 * t36 - 0.2593e4 / 0.1792e4 + t179
t183 = -8 * t138 * t124 + 0.24e2 / 0.5e1 * t142 * t124 + t94 * t157 + t32 * t162 + t34 * t165 + 0.64e2 / 0.35e2 * t37 * t167 + 9 * t170 * t124 + t68 * t180 + 0.9e1 / 0.2e1 * t106 - t44 + t108
t185 = t83 * t68
t188 = t36 - 0.31e2 / 0.6e1 * t30 - 0.1919e4 / 0.128e3 * t3 - 0.423e3 / 0.256e3
t191 = t6 * t183 + 0.128e3 / 0.35e2 * t188 * t185
t194 = t30 + 0.83e2 / 0.64e2 * t3 + 0.3e1 / 0.128e3
t195 = t194 * t95
t197 = (-0.1e1 / 0.4e1 + t3) * t15
t198 = t6 * t197
t201 = 9 * t99 * t97 + t120 * t119 - t122 * t119 + 0.9e1 / 0.2e1 * t131 * t125 + t10 * t191 - 0.128e3 / 0.5e1 * t198 * t195
t203 = -0.81e2 / 0.32e2 * t64 * t60 - 54 * t90 * t84 + t63 * t201
t204 = 0.1e1 / t63
t205 = t204 * t203
t208 = 0.1e1 / t6
t214 = t11 * t208
t215 = undefined if t14 else 0
t230 = undefined if t3 == 1 else 0
t232 = undefined if t3 == 0.3e1 / 0.4e1 else 0
t233 = t230 - t232
t234 = t233 * t67
t236 = 0 if t74 else 0
t240 = 0 if t77 else 0
t248 = t89 * t10
t253 = t11 * t3
t258 = 1 / t30
t264 = math.sqrt(t258 * t62)
t269 = t230 - t215
t280 = t15 * t230
t342 = t15 * math.pi * t67
t344 = 4 * t36
t345 = 27 * t30
t346 = 0.64e2 / 0.35e2 * t37 * t230 + 0.32e2 / 0.5e1 * t32 * t167 + 9 * t170 * t280 + 9 * t136 * t124 + 0.9e1 / 0.2e1 / t136 * t169 * t124 + t68 * (-0.8e1 / 0.3e1 * t36 + 0.439e3 / 0.4e1 * t3 + 0.7493e4 / 0.240e3 - 6 * t30 + 0.9e1 / 0.2e1 * t105) + t230 * t180 + t179 + 0.9e1 / 0.2e1 * t342 - t344 + t345
t356 = 3 * t30
t365 = t118 * t233
t384 = (t6 * (24 * t34 * t15 - 0.9e1 / 0.4e1 * t178 - 0.9e1 / 0.4e1 * t342 + t344 - t345 - 0.81e2 / 0.4e1 * t3 + 0.27e2 / 0.16e2) + 2 * t208 * t111 - 0.128e3 / 0.35e2 * t83 * (t356 - 0.52e2 / 0.3e1 * t3 + 0.433e3 / 0.128e3) - 0.128e3 / 0.35e2 * t36 + 0.3328e4 / 0.105e3 * t30 - 0.433e3 / 0.35e2 * t3 - 0.333e3 / 0.70e2) * t71
t392 = 9 * t99 * t6 * t15 * t269 + 18 * t99 * t208 * t96 + 18 * t98 * t253 * t97 + 0.9e1 / 0.2e1 * t131 * t6 * t280 + 9 * t131 * t208 * t124 + 9 * t129 * t11 * t126 * t125 + t10 * (t6 * (-8 * t138 * t280 + 4 * t137 * t124 - 12 * t136 * t134 * t124 + 0.24e2 / 0.5e1 * t142 * t280 + t94 * (0.200e3 / 0.3e1 * t34 + 18 * t31 - 0.64e2 / 0.5e1 * t32 + t15 * (-0.3e1 / 0.2e1 * math.pi - 9 * t39) - 53 * t3 + 0.16e2 / 0.3e1 * t36 - 24 * t30 - 0.34e2 / 0.15e2) + t215 * t157 + t32 * t230 * t160 + 0.5e1 / 0.2e1 * t34 * t162 - 24 * t34 * t230 + 0.3e1 / 0.2e1 * t31 * t165 + t346) + 2 * t208 * t183 + 0.128e3 / 0.35e2 * t188 * t83 * t230 + 0.128e3 / 0.35e2 * t188 * t68 + 0.128e3 / 0.35e2 * (t356 - 0.31e2 / 0.3e1 * t3 - 0.1919e4 / 0.128e3) * t185) + 2 * t11 * t191 + t120 * t365 + t120 * t384 + t236 * t10 * t119 + 2 * t75 * t11 * t119 - t122 * t365
t407 = math.sqrt(4) * math.sqrt(t258 * t135)
t415 = t17 * t10
t420 = math.sqrt(t258 * t9)
t449 = -t122 * t384 - t240 * t10 * t119 - 2 * t78 * t11 * t119 - 0.128e3 / 0.5e1 * t6 * t15 * t195 + 0.9e1 / 0.2e1 / t407 * (-t258 * t127 + t21) * t85 * t80 * t6 * t124 + 18 / t420 * (t21 - t258 * t20 / 2) * t415 * t3 * t6 * t96 + 9 * t24 * t415 * t97 + 9 * t24 * t85 * t97 + 0.9e1 / 0.2e1 * t129 * t85 * t125 + 0.9e1 / 0.2e1 * t129 * t10 * t80 * t125 - 0.128e3 / 0.5e1 * t198 * t194 * t269 - 0.128e3 / 0.5e1 * t198 * (t19 + 0.83e2 / 0.64e2) * t95 - 0.256e3 / 0.5e1 * t208 * t197 * t195
t454 = -0.81e2 / 0.32e2 * t64 * t6 * t215 - 0.81e2 / 0.16e2 * t64 * t208 * t59 - 0.81e2 / 0.16e2 * t63 * t11 * t60 - 0.243e3 / 0.16e2 * t204 * t10 * t60 - 54 * t90 * t83 * (-t80 * t230 + t75 * t234 - t78 * t234 + t236 * t72 - t240 * t72 + t75 * t71 - t78 * t71 - t68) - 54 * t248 * t3 * t82 - 54 * t248 * t84 - 108 * t89 * t253 * t84 - 108 / t264 * (t21 - t258 * t86 / 2) * t10 * t3 * t84 + t63 * (t392 + t449) + 6 * t204 * t201
t466 = 0 if t <= 0 else 2 * (0.64e2 / 0.81e2 * t205 * t11 / t6 / t5 + 0.64e2 / 0.81e2 * t205 / t10 / t9 * t208 - 0.32e2 / 0.81e2 * t204 * t454 * t214 + 0.64e2 / 0.27e2 / t63 / t62 * t203 * t214) * t if 0 < t else 0


Or enlarge this image

Customize as you want

NULL

 

NULL

``

restart

with(Statistics):

X := RandomVariable(Normal(0, 1)):

f := unapply(PDF(X, x), x):

g := unapply(CDF(X, x), x):

``

h := proc(t)
if t::numeric then
  plots:-display(
    plot(f(x),  x = -5 .. t, gridlines = true, filled=true, color="Chartreuse")
    ,
    plot(g(x), x = -5..t, gridlines=false, color=blue)
    ,
    plot([[t, 0], [t, g(t)]], linestyle=3, color="DarkGray")
    , title=typeset(Int('f'(x), x=-5..evalf[2](t)) = evalf[4](g(t)))
  )
else
  'procname(_passed)'
end if:
end proc:

plots:-animate(h, [t], t=-5.0 .. 5.0, paraminfo=false)

 

 

``

``

NULL

 

 

NULL

``

 

Download Q_CDF_PDF_sand15.mw


A variant using dualaxisplot: (look how to suppress the info about the parameter value [Maple 2015]) 

h := proc(t)
if t::numeric then
  plots:-dualaxisplot(
    plot(f(x),  x = -5 .. t, gridlines = true, filled=true, color="Chartreuse", legend=:-PDF('X'))
    ,
    plots:-display(
      plot(g(x), x = -5..t, gridlines=false, color=blue, legend=:-CDF('X'))
      ,
      plot([[t, 0], [t, g(t)]], linestyle=3, color="DarkGray")
      , axis[2]=[color=blue]
    )
    , title=typeset(Int('f'(x), x=-5..evalf[2](t)) = evalf[4](g(t)))
  )
else
  'procname(_passed)'
end if:
end proc:

plots:-animate(h, [t], t=-5.0 .. 5.0, paraminfo=false)


With restrictions  y > 1 and _C1::real one gets this solution
Asolution.mw


Given we assumed y > 1, the solution is A_solution_gen[1].

I don't think that using the Horner representation is judicious.
What we generally want to know are the coefficients (that is the effects, and often the statistical significances) of xn terms, an information the Horner factorization does not provide directly, which makes it useless.

Note: be careful in using statistics returned by the output=solutionmodule option in the case you use PolynomialFit or ExponentialFit because some are wrong (I believe I posted something about that, either under this login or under the login @mmcdara, I don't remembe exactly).

Nevertheless "leastsquaresfunction" and "residuals" output, and a few more, are correct.

restart;

X   := [seq(0..7,0.5)]:
phi := x -> sin(x):
Y   := phi~(X):


Only necessary outputs

fit, residuals := op(Statistics:-PolynomialFit(10, X, Y, x, output=["leastsquaresfunction", "residuals"])):

Predictor := unapply(fit, x):

plot(
  [phi, Predictor](x)
  , x=(min..max)(X)
  , color=[blue, red]
  , transparency=[0, 0.7]
  , thickness=[2, 11]
  , legend=[typeset('phi'(x)), typeset('Predictor'(x))]
);

plot(
  [ seq([X[i], residuals[i]], i=1..numelems(X)) ]
)

 

 


You might be interested to know a lot of other outpurts are avaliable

# Yo get all of them:

P := Statistics:-PolynomialFit(10, X, Y, x, output=solutionmodule):

# To see them

P:-Results():

# Extract residuals alone

residuals := P:-Results("residuals"):

dataplot(residuals);

 

Predictor := unapply(P:-Results("leastsquaresfunction"), x):

plot(
  [phi, Predictor](x)
  , x=(min..max)(X)
  , color=[blue, red]
  , transparency=[0, 0.7]
  , thickness=[2, 11]
  , legend=[typeset('phi'(x)), typeset('Predictor'(x))]
)

 
 

 

Download regr_fun_sand15.mw

I suggest you not to load simultanously Student[VectorCalculus] and LinearAlgebra to avoid conclicts (there are some in Maple 2015)
Once done (Maple 2015) I can't see any error.

Do you mean instead that the result Maple produces does not suit you?
If it is so maybe you should verify twice what you did by hand...

Question_1_sand15.mw

EDITED

Because the formal integration of f,  not even to speak of the convolution product f  f, does not seem possible, I propose you two alternatives:

  1. A pointwise approximation of f followed by a discrete approximation (in fact a Riemann sum) of the integral which represents f  f
  2. A spline reconstruction of f  followed by an exact calculation of  f  f

Both approaches give equivalent results provided enough points are considered in both of them (I used 41 points in approach (1) and 21 in approach (2) but those numbers were chosen arbitrarily and it is likely that already good results could be obtained with smaller numbers).

Details are given in the attached worksheet, the animation below is the ice on the cake (if you like iced cakes of course).

Here is an updated version of my initial worksheet which contains a numerical estimation of  f  f  using evalf/Int (cpu times seems prohibitive).

lternatives_to_formal_integration_updated.mw

add option echoexpression=false

Maple 2015 help:
help(Explore) > section "See Also" > Explore example worksheet > Exploring a mathematical expression > example 2

Your question can be highly simplified if you ask if the arguments of the log functions in G0 and G01 are equal.

NULL

restart

with(plots):

H0 := -S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2;

-S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2

(1)

NULL

Z0 := exp(-beta*H0);

exp(-beta*(-S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2))

(2)

Z0 := add(Z0, S1 = [-2, -1, 0, 1, 2]):

# Your problem will be simpler to solve if you consider Z0 instead of G0

Z01 := (2*exp(4*beta*eta1)*cosh(2*beta*gamma1)+2*exp(beta*eta1)*cosh(beta*gamma1)+1)^((1/2)*N)*(2*exp(4*beta*eta2)*cosh(2*beta*gamma2)+2*exp(beta*eta2)*cosh(beta*gamma2)+1)^((1/2)*N)

(2*exp(4*beta*eta1)*cosh(2*beta*gamma1)+2*exp(beta*eta1)*cosh(beta*gamma1)+1)^((1/2)*N)*(2*exp(4*beta*eta2)*cosh(2*beta*gamma2)+2*exp(beta*eta2)*cosh(beta*gamma2)+1)^((1/2)*N)

(3)

# convert(Z01, exp) replaces hyperbolic trigs by exponentials

T := simplify(convert(Z01, exp));

(exp(2*beta*(2*eta1+gamma1))+exp(beta*(eta1+gamma1))+exp(2*beta*(2*eta1-gamma1))+exp(beta*(eta1-gamma1))+1)^((1/2)*N)*(exp(2*beta*(2*eta2+gamma2))+exp(beta*(eta2+gamma2))+exp(2*beta*(2*eta2-gamma2))+exp(beta*(eta2-gamma2))+1)^((1/2)*N)

(4)

# replace T = AN/2 * BN/2 by (A*B)N/2
# and discard the exponent already in common with G0

T0 := op([1, 1], T)* op([2, 1], T)

(exp(2*beta*(2*eta1+gamma1))+exp(beta*(eta1+gamma1))+exp(2*beta*(2*eta1-gamma1))+exp(beta*(eta1-gamma1))+1)*(exp(2*beta*(2*eta2+gamma2))+exp(beta*(eta2+gamma2))+exp(2*beta*(2*eta2-gamma2))+exp(beta*(eta2-gamma2))+1)

(5)

# the question becomes: Is T0 equal to Z0?

simplify(T0 - Z0)

0

(6)
 

 

Download TESTE_MAPLE_sand15.mw

You cannot simply use the commands solve(cot(y)=A, y) and solve(cot(y)=-A, y) to get the inverse function of cot function: arccot is made of an infinity of branches and the expression of cot(-1) depends on the branch you consider.

Using solve(cot(y)=A, y) and solve(cot(y)=-A, y) gives you the impression that cot(-1) is made of two branches, one for A > 0, the other for A < 0.
(BTW how can Maple distinguish these two cases without explicitely using assumptions 

solve(cot(y) =  A, y, useassumptions=true) assuming A > 0;
solve(cot(y) = -A, y, useassumptions=true) assuming A > 0;

?)

Use solve(cot(y)=A, y, allsolutions) instead to build cot(-1) on the branch of interest.

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

branches(arccot)

 

# The inverse function of 'cot' depends on the branch you consider

`#msup(mo("cot"),mo("-1"))` := solve(cot(y)=A, y, allsolutions);

about(_Z1);

arccot(A)+Pi*_Z1

 

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

 

# Main branch

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=0)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)

(2)

# branch (-Pi, 0)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-1)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)-Pi

(3)

# branch (-2*Pi, -Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-2)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)-2*Pi

(4)

# branch (Pi, 2*Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=1)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)+Pi

(5)
 

 

Download arccot.mw

You say that other CAS answer differently? I advice you to ree you to read carefully how the inverse trig functions are (by default) defined and how to get them on any branch. 
Maybe Maple has its own strategy about what the "default" solve(cot(y)=A, y) must return,  but I doubt using Maple, MMA or SageMath give another answer if they are used cautiously.

@FDS 

A question: your second graph shows your material experiences some kind of strain-hardening and does not have an elastic regime.
As there is not linear part in the charge and discharge branches of the "stress-strain curve" (between quotes, see the end of the next paragraph), how do you define the slopes you are looking for?

By the way I'm quite surprised by the unusual concavity of the loading branches of the stress-strain curve: Shouldn't have an opposite concavity to those of the unloading branches?  What is the constitutive law which governs the behaviour of your material  (it looks like some kind of dilatant fluid)?
Would the horizonal axis represent the shear rate instead of the strain?

Please have a look to Stress_Strain_Curve.mw and tell me in what direction go.
(The reason I fit charge/discharge branches with polynomial models comes from my intuition you are working with a shear-thickening fluid, fluids whose behaviours are usually represented by powerlaw models... but maybe I am on the wrong track).

With "my" Maple 2015 (but this would work for more recent versions too), and using MathML:

restart

# MathML hexa code of symbol x

`#mo("&#x2243;")`

`#mo("&#x2243;")`

(1)

# Use this symbol as a binary operator.. version 1

`&#x2243;`(Pi, 3.14);  # not very pretty

`&#x2243;`(Pi, 3.14)

(2)

# Version 2;: define operator 'ae' (almost equal)
# see acer's comment
`print/ae` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("&InvisibleTimes;"),
        mo("&#x2243;",fontweight="bold"), mo("&InvisibleTimes;"),
        Typeset(y));
end proc:

ae(Pi, 3.14);  # prettier

ae(Pi, 3.14)

(3)

`print/ra` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("&InvisibleTimes;"),
        mo("&#x2192;",fontweight="bold"), mo("&InvisibleTimes;"),
        Typeset(y));
end proc:

plot(cos(x), x=0..2, title = typeset(ra(cos(x)=0, ae(x, 1.6))), titlefont=[Times, bold, 16])

 
 

 

Download AlmostEqual.mw

Understandable, not "elegant" and far from being a one-liner command (unless you wrapped this intoa procedure... if course).

ex:=(a+b*c)^2

(b*c+a)^2

(1)

shorten := table([_Inert_POWER="`^`", _Inert_SUM="`+`", _Inert_PROD="`*`", "("="[", ")"="]"]):

inertex := ToInert(ex):

shortstr := convert( eval(inertex, {_Inert_NAME = (u -> parse(u)), _Inert_INTPOS = (u -> u)}), string):

shortstr := copy(shortstr):
for i in [indices(shorten, nolist)] do
  shortstr := StringTools:-SubstituteAll(shortstr, i, short[i])
end do:

cat("[", shortstr, "]")

"[`^`[`+`[`*`[b,c],a],2]]"

(2)
 

 

Download proposal.mw

UPDATED


The bug you observed is more subtle than what you think and does not reduce to a regression from one version to another.

Indeed, even if some Maple versions prior to Maple 2024 may return a no-error result (see @acer for instance), this does not guarantee this result is correct.
For instance Maple 2015 generates no error in calculating the pdf of sin(theta) but the resul it provides is wrong.
 

This latter claim is easy to diagnose:

  • If you know a few of Statistics it is obvious that both cos(theta) and  sin(theta) must have the same PDF... and clearly Maple does not provide identical results.
    Note that Maple 2015 computes correctly the pdf of cos(theta), but this is only by chance and Maple 2015 is not that lucky in computing the one of sin(theta)).
     
  • If you know a little bit more in Statistics, sin(theta) has, by definition, an arsine distribution , and if you know what thecorresponding pdf looks like, you can easily verify that Maple 2015's result is wrong.

The fact Maple 2015 (and likely all other Maple versions) returns (may return?) a wrong result is the consequence of a flawed mechanism for calculating PDF(𝜑(theta)) when 𝜑 is not a one-to-one map over the whole support of theta (see my comment ahead).
 

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

with(Statistics)

theta := RandomVariable(Uniform(0, 2*Pi))

_R

(2)

phi := arccos(-1+2*RandomVariable(Uniform(0, 1))); phi := arccos(RandomVariable(Uniform(-1, 1)))

arccos(-1+2*_R0)

 

arccos(_R1)

(3)

r := RandomVariable(Uniform(0, 1))^(1/3)

_R2^(1/3)

(4)

SinPhi := sin(phi)

(-_R1^2+1)^(1/2)

(5)

f__sp := unapply(PDF(SinPhi, t), t); f__sp(t)

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

(6)

SinTheta := sin(theta)

sin(_R)

(7)

f__st := unapply(PDF(SinTheta, t), t); f__st(t)

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

(8)

CosPhi := cos(phi)

_R1

(9)

f__cp := unapply(PDF(CosPhi, t), t); f__cp(t)

piecewise(t < -1, 0, t < 1, 1/2, 0)

(10)

CosTheta := cos(theta)

cos(_R)

(11)

f__ct := unapply(PDF(CosTheta, t), t); f__ct(t)

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

(12)

 

Supp := `~`[Support]([SinPhi, CosPhi, SinTheta, CosTheta], output = 'range')

[0 .. 1, -1 .. 1, 0 .. 1, -1 .. 1]

(13)

plots:-display(plot(f__sp(t), t = Supp[1], color = red, legend = typeset(sin('Phi'))), plot(f__cp(t), t = Supp[2], color = red, linestyle = 3, legend = typeset(cos('Phi'))), plot(f__st(t), t = Supp[3], color = blue, legend = typeset(sin('Theta'))), plot(f__ct(t), t = Supp[4], color = blue, linestyle = 3, legend = typeset(cos('Theta'))))

 

DocumentTools:-Tabulate(
  [
    plots:-display(
      Histogram(Sample(SinPhi  , 10^4), title=typeset(sin('Phi')), style='polygon', transparency=0.8),
      plot(f__sp(t), t = Supp[1], color = red)
    ),
    plots:-display(
      Histogram(Sample(CosPhi  , 10^4), title=typeset(cos('Phi')), style='polygon', transparency=0.8),    
      plot(`f__cp`(t), t=Supp[2], color=red)
    ),
    plots:-display(
      Histogram(Sample(SinTheta, 10^4), title=typeset(sin('Theta')), style='polygon', transparency=0.8),    
      plot(`f__st`(t), t=Supp[3], color=red),
      plots:-textplot([0, 10, "WRONG PDF"], font=[Tahoma, bold, 20], color=red)
    ),
    plots:-display(
      Histogram(Sample(CosTheta, 10^4), title=typeset(cos('Theta')), style='polygon', transparency=0.8),    
      plot(`f__ct`(t), t=Supp[4], color=red)
    )
  ]
)

# Correct computation of PDF(SinTheta)

Tneg := RandomVariable(Uniform(-Pi, 0)):
Tpos := RandomVariable(Uniform(0, Pi)):

fneg := unapply(PDF(sin(Tneg), t), t);
fpos := unapply(PDF(sin(Tpos), t), t);

plots:-display(
  plot((1/2)*~[fneg, fpos](t), t=-1..1, color=[red, blue]),
  Histogram(Sample(SinTheta, 10^5), title=typeset(sin('Phi')), style='polygon', transparency=0.8, minbins=100)
  , gridlines=true
)

proc (t) options operator, arrow; piecewise(t <= -1, 0, t < 0, 2/(Pi*(-t^2+1)^(1/2)), 0 <= t, 0) end proc

 

proc (t) options operator, arrow; piecewise(t < 0, 0, t < 1, 2/(Pi*(-t^2+1)^(1/2)), 1 <= t, 0) end proc

 

 

# Or, directly from symmetry considerations

Tpos := RandomVariable(Uniform(0, Pi)):
fpos := unapply(PDF(sin(Tpos), t), t);

f__SinTheta := 1/2 * (fpos(-t) + fpos(t));

plots:-display(
  plot(f__SinTheta, t=-1..1, color=[red, blue]),
  Histogram(Sample(SinTheta, 10^5), title=typeset(sin('Phi')), style='polygon', transparency=0.8, minbins=100)
  , gridlines=true
)

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

 

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

 

 


Why doesn't Maple compute correctly the PDF of  SinTheta?

d

# Let X some random variable and f : X --> Y a C1 function of X.
# Let Supp(X) the support of X.
#
# In the simplest case where phi has a unique inverse for each x Supp(X), then
# an infinitesimal neihghborhood V(x) of diameter d(x) of x has the image
# d(y) = |f'(x)|d(x) (d(x) assumed > 0 and f'(x) denotes the derivative of f(x)).
#
# The probability that an outcome of X belongs to the interval [x-d(x), x+d(x)]
# has the value fX(x)d(x), where fX denotes the PDF of X.
# This can be rewritten fX(x)d(x) = fX(x)/|f'(x)|d(y).
# So the PDF of Y is simply fY(y) = fX(x)/|f'(x)|,
#
# Introducing the reciprocal function y = f(-1) (which exists given the assumptions
# above) one has x=y(y) and thus
#         fY(y) = fX(y(y))|y'(y)|
#               = ( fX(y(y)) )'
#               = ( FX(y(y)) )'   (FX denotes the CDF of X)
#               = ( FY(y) )'      (= derivative of the CDF of Y = f(x))
#
#
#
# For your problem it is worth noting that neither the sine nor cosine f functions have
# a unique inverse on the support of tandom variable theta.
#
# When f does not have a unique inverse one must proceed more carefully:
#    (1) split the whole support WX of X into sub-intervals wn such that the restriction
#        of f has a unique inverse over each wn,
#    (2) compute the contributions Cn = (FX+y)'(y) over each wn and sum them to get the
#        final result.
#
# Here is a graphical illustration for the sine function.
# The probability that Y belongs to the interval [0.45, 0.55] is equal to the probability
# that X belongs to interval I1 OR that X belongs to interval I2.


I1 := sort([fsolve(sin(x)=0.45, x=0..Pi/2.), fsolve(sin(x)=0.55, x=0..Pi/2.)]);
I2 := sort([fsolve(sin(x)=0.45, x=Pi/2...Pi), fsolve(sin(x)=0.55, x=Pi/2...Pi)]);


plots:-display(
  plot(sin(x), x=0..2*Pi, color=blue, labels=["X", "Y"], labelfont=[Tahoma, bold, 12])
  , plot(0.5, x=0..2*Pi, color=black)
  , plot(0.55, x=0..2*Pi, color=black, linestyle=3)
  , plot(0.45, x=0..2*Pi, color=black, linestyle=3)
  , plottools:-rectangle([-0.05, 0.45], [0.05, 0.55], color=black)
  , plots:-textplot([2*Pi, 0.55, y+dy], align={left, above})
  , plots:-textplot([2*Pi, 0.45, y-dy], align={left, below})
  , plot([[I1[1], 0.45], [I1[1], 0]], color=black, linestyle=3)
  , plot([[I1[2], 0.55], [I1[2], 0]], color=black, linestyle=3)
  , plottools:-rectangle([I1[1], -0.05], [I1[2], 0.05], color=black)
  , plot([[I2[1], 0.45], [I2[1], 0]], color=black, linestyle=3)
  , plot([[I2[2], 0.55], [I2[2], 0]], color=black, linestyle=3)
  , plottools:-rectangle([I2[1], -0.05], [I2[2], 0.05], color=black)
  , plots:-textplot([I1[1], 0.1, x-dx], align={left})
  , plots:-textplot([I1[2], 0.1, x+dx], align={right})
  , plots:-textplot([I2[1], 0.1, x-dx], align={left})
  , plots:-textplot([I2[2], 0.1, x+dx], align={right})
  , size=[800, 400]
)

[.4667653390, .5823642379]

 

[2.559228416, 2.674827315]

 

 

# A domain where phi (aka 'sin') has a unique inverse:

T := RandomVariable(Uniform(-Pi/2, Pi/2)):

PDF(sin(T), t);
plot(%, t=-1..1);
 

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

 

 

# But on the same domain phi (aka 'cos') does not have a unique inverse... so a wrong result

PDF(cos(T), t);
plot(%, t=-1..1);

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

 

 

# To get a correct result one must translate T in such a way that Ttranslated has a unique
# 'cos' inverse once the translation is done.
# Here, simply:

PDF(cos(T + Pi/2), t);
plot(%, t=-1..1);

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

 

 

# It should be clear to you that cos(theta) and sin(theta) (once computed correctly)
# must have the same distribution.

 

 

Download Basics_Maple2015.mw

f(x) is singular at x=0 meaning its graph C is not continuous.
It's clear hat max(f(x) = +oo and min(f(x)) = -oo
So there is obviously no circle which answers your problem.

So i suppose that you saying "where the curve attains its maximum and minimum" has to be interpreted as "where diff(f(x), x) = 0 and diff(f(x), x, x) <> 0", am I right?
If it is so, because f(x) is antisymmetric if the point (a, f(a)), a being assumed > 0, is the point where f'(a)=0, then we have f'(-a)=0 too.
If there were a circle tangent to C at points (a, f(a)) and (-a, -f(a)), then, by symmetry condition, its center would br the point (0, 0).
So this circle could be tangenf to the curve C at points (a, f(a)) and (-a, -f(a)) if and only if a=0... which is impossible because f(0) is undefined.
So there is no circle solution of your problem.

@acer has provided the correct answer  but I am not pleased with it because it is a workaround to a Maple's weakness (unfortunately a workaround you must use if x and y have less "sympathetic" distrbutions).

As a rule I personnaly prefer to try and apply firstly a log transformation (logz = log(x) + log(y)) and compute next the pdf of z=exp(logz).
This trick does not work for any couple of distributions, but when it does, and it happens to be case here, it avoids using artificial (because a pdf is defined over the whole real line and there is no real justification to limit it to some interval) assumptions.

restart:

with(Statistics):

x := RandomVariable(Uniform(1, 2)):
y := RandomVariable(Uniform(2, 3)):

LX := log(x):
LY := log(y):

LZ := LX + LY:

PDF(LZ, zeta)

piecewise(Zeta <= ln(2), 0, Zeta <= ln(3), -exp(Zeta)*(-Zeta+ln(2)), Zeta <= 2*ln(2), exp(Zeta)*(ln(3)-ln(2)), Zeta <= ln(2)+ln(3), exp(Zeta)*(-Zeta+ln(2)+ln(3)), ln(2)+ln(3) < Zeta, 0)

(1)

Z := exp(LZ):

PDF(Z, z)

piecewise(z <= 2, 0, z <= 3, ln(z)-ln(2), z <= 4, ln(3)-ln(2), z <= 6, ln(2)+ln(3)-ln(z), 6 < z, 0)

(2)
 

 

Download MultiplyPDFFunctions_sand15.mw

1 2 3 4 5 6 7 Last Page 1 of 11