sand15

Explanations...

Answered: sand15 1067

October 03 2025

1 0

@JoyDivisionMan

Using "blindly" Maple (2015) provides a wrong result (see attached file) for, I guess, maple does not handle correctly the absolute value.

So here is he trick.
Let f_Z(z) the pdf of Z = |Y₁- Y₂| and f_U(u) the pdf of U = Y₁- Y₂.
Of course the support of Z is (0, 2) while U's is (-2, 2). But Y₁ and Y₂ both having the same symmetic distribution it is easy to see that f_U(u) is symmetric too.
So consider only the restriction of f_U(u) the its half support (0, 2) and twice this latter: then 2f_U(u) = f_Z(u).

With that trixk Maple provides the corret result.

Note that f_U(u) can be computed using a convolution product.
Indeed U = Y₁ - Y₂ = Y₁ + (-Y₂) = Y₁ + Y₂ because f_Y(y) is symmetrc wrt y=0.
Then f_U(u) = (f_Y * f_Y)(y) where '*' denotes the convolution produce.

Please, if that's okay with you, don't forget (as it was the case with your last question) to let me know.
(Doing work for someone without knowing whether they appreciated it or not is never pleasant.)

>

restart

>	kernelopts(version)

(1)

>	with(Statistics):

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

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);

proc (x) options operator, arrow; piecewise(x <= -1, 0, x <= 1, (1/2)*(2*x*(1-x^2)^(1/2)+2*arcsin(x)+Pi)/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:

>	# Note the ∞ term in the case 0 < z < 1 !!! # This is obviously a bug. PDF(Z, z);

(4)

>	# 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);

(5)

>	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)

>	# 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

(6)

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

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

(7)

>	:-Median('Z') = Median(Z); :-Median('Z') = Median(Z, numeric);

(8)

>	:-Median('Z') = fsolve('F__U'(u)=0.5, u=0..2); med := rhs(%):

(9)

Spline approximation of F__U

>	S__U := unapply(CurveFitting:-Spline([seq([u, F__U(u)], u=0..2, 0.25)], u), u):

Procedure to compute quantiles of Z

>

quantile := proc(p)
  local guess:
  if F__U(1) > p then
  # guess := fsolve(S__U(u)=p, u=0..1)
    guess := fsolve('F__U'(u)=p, u=0..1)
  else
    guess := fsolve(S__U(u)=p, u=1..2)
  end if:
  return [guess, ``(F__U(guess))]:
end proc:

:-Quantile('Z', 0.05) = quantile(0.05); # Exact, based on F__U
:-Quantile('Z', 0.80) = quantile(0.80); # Exact, based on F__U
:-Quantile('Z', 0.85) = quantile(0.85); # Approximation based on S__U
:-Quantile('Z', 0.95) = quantile(0.95); # Approximation bBased on S__U

(10)

>	L := limit(f__U(u), u=0, right); if false then plots:-display( plot(f__U(u), u=0..2, color=red), plottools:-transform((x, y) -> [2-x, L-y])(plot(f__U(u), u=0..2, color=blue)) ); end if;

(32/3)/Pi^2

(11)

>

Download PDF_Z_sand15.mw

KKT doesn't seem necessary here...

Answered: sand15 1067

September 30 2025

1 11

𝛑_m is a priori quadratic wrt i₁ so 𝛑_m has only one extremum swhose status (minimum vs maximum) depends on the sign of the coefficient of i₁².

Note: I say "𝛑_m is a priori quadratic wrt i₁" for 𝛑_m writes P . i₁² + Q . i₁ + R where P, Q, R are expressions which depend on several parameters. It could therefore happen that, depending on the values of these parameters, P, or even P and Q, might be equal to 0.
The case P = 0, Q <> 0 is even simpler because 𝛑_m necessary reaches its the minimum at one of the two boundaries defined by A and B (see below the meaning of these two quantities)

As you did not provide more informations I considered arbitrarily that P was different from 0.

Let x* the location of this extremum and D the domain defined by the two primal feasibility conditions

f₁ : A <= i₁
f₂ : i₁ <= B (B supposedly larger than A of course)

Assuming for instance that P <> 0 and that 𝛑_m reaches its minimum at some point i₁= x*, three cases are to be examined:

x* < A:
The constrained minimum is x* = A
x* > B:
The constrained minimum is x* = B
A < x* > B:
The constrained minimum is x*

Assuming now that P = 0 and Q <> 0, we find the location of the minimizer i₁= x* of 𝛑_m corresponds to one of these two cases:

𝛑_m(A) < 𝛑_m(B)
The constrained minimum is x* = A
𝛑_m(A) > 𝛑_m(B)
The constrained minimum is x* = B

I don't really think it's necessary to write KKT conditions and solve grad(Lagrangian(i₁, 𝛍₁, 𝛍₂)) = 0 wrt i₁, 𝛍₁, 𝛍₂ to get this result.

Basically your problem can be parameterized by the four quantities P, Q, A, B.

Why not simply use solve/parametric to get all the possible solutions?

>

restart

>	SC_CS_sols := solve( { diff(Pi1^2 + Qi1 + R, i1) - mu[1] + mu[2] = 0, mu[1](A-i1) = 0, mu[2](i1-B) = 0 } , {i1, mu[1], mu[2]} , 'parametric'='full', 'parameters'={P, Q, A, B} )

-

(1)

