sand15

1479 Reputation

15 Badges

11 years, 90 days

MaplePrimes Activity


These are replies submitted by sand15

@FDS 

The updated worksheet  Stress_Strain_Curve_2.mw
And an illustraton of what you can do with it in the spirit of what you show in your pptx

About polynomial regression... there is really no mistery behind it.
The first question is: Are you familiar with the "multiple" linear regression in the case of Q regressors X1, ...XQ?
If it is so, "regressing" a dense polynomial of degree Q-1 with a single indeterminate Z, is nothing but using the "multiple" linear regression while setting Xn = Z(n-1), n=1..Q+1.
To give an example, suppose you want to fit a dense quadratic polynomial a+bZ+cZ2: simply introduce auxiliary variables X1 = 1, X2 =Z, X3 =Z2 and regress your dependent quantity no longer on  a+bZ+cZ2:  but on aX1+bX2+cX3...and that's it!

If you are not familiar with "multiple" regression I advice you to read any elementary textbook/course/mooc... on the subject.

Beyond that an important question often appears: what is the degree of the polynomial to use in order that the fitted polynomial fulfills two contradictory requirements:

  1. Fidelity to the data: otherwise said the fitted polynomial must be "close" to the data or, mathematically speaking, the residual sum of squares between the data and their restitution by the fitted polynomial must be small.
    Obviously, as a dense polynomial of degree D-1 depends on D parameters, a set od data of size D sill be exactly fitted by this dense polynomial.
    We use to say here that the model is over-parameterized
  2. Stability (or robustness, or generalization capability...): let us take this previous sample of size D and this same polynomial of degree D-1.
    This polynomial may have D-1 real roots and if, unhappily, some of them are located within the range of the regressor (Z), then the fitted polynomial will present undesired overshoots in between consecutive values of the regressor.
    In fact it is extremely likely that the fitted polynomial will behave this way.
    If you want to avoid this phenomenon use a polynomial of degree 0 or 1.

So, to satisfy the fidelity to the data requirement you must take a large degree, close to the sample size, but if you are more concern by "stability" or "smooth fit", then you must take low degree polynomials.
There exist a lot of litterature about strategies which balance these two contradictory requirements. I already gave you a few Wiki references on the subject.

What I wrote in my worksheet is that each branch is a sample of size about one thousand.
So a dense polynomial of degree 5, for instance, is very far ffrom being an  over-parameterized polynomial. This balance stuff I spoke abour above matters only when the number of parameter of the model (not necessarily a polynomial one) is of the order of the sample size, let us say 20% or 25% of this latter.
You would be right considering model selection (here degree selection) if you were investigating polynomials of degree 200 for example. But it is very farfrom being the case.

As I wrote in this last worksheet a mor important criterion is that two polynomial models intersect within L5 and EP5 ranges if you do not want to create unealistic approximated cycles.

@FDS 

Feel free to tell me if this goes to the right direction
Stress_Strain_Curve_2.mw

I understand the notion of cyclic loading-unloading , but you want to assess the "surface area" of each cycle, and so you need to define more precisely what characterizes a cycle in terms of geometry: where does it start (likely the intersection point between a charge path and the next/previous discharge path) and where it ends.
So I propose you my interpretation for the first cycle only... I am not going to go further if I am not on the right track in me doing this.

My Maple version is too old to manage workbooks. Could you simply upload the stress-strain data?

@janhardo 

Beautiful, I vote up

Error messages are pretty clear: you have to specify an approximate solution.

For more detais execute help(dsolve,numeric_bvp,advanced) and help(dsolve,numeric,BVP) in a worksheet and search for the word approxsoln.

Defining an approximate solution can be quite tricky and the best way to do this is to construct a simplified version your problem using a simplified physics (do a dimention analysis, identify second order effects, discard them), to solve it, and to use this solution as approxsoln (which can be either a procedure or an array).
Your are likely the better placed to do this because it is you who knows the physics behind your equations.

For instance, what is order of magnitude of
compared to other terms? Can you simplify it?

Or can you approximate BCs like?
Can you simplify them from a physical point of view?
As a rule a BC on the first derivative which nonlinearly depends on the second derivative makes a BC problem complicated to solve.

My advice is: begin to do this by yourself and come back to us if you get Maple difficulties to sove your problem.

@JoyDivisionMan 

By the way, the attached worksheet explains in more details how to get the pdf of Y = cos(X) where X  is a uniform random variable with support (-𝜋, +𝜋).
Two approaches are used: a direct construction of PDF(Y) and an indirect one by derivation of CDF(Y).

