5284 Reputation

17 Badges

7 years, 121 days

MaplePrimes Activity

These are replies submitted by mmcdara

I agree with @Rouben Rostamian : some clarification is needed.

A few remarks:

  • Do not load useless packages, this is the best way to generate conflicts.
  • Do not try to use odeplot if you didn't use successfuly dsolve/numeric before.
  • If you use dsolve/numeric, you must provide initial and/or boundary conditions.
  • Before using a package, CurveFitting for instance, look at the corresponding help page to see what procedures it contains instead of imagining them (NonlinearFit is part of the Statistics package).

You will find in the attached file some elements which lead to a solution of the differential system.
I took some liberties in choosing arbitrary the initial conditions and the values of the parameters (which I all rewrote in, IMO, a more convenient form).

I discuss your "fit to data problem", but this is not conclusive at all: you need to be more precise in your request.


@Rouben Rostamian : I hope I didn't step on your toes.


but it"s likely there are som conflicts:

  • between 2D input mode and 1D input more,
  • the Maple versions (I use Maple 2015).
    When I open your file I get a warning message, something like " this code has been written in a newer version and could not run correctly".
    If I agree to "activate" the worksheet, then when I save it under another name, the m file is a Maple 2015 mw file.
    When you will open it you won't get any warning on your side. But the fact is hat there is also a version mixing between 2023 and 2015.

Here is corrected, at least in the sense that the error is removed.



That's would be to much work for me, not to mention the fact that writing this code requires an even more description of the DNSE method.
Hope you will find someone else to help you.


when you write "Can someone guide me on how to write the code for this program?", what do you mean exactly?

  • Do you want to write a code which "proves" the equations (14) given equation (9) and all the stuff in between?
    if it is so, here is a code
  • Or do you want to solve equations (14) and prove a particular solution as the form you give?
  • Or maybe something else?


I'm going to look more closely to the paper.


My mistake, the determinaant is 


and I read it 


when I wrote my re^ly.

So the determinant is always negative and the hessian matrix is never positive definite.

To check if it negative definite the first minor of this matrix must be negative.
Look at the end of this file

So, the hessian matrix never being positive  definite, there is no point in R3 at which TPpm reaches a minimum.

relative to the docs reference:

  • What happened to index n in equation (4) ?
    As you write it the formula means nothing.
  • In equation (9), is "i" the imaginary unit?
  • Next line, what does DNSE mean?
  • How can you transfor the single equation (9) into a system of 2 equations (12)?
    Is it assumed that u(t) = UA(t)+I*V(t)  and that phi(t) = P(t)+I*Q(t) with both A, B, P, Q functions with values in IR ?
    In which case (12.1) could be Re(eq(9)) and (12.2) Im(eq(9)) ?

Can you provide a more solid reference for what you want to achieve for, at least for me, it is incomprehensible.


You never calculate the Hessian of TPpm wrt Q and Qr.
You define a, b, c this way

a := diff(TPpm(Q, Qr), Q$2):
b := diff(TPpmg(Qr), Qr$2):
c := diff(TPpm(Q, Qr), Q, Qr):

On the other way TPpmg(Qr) is given by

TPpmg := Dr -> subs([g1], TPpmm(Q, Qr))

and then IS NOT identical to TPpm(Q, Qr).

Then the matrix 

Matrix(2, 2, {(1, 1) = a, (1, 2) = c, (2, 1) = c, (2, 2) = b})

is not the hessian matrix of TPpm(Q, Qr).

The Hessian of TPpm(Q, Qr) is 

H, det := Student:-VectorCalculus:-Hessian(TPpm(Q, Qr), [Q, Qr], determinant):

H := map(simplify, H);

H := Matrix(2, 2, {(1, 1) = (D-P)*Ch/(P*D), (1, 2) = -(1/2)*Ch/D, (2, 1) = -(1/2)*Ch/D, (2, 2) = 0})






Then the conditions for H to be definite positive are:

  • det > 0 ==> C < 0
  • H[1, 1] > 0 thus
    solve({H[1, 1] > 0) assuming C < 0:
                         {0 < Ch, 0 < P, D < 0}
                         {0 < Ch, 0 < P, P < D}
                         {0 < D, Ch < 0, P < 0}
                         {Ch < 0, D < P, P < 0}
                     {0 < Ch, D < 0, P < 0, P < D}
                     {0 < D, 0 < P, Ch < 0, D < P}


Complete file here




A simple way is to insert a 


command BEFORE the portion you want to debug.
This command MUST be located AFTER the last local and global declarations.

This means that instructions like

local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}:

cannot be debugged.

You must change the code this way