>	pi__m := (w-i1)(1/2+(i1-i2)/(2tau))(1-(Pn-Pr)/(1-delta))-C0-(1/2)Cr(1/2+(i1-i2)/(2tau))^2(1-(Pn-Pr)/(1-delta))^2+Cerho0*(1-(Pn-Pr)/(1-delta)): collect(pi__m, i1): PQR := [R, Q, P] =~ [coeffs(%, i1)]: print~(PQR):

R = w*(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))-C0-(1/2)*Cr*(1/2-(1/2)*i2/tau)^2*(1-(Pn-Pr)/(1-delta))^2+Ce*rho0*(1-(Pn-Pr)/(1-delta))

Q = ((1/2)*w/tau-1/2+(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))-(1/2)*Cr*(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))^2/tau

P = -(1/2)*(1-(Pn-Pr)/(1-delta))/tau-(1/8)*Cr*(1-(Pn-Pr)/(1-delta))^2/tau^2

(2)

>	AB := [A = (2rho0-1)tau+i2, B = tau+i2]

(3)

>	ind := indets(rhs~(PQR)) union indets(rhs~(AB));

(4)

>

randomize();

data := { seq(i = rand(0. .. 1.)(), i in ind) };
PQRAB := { eval(PQR, data)[], eval(AB, data)[] };

all_solutions := allvalues~(eval(SC_CS_sols, PQRAB)):
all_pi        := map(s -> eval(eval(pi__m, data), s), all_solutions):
here          := ListTools:-Search(min(all_pi), all_pi):
SOLUTION      := [all_pi[here], all_solutions[here]];

plot(eval(pi__m, data), i1=eval(eval(A..B, AB), data), gridlines=true)

>

Download proposal_sand15.mw

elementwise comparisons...

Answered: sand15 1067

September 28 2025

0 2

See the Maple help page for the meaning of the elementwise operator '~'

a := Matrix([[1, 2], [1, 2]]):
is~(a[1] =~ a[2]);
    
b := [seq(a[i], i=1..2)]:
is~(b[1] =~ b[2]);

Suggestion...

Answered: sand15 1067

September 28 2025

0 1

Look HERE to my reply to a recent similar question.

Simply...

Answered: sand15 1067

September 27 2025

3 1

The lazzy attitude: simply ask Maple the result

>

restart

>	with(Statistics):

>	f[1] := t -> piecewise(t <= 0, 0, 0 < t and t < 2, 1/(Pisqrt(1-(1-t)^2)), t >= 2, 0); f[2] := t -> piecewise(t <= 0, 0, 0 < t and t < 1, 2t, t >= 1, 0);

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

(1)

>	Dist1 := Distribution(PDF = (t -> f[1](t))): Dist2 := Distribution(PDF = (t -> f[2](t))):

>	X := RandomVariable(Dist1): Y := RandomVariable(Dist2): Z := X*Y:

>	`PDF(Z)` := PDF(Z, t); `CDF(Z)` := CDF(Z, t);

`PDF(Z)` := piecewise(t < 0, 0, t <= 2, -(2/3)*(t^2-t-2)/(sqrt(t)*sqrt(2-t)*Pi), 2 < t, 0)

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

(2)

>	plot(`PDF(Z)`, t=0..2)

>

Download Inquiry_sand15.mw

A classical result (under some restricions)
Let X and Y two independent random variables whose both supports belong to (0, +∞) and with respective pdf f_X and f_Y then the pdf f_Z of Z = X Y is given by

For more details you can see for instance Glen at al.or go to Google Scholar and search for "product of two independent random variables".

If you want to understand how to get the result while not knowing about this formula...
As a rule
To obtain the probability density function of the function h of one (or more) continuous random variable(s) it is always simpler to compute its cumulative function first and next to differentiate it (assuming this cumulative function is derivable).
Thus is all the more simple if h has an unique inverse over the support of this (these) random variable(s) (if not some precautions have to be taken).
This your case: Z = h( X , Y ) = X Y.
Worksheet PROCEDURE.mw describes the step-by-step procedure to get f_Z :

>

restart

Let X and Y two independent random variables with respective probability functions

>	f__X := x -> piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pisqrt(1-(1-x)^2)), x >= 2, 0); f__Y := y -> piecewise(y <= 0, 0, 0 < t and y < 1, 2y, y >= 1, 0);

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

(1)

Their cumulative functions write

>	F__X := unapply(int(f__X(t), t=-infinity..x), x); F__Y := unapply(int(f__Y(t), t=-infinity..y), y);

(2)

Let C = X Y.

>

Prob(Z < z) = Prob(XY < z);
# Independency of X and Y gives

Prob(Z < z) = Prob(X < x and Y < z/x);

# Prob(Z < z | X < x) reads: Prob(Z < z) given that X < x

Prob(cat(Z < z, ` | `, X = x)) = Prob(Y < z/x);

# By definition of F_Y

Prob(cat(Z < z, ` | `, X = x)) = F__Y(z/x);

# To get Prob(Z < z) integrate Prob(Z < z | X < x) for all x values (accounting for the PDF of X of course)
Prob(Z < z) = Int(F__Y(z/x)*f__X(x), x=0..+infinity);

# By definition of F_Z

F__Z(z) = rhs(%);

# Derivation both sides wrt z gives the formal expression of f_Z

f__Z(z) = Diff(rhs(%), z);

# Exchange Diff and Int (assuming this operation is licit)