To compute the cdf of, let's say, X = sin(X) (which is obviously the same as Y) without doing all the stuff contained in the attached file (solve/allsolutions, branch identification, aggregation of partial pdf/cdf over the different branches, ...) it is much simpler to observe that sin(X) = cos(𝜋/2X) then define an auxiliary random variable A = 𝜋/2and replace everywhere it occurs  fX(something) by fA(some other thing) which is equal to fX(𝜋/2something)

restart

with(Statistics):

X := RandomVariable(Uniform(-Pi, Pi)):

Y := cos(X):

f__Y := CodeTools:-Usage( unapply(PDF(Y, y), y) )

memory used=17.78MiB, alloc change=34.00MiB, cpu time=273.00ms, real time=2.38s, gc time=0ns

 

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

(1)

restart

with(Statistics):

X := RandomVariable(Uniform(-Pi, Pi)):

Y := cos(X):

F__Y := CodeTools:-Usage( unapply(CDF(Y, y), y) )

memory used=41.23MiB, alloc change=70.01MiB, cpu time=618.00ms, real time=2.02s, gc time=56.15ms

 

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

(2)

f__Y := unapply(diff(F__Y(y), y) , y)

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

(3)


FIRST WAY: BUILD DIRECTLY PDF(Y)

plot(cos(x), x=Support(X, output='range'))

 


This graph is made of two "branches" where ths cosine function has a monotone variation:
          (1)  from -Pi to 0 cosine is monotonously increasing
          (2)  from 0 to Pi  cosine is monotonously decreasing


Watchout
The inverse cosine function is the arccosine function only on intervals where cosine is monotone.
So the reciprocal of  y = cos(x) over -p..p is not x = arccos(y).
Indeed:

solve(cos(x)=y, x, allsolutions)

2*Pi*_Z3-2*arccos(y)*_B3+arccos(y)

(4)

indets((4), name) minus {y, Pi}:
about~(%):

Originally _B3, renamed _B3~:

  is assumed to be: OrProp(0,1)

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

 


Choosing _Z3 = 0 gives
         

reciprocal := eval((4), _Z3=0)

-2*arccos(y)*_B3+arccos(y)

(5)


And the two branches are
         

phi[1] := unapply( eval(reciprocal, _B3=0), y);
phi[2] := unapply( eval(reciprocal, _B3=1), y);

proc (y) options operator, arrow; arccos(y) end proc

 

proc (y) options operator, arrow; -arccos(y) end proc

(6)


pdf  fX(x) of X
   

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

proc (x) options operator, arrow; piecewise(x < -Pi, 0, x < Pi, (1/2)/Pi, 0) end proc

(7)


The pdf  fY(y) of Y is given by this formula

result := Sum('f__X'(phi[i](y)) * abs(diff(phi[i](y), y)), i=1..2)

Sum(f__X(phi[i](y))*abs(diff(phi[i](y), y)), i = 1 .. 2)

(8)

result := value(result)

result := piecewise(arccos(y) < -Pi, 0, arccos(y) < Pi, 1/(2*Pi), 0)/sqrt(abs(y^2-1))+piecewise(-arccos(y) < -Pi, 0, -arccos(y) < Pi, 1/(2*Pi), 0)/sqrt(abs(y^2-1))

(9)

f__Y := unapply(simplify(result), y):
f__Y(y);

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

(10)

plot(f__Y(y), y=-1..1, title="Arcsine law")

 


ALTERNATIVE WAY: BUILD CDF(Y) AND DERIVE IT TO GET PDF(Y)


cdf  FX(x) of X
   

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

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

(11)


If all the fi(y) were increasing functions of y, the cdf  FY(y) of Y would write

result := Sum(chi('F__X'(phi[i](y))), i=1..2);

Sum(chi(F__X(phi[i](y))), i = 1 .. 2)

(12)


If some fi(y) are decreasing functions of y, the cdf  FY(y) of Y is defined by this formula

result := Sum(chi('F__X'(phi[i](y))), i=1..2), ``(chi('F__X'(phi[i](y))) = piecewise(diff(phi[i](y), y) > 0, 'F__X'(phi[i](y)), 1-'F__X'(phi[i](y))))

result := Sum(chi(`#msub(mi("F"),mi("X"))`(phi[i](y)))*`,`(chi(`#msub(mi("F"),mi("X"))`(phi[i](y))) = piecewise(0 < diff(phi[i](y), y), `#msub(mi("F"),mi("X"))`(phi[i](y)), 1-`#msub(mi("F"),mi("X"))`(phi[i](y)))), i = 1 .. 2)