TriangularTypeII := proc(
                        P::list(numeric):=[0.25, 0.75],
   description "Q is a pair of quantiles associated, by default, to probabilities 1/4 and 3/4 (in the case you have other values for these probabilities, use a third argument P to mention them, c is the mode";
   local TRI, F, eqs, AB:

   if debugmode then DEBUG() end if:

   TRI := RandomVariable(Triangular(a, b, c)):
   F   := unapply(CDF(TRI, t), t):
   eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}:
   AB  := solve(eqs, {a, b}):

   Specialize(TRI, convert(AB, list));
end proc:

# the sequel is in the attached file



When you run the command

T := TriangularTypeII([3, 5], 4):

the debugger is not invoked.
But if you run this one

T := TriangularTypeII([3, 5], 4, debugmode=true):

the debugger opens a pop window in which you can execute step by setp the procedure and see what happens.

Type help(DEBUG) for more information.

Another way I'm less familiar with is to use the Maple debugger directly from your worksheet.
Type help(debugger) for more information.

As I use to use interactive debuggers in other languages, I prefer using DEBUG() in the way I present in the attached file (this mimic a compilation directive).
But I'm not competent enough to say which of the two methods is more efficient.


Let me know if you have any problems.

Maple offers the possibility to define optional arguments for a procedure.
Optional arguments are gathered in a set of the form

{OptArg1::type:=DefaultValue1, ...OptArg?::type:=DefaultValueN}

If the procedure is regularly called with an argument which has a given value, you can declare this argument as optional and define this value as its default value.
In the attached file I assumed that  you generally call TriangularTypeII with Q=[Fitst Quantile, Third Quantile]:

Simplifyng g1 seems a dead end to me.
Getting rid of the "1." you mention is easy.

Can you precise your last question?

Minor1 := simplify(H[1, 1]);
Minor2 := d;
              /    2      3      2            2     3\
           Ch \P Qr  T - Q  - 3 Q  Qr - 3 Q Qr  - Qr /
         - -------------------------------------------
                           P (Q + Qr)                 
                              2    2  
                            Ch  T Q   
                         - -----------
                           P (Q + Qr) 
# H is positive definite if Minor1 and Minor2 are both > 0
# If H is positive definite at a given point (Ch, P, Q, Qr, T) then 
# this point is a minimum of TPpm.



How Maple realize that the pairs of Q and P are quantiles and their probabilities (the 2 pairs) in this equation:  local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}

This is the definition I adopt when I build the distribution: the first argument Q is a list of 2 quantiles and the second P is the list of their associated probabilities.
I hesitated between this definition and this one:

  • the first argument is a couple (quantile, probability), let's say A,
  • the second argument is another couple (quantile, probability), let's say B.

(or A and B are 2 couples (probability, quantile), maybe)

Whith this choice

eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}

would have been written

TriangularTypeII := proc(A, B, c)
eqs := {F(A[1])=A[2], F(B[1])=B[2]}

By the way, a safer construction of the calling arguments is

TriangularTypeII := proc(A::list(numeric), B::list(numeric), c::numeric)

This construction returns an error if you call TriangularTypeII with formal arguments (see my second reply), for instance

TriangularTypeII([q1, p1)], [q2, p2], c)

produces an error as soon as at least one of these five quantities is not a number.

Here is an example where I used the option description as an operand of the procedure (see help(procedure)) and the function Describe to get some information about this procedure.



TriangularTypeII := proc(Q::list(numeric), P::list(numeric), c::numeric)
   description "Q is a pair of quantiles, P the pair of their associated probabilitues, c the mode";
   local TRI := RandomVariable(Triangular(a, b, c)):
   local F   := unapply(CDF(TRI, t), t):
   local eqs := {F(Q[1])=P[1], F(Q[2])=P[2]}:
   local AB  := solve(eqs, {a, b}):
   Specialize(TRI, convert(AB, list));
end proc:


# Q is a pair of quantiles, P the pair of their associated probabilitues, c the
# mode
TriangularTypeII( Q::list(numeric), P::list(numeric), c::numeric )


test := TriangularTypeII([q, 5], [1/4, 3/4], 4)

Error, invalid input: TriangularTypeII expects its 1st argument, Q, to be of type list(numeric), but received [q, 5]


T := TriangularTypeII([3, 5], [1/4, 3/4], 4)






At the very end, if you have several reparameterized distributions to build, or new ones to construct, it worth wtitting help pages for each of them (see help(Templates)) and gather all these distribution in a specific package.


Thaks for your kind reply.

As you probably understood the trick (as in the definition of the SubjectiveBeta distribution --at least the revised version of my answer--) is to transform the 5 parameters (Q, P, c) into 3 parameters (a, b, c), and then Use the Triangular(a, b, c) distribution.

Based on the use of the function Statistics:-Specialize, there is a shorter version of my previous file.