f__Z(z) = Int(Diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# Evaluate the derivative of the integrand

f__Z(z) = Int(diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# Evaluate the integral of the previous result

f__Z(z) = int(diff(F__Y(z/x)*f__X(x), z), x=0..+infinity);

# (You see that the Mellin transform is absolutely useless here.)

Prob(`Z < z | X = x`) = piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)

Prob(Z < z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

`#msub(mi("F"),mi("Z"))`(z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

`#msub(mi("f"),mi("Z"))`(z) = Diff(Int(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity), z)

`#msub(mi("f"),mi("Z"))`(z) = Int(Diff(piecewise(z/x <= 0, 0, z/x <= 1, z^2/x^2, 1 < z/x, 1)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), z), x = 0 .. infinity)

`#msub(mi("f"),mi("Z"))`(z) = Int(piecewise(z/x <= 0, 0, z/x <= 1, 2*z/x^2, 1 < z/x, 0)*piecewise(x <= 0, 0, 0 < x and x < 2, 1/(Pi*sqrt(1-(1-x)^2)), 2 <= x, 0), x = 0 .. infinity)

f__Z(z) = piecewise(z < 0, 0, z < 2, -(2/3)*(z^2-z-2)/((z*(2-z))^(1/2)*Pi), 2 <= z, 0)

(3)

>	plot(rhs(%), z=0..2)

>

Worksheet PROCEDURE_formal.mw contains the demonstration of the classical result when X and Y are independent random variables with supports in (0, +∞).

As you already asked a question about the convolution product, here is a useful trick
Let A=ln( X ) and B=ln( Y ), the pdf of C=ln( Z ) is simply the convolution product of the pdfs of A and B.
Once you know the pdf of C it remains simply to compute the one of exp( C ) .
See worksheet TRICK.mw for more details.

No solution...

Answered: sand15 1067

September 26 2025

3 3

Unless Maple is wrong or if rewritting the system in a more concise form removes some simplifications or factorizations, it doesn't seem a solution exists:

>

K2 := -(Pn-Pr)/(1-delta)+(-beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+Pr)/delta+(Pr-w-Cr)*beta*(-2+2*delta)*tau*upsilon/(((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))*delta) = 0:

>

K3 := 1-(Pn-Pr)/(1-delta)-(Pn-Cn)/(1-delta)+(Pr-w-Cr)*(1/(1-delta)-(-beta*(Cr*i2-Cr*tau)*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon*Cr/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))^2)/delta)+Ce*rho0/(1-delta) = 0:

>

Error, (in LinearAlgebra:-Determinant) time expired

A good strategy almost always consits in rewritting the equations in a more concise form

>

rewrite := proc(e, v)
  local f1, f2, C, M, f3:
  # assuming that the roots of the numerator do not cancel out the denominator:
  f1 := numer(lhs(e));

  C := [coeffs(f1, [Pn, Pr, w], 'M')];
  M := [M];
  f2 := [seq([seq(degree(t, i), i in [Pn, Pr, w])], t in M)];
  f3 := add(v[f2[i]]*M[i], i=1..nops(M));
  return f3=0, [seq(v[f2[i]] = C[i], i=1..nops(M)) ];
end proc:

newK2, rw2 := rewrite(K2, U):
newK3, rw4 := rewrite(K3, V):
newK4, rw4 := rewrite(K4, W):

Example: K2 rewritten and the corresponding values of the U coefficients:

>	newK2; print(): print~(rw2):

U[[0, 0, 0]] = -Cr*beta*delta^2*i2*upsilon-Cr*beta*delta^2*tau*upsilon-2*beta*delta^2*i2*tau*upsilon+2*beta*delta^2*tau^2*upsilon+2*Cr*beta*delta*i2*upsilon+2*Cr*beta*delta*tau*upsilon+4*beta*delta*i2*tau*upsilon-4*beta*delta*tau^2*upsilon-Cr*beta*i2*upsilon-Cr*beta*tau*upsilon-2*beta*i2*tau*upsilon+2*beta*tau^2*upsilon

(1)

Here is the rewritten system

>	sys := {newK2, newK3, newK4}: for e in sys do e; print(); end do;

Pn^3*V[[3, 0, 0]]+Pn^2*Pr*V[[2, 1, 0]]+Pn^2*w*V[[2, 0, 1]]+Pn*Pr^2*V[[1, 2, 0]]+Pn*Pr*w*V[[1, 1, 1]]+Pr^3*V[[0, 3, 0]]+Pr^2*w*V[[0, 2, 1]]+Pn^2*V[[2, 0, 0]]+Pn*Pr*V[[1, 1, 0]]+Pn*w*V[[1, 0, 1]]+Pr^2*V[[0, 2, 0]]+Pr*w*V[[0, 1, 1]]+w^2*V[[0, 0, 2]]+Pn*V[[1, 0, 0]]+Pr*V[[0, 1, 0]]+w*V[[0, 0, 1]]+V[[0, 0, 0]] = 0

Pn^3*W[[3, 0, 0]]+Pn^2*Pr*W[[2, 1, 0]]+Pn^2*w*W[[2, 0, 1]]+Pn*Pr^2*W[[1, 2, 0]]+Pn*Pr*w*W[[1, 1, 1]]+Pr^3*W[[0, 3, 0]]+Pr^2*w*W[[0, 2, 1]]+Pn^2*W[[2, 0, 0]]+Pn*Pr*W[[1, 1, 0]]+Pn*w*W[[1, 0, 1]]+Pr^2*W[[0, 2, 0]]+Pr*w*W[[0, 1, 1]]+w^2*W[[0, 0, 2]]+Pn*W[[1, 0, 0]]+Pr*W[[0, 1, 0]]+w*W[[0, 0, 1]]+W[[0, 0, 0]] = 0

