mmcdara

7891 Reputation

22 Badges

9 years, 55 days

MaplePrimes Activity


These are replies submitted by mmcdara

@delvin 

Hi,

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.

@delvin 

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  From_9_to_14.mw
     
  • Or do you want to solve equations (14) and prove a particular solution as the form you give?
     
  • Or maybe something else?

@delvin 

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

@MANUTTM 

My mistake, the determinaant is 

-(1/4)*Ch^2/D^2

and I read it 

-(1/4)*C*h^2/D^2

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 q2b_mmcdara.mw

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.
TIA

@MANUTTM 

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

(1)

det

-(1/4)*Ch^2/D^2

(2)

 


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:
    print~({%}):
                         {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 q2_mmcdara.mw

@NanaP 

Thanks

@NanaP 

A simple way is to insert a 

DEBUG():

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

restart:

with(Statistics):

TriangularTypeII := proc(
                      Q::list(numeric),
                      c::numeric,
                      {
                        P::list(numeric):=[0.25, 0.75],
                        debugmode::boolean:=false
                      }
                    )
   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

 

Download TriangularDistribution_V2d.mw

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.

@NanaP 

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]:
TriangularDistribution_V2c.mw

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

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 /
         - -------------------------------------------
                                     3                
                           P (Q + Qr)                 
                              2    2  
                            Ch  T Q   
                         - -----------
                                     3
                           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.

@NanaP 

 

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.
 

restart:

with(Statistics):

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:

Describe(TriangularTypeII)


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

_R0

(1)

 

 

Download TriangularDistribution_V2b.mw

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.

@NanaP 

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.
TriangularDistribution_V2.mw

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:

restart:

with(Statistics):

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

 

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

(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

(2)

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

(3)

# 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)]:
numelems(LR_all);

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

4

 

[{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}]

(4)

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

(5)

# 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

LR_ab;

Q := 3, 5:

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

[{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}]

(6)

 

Download TooManySolutions.mw

@Aung 

NULL

restart:

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

(1)

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

 

plots:-display(
  seq(
    plot(
      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)]
  )
);

 

 

Download DoYouWantThis.mw

@Anthrazit 

You should have attached this file when you asked your question, this would have saved me from taking off in the wrong direction and wasting my time.

@Anthrazit 

That doesn't change anything: your code still runs without errors with Maple 2015:

restart:

kernelopts(version)

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

(1)

Segment2Arrow := proc(segm)

    uses geometry;

    description "Annotation dimension";

    local a, b, c;

 

    try

        a := map(coordinates, DefinedAs(segm))[1];

        b := map(coordinates, DefinedAs(segm))[2];

    catch "wrong type of arguments":

        Alert(cat("Error, (in geometry:-DefinedAs) wrong type of arguments, variable segm: ", whattype(segm), segm), table(), 5)

    end try;

    c := (a + b) / 2;

    return plots:-arrow(c, [a - c, b - c], width = 0.5, color = grey)

end proc:

geometry:-point(A,0,0):
geometry:-point(B,1,1):
geometry:-segment(segm, [A,B]):

Segment2Arrow(segm)

 

 

Download 2015_a.mw

Maybe you should add this line after the last local declaration

print(detail(segm));

to see what segm really is.
You should get this`

[name of the object          segm
form of the object           segment2d
the two ends of the segment  [[0,0],[1,1]]
First 44 45 46 47 48 49 50 Last Page 46 of 154