(13)

plot([phi[1](y), phi[2](y)], y=-1..1, color=[red, blue], legend=[typeset(phi__1(y)), typeset(phi__2(y))])

 


Thus

result := ``(1 - F__X(phi[1](y))) + F__X(phi[2](y))

result*`:=`(1-piecewise(arccos(y) < -Pi, 0, arccos(y) < Pi, (1/2)*(arccos(y)+Pi)/Pi, 1))+piecewise(-arccos(y) < -Pi, 0, -arccos(y) < Pi, (1/2)*(Pi-arccos(y))/Pi, 1)

(14)

F__Y := unapply(simplify(expand~(result)), y):
F__Y(y);

piecewise(y < -1, 1, y <= 1, (Pi-arccos(y))/Pi, 1 < y, 1)

(15)


So PDF(Y) is:

diff(F__Y(y), y);

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

(16)


Verify that fY(y) obtained using the direct way (equation (10)) is almost everywhere equal to expression (16) 

"almost everywhere" means "out of the points of null probability measure (y=±1) where (16) is undefined.
   Note that the existence of these points come from the discontinuity of PDF(X) and that they could be removed from
   formulas (10) and (16) (because of their null probability measure).

simplify(f__Y(y) - diff(F__Y(y), y));

piecewise(y = -1, undefined, y = 1, undefined, 0)

(17)

plot(
  [f__Y(y), (16)]
  , y=-1..1
  , color=[red, blue]
  , thickness=[1, 5]
  , transparency=[0, 0.7]
  , legend=[typeset('f__Y'(y)), typeset('diff(F__Y(y), y)')]
)

 


NOTE THAT MAPLE IS CAPABLE TO COMPUTE PDF(Y) DIRECTLY

PDF(cos(X), y)

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

(18)


But I am pretty sure this is only by chance for, from what I have observed asking Maple to compute the PDF of a polynomial
function of X, Maple does not seem to identify "branches" nor to use formulas like (8) or (13).

In the cos(X) case, without no particular attention to the correctness of the formulas and by invoking wrongly a symmetry argument 
we would get:

WrongFormulaButCorrectResult := diff(simplify(2*F__X(phi[2](y))), y);

WrongFormulaButCorrectResult := diff(simplify(2*F__X(phi[1](y))), y);

WrongFormulaButCorrectResult := piecewise(y < -1, 0, y = -1, undefined, y < 1, 1/(sqrt(-y^2+1)*Pi), y = 1, undefined, 1 < y, 0)

 

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

(19)

V := RandomVariable(Uniform(-5*Pi/3, Pi/3)):

PDF(cos(V), y);    # Should be equal to (10) or (16)

plot(%, y=-1..1, title="Obviously wrong", titlefont=[Tahoma, bold, 14]);

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

 

 
 

 

Download cosine.mw

@JoyDivisionMan 

You're welcome!

@janhardo 

I wrote "at points (a, f(a)) and (-a, -f(a)) if and only if a=0... which is impossible because f(0) is undefined" because I didn't want to write explicitely a=1 and f(a)=2... or more exactly that a was strictly larger than zero.
Indeed, if a > 0 there is no circle tangent to the curve C at points at points (a, f(a)) and (-a, -f(a)).
The only possibility for such a circle to exist is a = 0, which cannot be accepted for the function f is not defined at 0.

@acer 

Thank you for your intervention.

If this can be of any interest here is a worksheet which shows when Maple succeeds in computing the PDF and CDF of a function 𝜙 of a random variable X :

  • Maple successes are illustrated by two examples:
    • the first one is trivial because 𝜙 is linear,
    • the second one is a priori more complex because 𝜙 is not monotonous... but the "branches" have simple analytic expressions and Maple can do the job.
  • Maple failure is illustrated by the  non monotonuous function 𝜙 : X ⟼X3/3 + X 
    (note that Maple succeeds if 𝜙 : X ⟼k*X3 because  𝜙 is now monotonuous).

I know how to compute the CDF of 𝜙(X) in an ad hoc way for arbitrary functions 𝜙, but I'm not capable to know how to do it programatically in a general case.
And I guess this is the lack of such an algorithm which pushes Maple to return FAIL for almost all non monotonuous functions 𝜙 (not to speak of no-failure situations where wrong results are returned... likely for the same lack of algorithm).

When_Maple_fails_in_computing_PDF_ot_CDF.mw