(2)

newK2 being linear wrt w just solve it wrt w:

>	infolevel[solve] := 4: wsol_1 := solve(newK2, w);

Main: Entering solver with 1 equation in 1 variable

Transformer: solving the uncoupled linear subsystem in w

Main: solving successful - now forming solutions
Main: Exiting solver returning 1 solution

-(Pn^2*U[[2, 0, 0]]+Pn*Pr*U[[1, 1, 0]]+Pr^2*U[[0, 2, 0]]+Pn*U[[1, 0, 0]]+Pr*U[[0, 1, 0]]+U[[0, 0, 0]])/U[[0, 0, 1]]

(3)

Form a 2-by-2 subsystem while plugging the solution above into {newK3 and newK4}

>	sys34 := eval({newK3, newK4}, w=wsol_1): newnewK3, newrw3 := rewrite(sys34[1], X): newnewK4, newrw4 := rewrite(sys34[2], Y):

>	newsys := {newnewK3, newnewK4}: for e in newsys do e; print(); end do;

Pn^4*X[[4, 0, 0]]+Pn^3*Pr*X[[3, 1, 0]]+Pn^2*Pr^2*X[[2, 2, 0]]+Pn*Pr^3*X[[1, 3, 0]]+Pr^4*X[[0, 4, 0]]+Pn^3*X[[3, 0, 0]]+Pn^2*Pr*X[[2, 1, 0]]+Pn*Pr^2*X[[1, 2, 0]]+Pr^3*X[[0, 3, 0]]+Pn^2*X[[2, 0, 0]]+Pn*Pr*X[[1, 1, 0]]+Pr^2*X[[0, 2, 0]]+Pn*X[[1, 0, 0]]+Pr*X[[0, 1, 0]]+X[[0, 0, 0]] = 0

Pn^4*Y[[4, 0, 0]]+Pn^3*Pr*Y[[3, 1, 0]]+Pn^2*Pr^2*Y[[2, 2, 0]]+Pn*Pr^3*Y[[1, 3, 0]]+Pr^4*Y[[0, 4, 0]]+Pn^3*Y[[3, 0, 0]]+Pn^2*Pr*Y[[2, 1, 0]]+Pn*Pr^2*Y[[1, 2, 0]]+Pr^3*Y[[0, 3, 0]]+Pn^2*Y[[2, 0, 0]]+Pn*Pr*Y[[1, 1, 0]]+Pr^2*Y[[0, 2, 0]]+Pn*Y[[1, 0, 0]]+Pr*Y[[0, 1, 0]]+Y[[0, 0, 0]] = 0

(4)

Note that the lhs of each equations are complete polynomials of total degree 4 in the two indeterminates Pn and Pr

>	infolevel[solve] := 4: sol := timelimit(20, solve(newsys, {Pn, Pr}));

Main: Entering solver with 2 equations in 2 variables
Main: attempting to solve as a linear system
Main: attempting to solve as a polynomial system
Main: polynomial system split into 1 parts under preprocessing
Main: trying resultant methods
PseudoResultant: 4093009 [2 7000140160 Pn] 2 2 649 0 3 0

PseudoResultant: factoring all equations length = 649
PseudoResultant: factoring done length = 649
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512

LinearEq: subresultant determinant of length 92449 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512
LinearEq: subresultant determinant of length 92449 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 322 and 322

Error, (in convert/sqrfree/sqrfree) time expired

>	timelimit(20, SolveTools:-PolynomialSystem(lhs~(newsys), {Pn, Pr}));

Main: polynomial system split into 1 parts under preprocessing

Main: trying resultant methods
PseudoResultant: 4445613 [2 7000142560 Pn] 2 2 809 0 3 0
PseudoResultant: factoring all equations length = 809

PseudoResultant: factoring done length = 809
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850
LinearEq: subresultant determinant of length 127039 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850

LinearEq: subresultant determinant of length 127039 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 402 and 402
SolutionsLost: setting solutions lost flag
PseudoResultant: 0 solutions found, now doing backsubstitution

SolutionsLost: setting solutions lost flag

Seems to be a dead line

Download Q_solve_sand15.mw

Here it is...

Answered: sand15 1067

September 23 2025

1 2

>

restart

>	aux := B12 = -6*(p1 + p2)/(p1 - p2)^2;

(1)

>	F2 := (theta1theta2p1^2 + (-2p2theta1theta2 - 6)p1 + theta1theta2p2^2 - 6*p2)/(p1 - p2)^2;

(theta1*theta2*p1^2+(-2*p2*theta1*theta2-6)*p1+theta1*theta2*p2^2-6*p2)/(p1-p2)^2

(2)

>	collect(expand(F2), theta1):; map(factor, %); F2a := op(1, %) + simplify(op(2, %)+op(3, %));

(3)

>	F2b := eval(F2a, (rhs=lhs)(aux))

(4)

Download rewrite.mw

Linearize to avoid non convergence issue...

Answered: sand15 1067

September 23 2025

4 2