It gives exactly the same results than the first version, but I don't like it because the TriangularTypeII distribution is not really defined: what I do here is, starting from an abstract Triangular(a, b, c) distribution, to build a triangular random variable whose parameters are no longer  (a, b, c) but (Q, P, c).
Statistics:-Specialize plays a key role here but, I repeat, the TriangularTypeII distribution doesn(t really exist.

About your "4 solution issue".
There is an explanation:



TRI := RandomVariable(Triangular(a, b, c)):
F   := unapply(CDF(TRI, t), t):


piecewise(q < a, 0, q <= c, (q-a)^2/((b-a)*(c-a)), q <= b, 1-(b-q)^2/((b-a)*(-c+b)), 1)


# Each of the 2equalities you must solve are of the form

F(q) = p

piecewise(q < a, 0, q <= c, (q-a)^2/((b-a)*(c-a)), q <= b, 1-(b-q)^2/((b-a)*(-c+b)), 1) = p


# Assuming you are cautious enough not to define a quantile q such that q < a or q > b,
# you can simplify the relation above by considering only the case q in [a, b].
# More than this it will be useful to rewrite this new relation by splitting it into
# a "left relation" corresponding qo q < c and a "right relation" when q > c.

LeftRel  := (q, p) -> op(4, F(q)) = p:
RightRel := (q, p) -> op(6, F(q)) = p:

# You then have 3 sets of relations to consider

LL_eqs := {LeftRel (Q[1], P[1]), LeftRel (Q[2], P[2])};  # Q[1] < c & Q[2] > c
LR_eqs := {LeftRel (Q[1], P[1]), RightRel(Q[2], P[2])};  # Q[1] < c & Q[2] > c
RR_eqs := {RightRel(Q[1], P[1]), RightRel(Q[2], P[2])};  # Q[1] > c & Q[2] > c

{(Q[1]-a)^2/((b-a)*(c-a)) = P[1], (Q[2]-a)^2/((b-a)*(c-a)) = P[2]}


{(Q[1]-a)^2/((b-a)*(c-a)) = P[1], 1-(b-Q[2])^2/((b-a)*(-c+b)) = P[2]}


{1-(b-Q[1])^2/((b-a)*(-c+b)) = P[1], 1-(b-Q[2])^2/((b-a)*(-c+b)) = P[2]}


# Here are their solutions

LL_sol := solve(LL_eqs, {a, b}):
LR_sol := solve(LR_eqs, {a, b}):
RR_sol := solve(RR_eqs, {a, b}):

# In your case (Q[1]=3) < (c=4) < (Q[2]=5).
# Thus select the LR_sol solution:

LR_all := [allvalues(LR_sol)]:

eval(LR_all, [Q[1]=3, Q[2]=5, P[1]=1/4, P[2]=3/4, c=4]):
LR_ab := evalf(%);
LR_ab := map(s ->(lhs=Re@rhs)~(s), LR_ab);  # these are the 4 solutions you talked about



[{a = 3.548583769+0.*I, b = 6.215250432+0.*I}, {a = 3.414213563+0.*I, b = 4.585786412+0.*I}, {a = 1.784749562+0.*I, b = 4.451416224+0.*I}, {a = .5857864375-0.*I, b = 7.414213560-0.*I}]


[{a = 3.548583769, b = 6.215250432}, {a = 3.414213563, b = 4.585786412}, {a = 1.784749562, b = 4.451416224}, {a = .5857864375, b = 7.414213560}]


# The fact you get 4 solutions is disturbing given TriangularTypeII returns only this one.

evalf([2-2^(1/2), 6+2^(1/2)])

[.585786438, 7.414213562]


# But, if look carefully to these 4 solutions, you will see that not all of them
# are correct because they must verify
#     a <= min(Q[1], Q[2])
#     b >= max(Q[1], Q[2])
# Trying to build LR_sol while accounting for these conditions seems to be
# a dead end.

# LR_sol_cond := solve({RR_eqs[], a <= Q[1], a <= Q[2], b >= Q[1], b >= Q[2]}, {a, b}): # dead end

# So I think it's better to select the solutions which verify the conditions above


Q := 3, 5:

map(s -> if
            verify(eval(a, s), -infinity..min(Q), interval)
            verify(eval(b, s), max(Q)..+infinity, interval)
         end if,

[{a = 3.548583769, b = 6.215250432}, {a = 3.414213563, b = 4.585786412}, {a = 1.784749562, b = 4.451416224}, {a = .5857864375, b = 7.414213560}]


[{a = .5857864375, b = 7.414213560}]







f := (b, e) -> 47.965-(-(3.08*0.2*b^(1-0.7))/(GAMMA(3-0.7)))*e

proc (b, e) options operator, arrow; 47.965-(-1)*3.08*.2*b^(1-.7)*e/GAMMA(3-.7) end proc


plot3d(f(z, e), z=1..100, e=-1..1, style=surface, color=blue)


      f(b, e)
      , b=1..100
      , legend=e
      , color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
    , e in [seq(-1..1, 0.5)]




3 4 5 6 7 8 9 Last Page 5 of 113