@JoyDivisionMan 

I don't know if you intend to follow @acer's suggestion of submitting an SCR (Software Change Request, see the screen shot below)


But if you do it, please inform the development team that Maple 2024 and beyond must not simply reproduce the result older versions provide (2015 or instance) because this result is wrong (for your "test" case) or might be wrong in case the 
function 𝜙 : X ⟼ Y is not a one-to-one map over the whole support of X.

There is indeed an error in the algorithm which computes PDF(Y) and this is it which must be corrected.

@JoyDivisionMan

Please look my pevious answer to understand why I limit my construction to a function 𝜙 : X ⟼Y which is a one-to-one C1 map over the support of X

restart

with(Statistics):


As a rule it's almost always simpler to compute the PDF of a function of a (several) random variables by computing
firstly their CDF and differentiating the result, instead than computing directly the PDF.

Here is an example: Computing PDF(sin(X)) where X ~ Uniform(-p/2, p/2).

Note that the support of X has been chosen in order that the sine function is a one-to-one map over this support. The
case where the X Y map is not one-to-one over the support of X, for instance if X ~ Uniform(0, 2p), is a little bit more delicate to treat (see my previous answer).

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

F__X := unapply(CDF(X, x), x):
F__X(x);

piecewise(x < -(1/2)*Pi, 0, x < (1/2)*Pi, ((1/2)*Pi+x)/Pi, 1)

(1)

phi := x -> sin(x);
psi := unapply(solve(phi(x)=y, x), y):

'psi'(y) = psi(y);

proc (x) options operator, arrow; sin(x) end proc

 

psi(y) = arcsin(y)

(2)

F__Y := unapply(F__X(psi(y)), y):
F__Y(y);

piecewise(arcsin(y) < -(1/2)*Pi, 0, arcsin(y) < (1/2)*Pi, (arcsin(y)+(1/2)*Pi)/Pi, 1)

(3)

f__Y := unapply(diff(F__Y(y), y), y);

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

(4)

plots:-display(
 Histogram(Sample(sin(X), 10^5), minbins=100, style=polygon, color="LightGray", transparency=0.5),
  plot(f__Y(y), y=-1..1, color=blue, thickness=3)
)

 
 

 

Download PDF_of_a_function_of_a_random_variable.mw

@Andiguys 

Your expressions are that complex that I am pretty sure they come from some calculus you do not give the details here, and whose results you copied-pasted into the worksheet you provide.

The reason why you get those complex expressions is likely a consequence of you doing something like
x1 := some expression containing quantities q1, .., qn
x2 := some expression containing some quantities among q1, ... ,qn plus extra quanrities r1, ..., rn
x3 := some expression containing some quantities among q1, ..., qn and r1, .., rn, plus extra quanrities s1, ..., sn
... and so on up to a quite simple expression of the form z = f(x1, x2, ..., xm)

Doing this will automatically result in an extremely complec form for z.

A better way to proceed if you want to get a concise expression for z is to write instead
z := f(x1, x2, ..., xm)
substitutions := {x1 = ?..., x2= ..., x3= ..., ...}

Now z appears in its simplest form and you don't have t to worry about how to "shorten" the expression of z.
If you run the command eval(z, substitutions), you will get an extremely complex expression that you will probably struggle to “shorten” significantly.

Here are two simple illustrations WhatYouMustNotDo.mw WhatYouMustNotDo_2.mw of what I said.

Before CAS existed, people thought twice before embarking on laborious calculations and tried to simplify the initial expressions as much as possible before undertaking these calculations, sometimes even assigning other names to sub-expressions during the calculation.
Today's CAS are so powerful that they eliminate the need to think, leading many users to want to simplify results that they have unnecessarily complicated out of laziness.
Try to avoid taking this way too.

@C_R 

I need to think about it before giving a solid answer.

Sorry.

@C_R 

Yes It is indeed a pity that Maple does not offer a function to draw (in a correct way) the cumulative distribution of discrete random variable or of a sample.

The CumulativeSumChart (see CUSUM) is essentially used in "quality control" and has nothing to do with the cumulative function (the name CumulativeSumChart can be misleading). Note that the latter is a non decreasing function while the former can have arbitrary variations. 

Cumulative Sum Charts should be used to plot stochastic process (the corresponding help page says it plots  "the cumulative sum chart for the specified data" but unfortunately gives no reference to indicate what this cumulative sum chart, which can misled those who use the CumulativeSumChart function.

For fun, look Here.mw .

@janhardo 

I vote up

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