Main remark:
Your model y = P(x)/Q(x) where P and Q are two polynomials of x can be linearized under a convenient
reparameterization: basically you discard the regressor x and introduce new ones which make the model linear).
Once done you can use  Statistics:-Fit instead of Statistics:-NonlinearFit (in fact Optimization:-LSSolve) and thus discard any convergence issue.
More of this you avoid getting a solution which depends on the initial values
(you already noticed that the closer to the values used in the generative model, the better the results [which supposes this
  generative model is noiseless])
and you get access to all the output statisctics than  Statistics:-Fit offers)

Awaiting your data Ycg and Sx, here is a notional example (same model than yours but data produced by a noisy process) which explains how to linearize the model and go back to its initial nonlinear form.

>

restart

>	with(Statistics):

>	Model__ini := y = (p1x^2 + p2x + p3)/(q1*x + x^2 + q2)

y = (p1*x^2+p2*x+p3)/(q1*x+x^2+q2)

(1)

It is more practical do rewrite Model₀ this way

>	Model__0 := y = (a[2]x^2 + a[1]x + a[0])/(b[2]x^2 + b[1]x + 1)

y = (x^2*a[2]+x*a[1]+a[0])/(x^2*b[2]+x*b[1]+1)

(2)

Assuming the denominator is never null in some range of x, this model can be linearized.
Indeed:

>	Model__1 := expand(lhs(Model__0)*denom(rhs(Model__0))) = numer(rhs(Model__0))

(3)

>	VarChange := {x^2yb[2] = Z__21b[2], xyb[1] = Z__11b[1], x^2a[2] = Z__20a[2], xa[1] = Z__10a[1]}

${a[1]*x = Z__10*a[1], a[2]*x^2 = Z__20*a[2], x*y*b[1] = Z__11*b[1], x^2*y*b[2] = Z__21*b[2]}$

(4)

>	Model__2 := subs(VarChange, Model__1)

(5)

>	Model__lin := isolate(Model__2, y)

(6)

Example.
I use as a generative process for Sx and Ycg (given the order you use in your question Ycg is a vector of
x values and Sx a vector of y values).
For instance

>	Seed := 175860914374622; randomize(Seed);

(7)

>

GenerativeProcess := proc(N, sigma::positive)
  local R, params, values, model, f, X, Y:

  R      := () -> rand(-10..10)() / rand(1..20)():

  params := convert(indets(Model__0) minus {x, y}, list):
  values := [seq(R(), n=1..numelems(params))]:
  model := eval(Model__0, params =~ values);

  print(model + Noemal(0, sigma)):

  f := unapply(rhs(model), x):
  X := Sample(Uniform(-1, 1), N)^+:
  Y := f~(X) + Sample(Normal(0, sigma), N)^+:
  return X, Y, params =~ values
end proc:

>	N := 40: sigma := 0.1: Ycg, Sx, params := GenerativeProcess(N, sigma);

y+Noemal(0, .1) = (-(10/9)*x^2+(10/17)*x-3/10)/(-(2/3)*x^2+(7/11)*x+1)+Noemal(0, .1)

Vector[column](%id = 18446744078182577086), Vector[column](%id = 18446744078182577686), [a[0] = -3/10, a[1] = 10/17, a[2] = -10/9, b[1] = 7/11, b[2] = -2/3]

(8)

Construct new regressors

>	Model__lin; z__10 := Ycg: z__11 := Ycg ~Sx: z__20 := Ycg^~2: z__21 := Ycg^~2 ~Sx:

(9)

Fit linear model

>	res := Fit( rhs(Model__lin) , `<\|>`(z__10, z__11, z__20, z__21) , Sx , [Z__10, Z__11, Z__20, Z__21] , output=[leastsquaresfunction, parametervalues] )

[-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*Z__10-HFloat(1.1869703402160423)*Z__20-HFloat(0.6608454314903524)*Z__11+HFloat(0.6374123072191089)*Z__21, [a[0] = HFloat(-0.2777298343064421), a[1] = HFloat(0.5262609535304735), a[2] = HFloat(-1.1869703402160423), b[1] = HFloat(0.6608454314903524), b[2] = HFloat(-0.6374123072191089)]]

(10)

Go back to original model with regressor x

>	Z_fit := res[1]

-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*Z__10-HFloat(1.1869703402160423)*Z__20-HFloat(0.6608454314903524)*Z__11+HFloat(0.6374123072191089)*Z__21

(11)

>	aux := indets(Model__lin) minus indets(Model__0): Z2xy := solve(VarChange, aux); X_fit := isolate(y = eval(Z_fit, Z2xy), y);

${Z__10 = x, Z__11 = x*y, Z__20 = x^2, Z__21 = x^2*y}$

y = (-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*x-HFloat(1.1869703402160423)*x^2)/(1+HFloat(0.6608454314903524)*x-HFloat(0.6374123072191089)*x^2)

(12)

Compare observations to model predictions

>	plots:-display( ScatterPlot(Ycg, Sx) , plot(rhs(X_fit), x=(min..max)(Ycg)) )

Compare the values of the parameters for the generative process and their estimations

>	`<\|>`( `<,>`("True values", "Fitted values") , convert([evalf(params), res[2]], Matrix) )

Matrix(2, 6, {(1, 1) = "True values", (1, 2) = a[0] = -.3000000000, (1, 3) = a[1] = .5882352941, (1, 4) = a[2] = -1.111111111, (1, 5) = b[1] = .6363636364, (1, 6) = b[2] = -.6666666667, (2, 1) = "Fitted values", (2, 2) = a[0] = -.277729834306442, (2, 3) = a[1] = .526260953530474, (2, 4) = a[2] = -1.18697034021604, (2, 5) = b[1] = .660845431490352, (2, 6) = b[2] = -.637412307219109})

(13)

>

Download Illustrative_example.mw

Add-on:
You could find instructive this comparison to NonlinearFit: Illustrative_example_bad_NonlinearFit.mw

A solution...

Answered: sand15 1067

September 21 2025

1 0

Based on Units but focuses only (as an example) on dimension LENGTH ans SI units (so "are" and "hectare" are not considered here).

>

restart

>	kernelopts(version)

(1)

>	with(Units):

>	randomize(175847324662552);

(2)

For lengths, areas and volumes only

>	Ulength := [seq(cat(u, m), u in [m, c, d, ``, da, h, k])];

(3)

>

Quizz := proc(N)
  local R, `Q&A`, n, r, u, U, e, x, X:

  R := proc()
    local n, e:
    n := rand(0..100)():
    e := rand(-6..6)():
    if e <= 0 then
      return evalf[2](n*10^e)
    else
      return n*10^e
    end if:
  end proc:

  `Q&A` := NULL:
  for n from 1 to N do
    r     := R():
    u, U := op(combinat:-randperm(Ulength)[1..2]):
    e     := rand(1..3)():
    x     := r*u^e:
    X     := convert(r*Unit(u^e), 'units', U^e):

    printf("\n%a = .......... %a\n", x, U^e):
    `Q&A` := `Q&A`, Printf("\n%a = %a %a\n", x, evalf[4](op(1, X)), U^e):
  end do:
  print():
  return [`Q&A`]
end proc:

>	answers := Quizz(5):

.14e-2*dam^3 = .......... mm^3

2.*dm^3 = .......... dam^3

.73*m^3 = .......... cm^3

3.*km^2 = .......... dam^2

89.*km^3 = .......... hm^3

(4)

>	# Correction map(evalf, eval(answers, Printf=printf))

.14e-2*dam^3 = .1400e10 mm^3

2.*dm^3 = .2000e-5 dam^3

.73*m^3 = .7300e6 cm^3

3.*km^2 = .3000e5 dam^2

89.*km^3 = .8900e5 hm^3

(5)

>

Download Units_Quizz_2.mw

An evolution which also enables quizz about velocities Units_Quizz_3.mw

discont...

Answered: sand15 1067

September 19 2025

1 1

discont(V, beta)
                            {-1, 0}

The disontinuity you see at x= 8 .10¹⁷ is an artifact due to number of Digits not large enough:

>

s := C*((1-theta/lambda)/(1+beta))^(1/beta):

>

V := (lambda*(n-m)*s*(1-(s/C)^beta)-k*s)/alpha;

16.50000000*(0.909090909e-1/(1+beta))^(1/beta)*(1-((0.909090909e-1/(1+beta))^(1/beta))^beta)-10.0*(0.909090909e-1/(1+beta))^(1/beta)

(1)

>	discont(V, beta)

(2)

>

plot(V, beta = 1 .. 10^18, color = [red]);

>	Digits := 20: plot(V, beta = 1 .. 10^18)

>	# Or better: plot(eval(V, beta=10^u), u=0..18)

Download MaplePrimes_Sep_19_sand15.mw

No solutions in the search domain you us...

Answered: sand15 1067

September 15 2025

1 1

>

restart;

>	with(plots):

>

q := 0.9:
M := 1.01:

points := []:

for tau from 0.01 by 0.01 to 0.09 do

    eq1 := M^2*(1 - sqrt(1 - 2*x/M^2))
           + f*tau*(1 - exp(x/tau))
           + (2*(1-f)/(3*q-1))*(1 - (1+(q-1)*x)^((3*q-1)/(2*q-2))) = 0:

    eq2 := 1/sqrt(1 - 2*x/M^2)
           - f*exp(x/tau)
           - (1-f)*(1+(q-1)*x)^((q+1)/(2*q-2)) = 0:


    sol := fsolve({eq1, eq2}, {f, x}, {f=0.01..0.9, x=0.001..2});

    if sol::function then
        printf("tau=%.2f -> no solution found\n", tau);
    else
        fval := eval(f, sol);
        xval := eval(x, sol);
        points := [op(points), [fval, xval]];
        printf("tau=%.2f -> f=%.6f, x=%.6f\n", tau, fval, xval);
    end if:

od:

tau=0.01 -> no solution found

tau=0.02 -> no solution found

tau=0.03 -> no solution found

tau=0.04 -> no solution found

tau=0.05 -> no solution found

tau=0.06 -> no solution found

tau=0.07 -> no solution found

tau=0.08 -> no solution found

tau=0.09 -> no solution found

>

# Change search ranges for fsolve or simply ignore them

points := []:

for tau from 0.01 by 0.01 to 0.09 do

    eq1 := M^2*(1 - sqrt(1 - 2*x/M^2))
           + f*tau*(1 - exp(x/tau))
           + (2*(1-f)/(3*q-1))*(1 - (1+(q-1)*x)^((3*q-1)/(2*q-2))) = 0:

    eq2 := 1/sqrt(1 - 2*x/M^2)
           - f*exp(x/tau)
           - (1-f)*(1+(q-1)*x)^((q+1)/(2*q-2)) = 0:


    sol := fsolve({eq1, eq2}, {f, x});

    if sol::function then
        printf("tau=%.2f -> no solution found\n", tau);
    else
        fval := eval(f, sol);
        xval := eval(x, sol);
        points := [op(points), [fval, xval]];
        printf("tau=%.2f -> f=%.6f, x=%.6f\n", tau, fval, xval);
    end if:

od:

tau=0.01 -> f=0.851356, x=0.000000

tau=0.02 -> f=0.856683, x=0.000000

tau=0.03 -> f=0.001070, x=0.032805

tau=0.04 -> f=0.001678, x=0.059488

tau=0.05 -> f=0.002478, x=0.088687

tau=0.06 -> f=0.003503, x=0.121100

tau=0.07 -> f=0.004804, x=0.158068

tau=0.08 -> f=0.006453, x=0.202525

tau=0.09 -> f=0.008590, x=0.264387

>

points

[[.8513555201, 0.6964830318e-12], [.8566830887, 0.3912817906e-13], [0.1069988964e-2, 0.3280452951e-1], [0.1678278783e-2, 0.5948758155e-1], [0.2477635868e-2, 0.8868702371e-1], [0.3503054692e-2, .1211004031], [0.4803501518e-2, .1580677255], [0.6453490291e-2, .2025245703], [0.8590449799e-2, .2643866253]]

(1)

>	# You can select a posterior the points which suit you: MyPoints := select(t -> is(t[1] > 0 and t[2] < 0.2), points); print(): MyPoints := select(t -> is(t[1] > 0 and t[2] < 0.2 and t[2] > 0.01), points);

[[.8513555201, 0.6964830318e-12], [.8566830887, 0.3912817906e-13], [0.1069988964e-2, 0.3280452951e-1], [0.1678278783e-2, 0.5948758155e-1], [0.2477635868e-2, 0.8868702371e-1], [0.3503054692e-2, .1211004031], [0.4803501518e-2, .1580677255]]

(2)

>	pointplot(points, symbol=solidcircle, color=blue, labels=["f","x"]):

Download solve_sand15.mw

No type declaration: Type depends on the...

Answered: sand15 1067

September 14 2025

2 3

The type depends on how you defne the distribution.
If some attributes are those of a discrete random variable than Maple understands that this random variable is discrete.

>

restart

>	with(Statistics):

PDF should be resereved to continuous random variable

>	dist := Distribution( PDF = unapply( piecewise(x > 0, exp(-6)*6^x/x!, 0), x) ): X := RandomVariable(dist): exports(attributes(X)[3]); PDF(X, x)

piecewise(0 < x, exp(-6)*6^x/factorial(x), 0)

(1)

For a discrete random variable use ProbabilityFunction instead

>	dist := Distribution( ProbabilityFunction = unapply( piecewise(x > 0, exp(-6)*6^x/x!, 0), x) ): X := RandomVariable(dist): exports(attributes(X)[3]); (attributes(X)[3]):-Type;

(2)

Look how Maple interpretes the probability function in terms of PDF

>	:-ProbabilityFunction('X') = ProbabilityFunction(X, x); :-PDF('X') = PDF(X, x)

ProbabilityFunction(X) = piecewise(0 < x, exp(-6)*6^x/factorial(x), 0)

PDF(X) = sum(exp(-6)*6^k*Dirac(x-k)/factorial(k), k = 1 .. infinity)

(3)

Download discrete.mw

BY THE WAY: A correct definition of the ProbabilityFunction is

ProbabilityFunction = unapply( piecewise(x >= 0, exp(-6)*6^x/x!, 0), x)

Indeed

# With your definition 

sum(PDF(X, x), x=-infinity..+infinity) = 1-e^-6# With my correction above

sum(PDF(X, x), x=-infinity..+infinity) = 1

Once correctly defined, the probability function is that of a Poisson random variable with parameter equal to 6.
I don't know if you were aware of this and if your question is to be interpreted in general terms as “How do you define a discrete random variable?” (the illustration provided being useless and only adding to the confusion), or if you didn't realize this, in which case there is no need to define a user distribution because the Poisson distribution already exists.

If you are interested to know how to define a parameterized distribution (assuming, of course, that you don't want to use Poisson 😉 , something I could understand given some general drawbacks concerning the way Maple defines discrete random variable, see below).

Look to this simplified version of the worksheet parameterized_discrete.mw, where a more complete definition of Distribution is made and which provides more details including a comparison to Statistics:-RandomVariable(Poisson(...)) :

>

restart

>	with(Statistics):

A parameterized Distribution

>	dist := proc(a) Distribution( ProbabilityFunction = unapply( piecewise(x >= 0, exp(-a)*a^x/x!, 0), x) ) end proc:

>	X := RandomVariable(dist(alpha)): :-ProbabilityFunction('X') = ProbabilityFunction(X, x);

ProbabilityFunction(X) = piecewise(0 <= x, exp(-alpha)*alpha^x/factorial(x), 0)

(1)

>	sum(ProbabilityFunction(X, x), x=-infinity..+infinity)

(2)

>	X := RandomVariable(dist(6)): :-ProbabilityFunction('X') = ProbabilityFunction(X, x);

ProbabilityFunction(X) = piecewise(0 <= x, exp(-6)*6^x/factorial(x), 0)

(3)

Funny: you can use PDF (even if it is ugly) and lure Maple this way: discrete.mw

An idea...

Answered: sand15 1067

September 12 2025

2 38

@AHSAN

>

restart

>	with(plots): with(plottools): with(Statistics):

>	Z := Sample(Normal(0, 1), 10):

>	A := Array([seq(2+sin(10Pii/15+2*Z[i]), i = 1..5)]):

>	B := Array([seq(2+sin(10Pii/15+3*Z[i]), i = 1..5)]):

>	C := Array([seq(2+sin(10Pii/15+4*Z[i]), i = 1..5)]):

>	p0 := ColumnGraph([A, B, C], legend = ["A", "B", "C"]):

>

NbGroups := 5:
NbConds := 3:

dy := 0.5:

p1 := display(
  transform((x, y, z) -> [x, z, y])(extrude(p0, 0..dy))

  , seq(PLOT3D(CURVES([[0, dy, z], [NbGroups-1+NbConds/(NbConds+1), dy, z]], COLOR(RGB, 0.8$3))), z = 0.25..2.5, 0.25)
  , seq(PLOT3D(CURVES([[0, 0, z], [0, dy, z]], COLOR(RGB, 0.8$3))), z = 0.25..2.5, 0.25)
  , seq(textplot3d([i+(NbConds/2)/(NbConds+1), 0, 0, convert(i+1, string)], align=below, font=[Helvetica, bold, 15]), i=0..NbGroups-1)

  , scaling=constrained
  , orientation=[-66, 76, 2]
  , axis[1]=[tickmarks=[]]
  , axis[2]=[tickmarks=[]]
  , labels=["", "", "H"]
):

p1;

>

CondNames := ["A", "B", "C"]:

if NbConds::odd then
  dx := [seq(NbGroups-1+(NbConds/2)/(NbConds+1) + 1/(NbConds+1)*i, i=-(NbConds-1)/2..(NbConds-1)/2)]
else
  dx := [seq(NbGroups-1+(NbConds/2)/(NbConds+1) + 1/(NbConds+1)*i, i=-NbConds/2..NbConds/2)]
end if:

p2 := display(
  p1,
  seq(textplot3d([dx[i], 0, 0.25, CondNames[i]], align=below, font=[Helvetica, bold, 15], color=white), i=1..NbConds)
):

p2;

Download Example_3D.mw

Beware of the ranges...

Answered: sand15 1067

September 09 2025

1 2

lhs(Eq1) > 0 in the range you provide.

Beware_of_the_ranges.mw

Some ways...

Answered: sand15 1067

September 09 2025

0 0

The first problem you face is that LaTeX is a sotware system for typesetting documents and not a computer algebra system.
Maple syntax relies upon maths while LaTeX syntax is most permissive: this means that if any Maple formula can be translated into a syntactically correct LaTeX expression the converse is not true.
For instance this LaTeX expression (where the "/" sign means, in the context of the paper, "conditionnally to") has no Maple meanings:

Another example taken at random from Goodle image.

So, assuming the LaTeX expression has a Maple sense, which would be the case if the author did some Maple work, converted some extression into LaTeX and imported it into its LaTeX development tool (Overleaf for instance), one may think that the LaTeX formula can be translated in Maple.

The LaTeX case

I assume below that you have the LaTeX source not only its compiled rendering.

The a priori best way:
(1) As you have Maple 2025 you must dispose of the MathML:-FromLaTeX function, see HERE.
As I am not that lucky I can't test it to tell you about its possible limitations.

(2) There exists some online translators (AI Chat for instance) but you must extremely very careful on the translation they provide (specifically if your expressions involve vectors, matrices, derivatives, integral).
To be confident in the online translator you could be tempted to use procced this way

In Maple:
1. write some expression expr
2. Convert it to LaTeX.
3. Copy the output.
In the online translator:
1. Paste the LaTeX expression.
2. Ask for its Maple translation
3. Copy it
In Maple:
1. Paste the experssion the translator provided
2. compare to the original maple expression expr.

(3) An old (2016) Mapleprimes question (note that SnuggleTeX no longer exists).

Nevertheless, even if you can use FromLaTeX, point (2) might be usefull depending on FromLaTeX limitations (if any).

The "no LaTeX source code" case

Without LaTeX source (pdf article for instance) a copy-paste of the expression into a Maple worksheet it possible but it is very likely that you will have to correct it.
For instance the illustration I provided at the top of this reply, copied-pasted in maple document mode gives this

The "hand-written" case

Is obviously even more complex than this latter because a first step to convert a handwritten note into LaTeX seems mandatory (see for instance NotesToLaTeX, MathWrite, or others)

In My Opinion,

Even i this work can seem tedious to you, I believe that converting handwritten notes or LaTeX compiled outputs (not a source code) into Maple expressions is much safer if you do it carefully without looking forthe help of some automatic translator.

Last point:

If you are interesting in writting something like LHS = RHS1 = RHS2 = ....RHSN in Maple, just browse the Mapleprimes site: there has been a recent question about whose @acer replied.

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

1067 Reputation

15 Badges

MaplePrimes Activity

These are answers submitted by sand15

Explanations...

KKT doesn't seem necessary here...

elementwise comparisons...

Suggestion...

Simply...

No solution...

Here it is...

Linearize to avoid non convergence issue...

A solution...

discont...

No solutions in the search domain you us...

No type declaration: Type depends on the...

An idea...

Beware of the ranges...

Some ways...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again