1781 Reputation

14 Badges

4 years, 233 days

MaplePrimes Activity

These are answers submitted by mmcdara

If I'm not mistaken, here is the "first order" aymptotic expansion.
(This could be done by hand using the Laplace's Method)

The result is of the form

G(x) * N*exp(-N*sqrt(x^2+1)+N*arcsinh(1/x))

(G is given below and in the attached file).

Note: this asymptotic expansion goes to 0 if x is larger than a critical value equl to 0.6627 and goes to +oo is x < 0.6627.
This means this asymtotic expansion is wrong (at least) for x < 0.6627 (the integral vanishes for any positive x as N --> +oo)... either one should push the expansion to larger orders or I did something wrong

# the critical value of x verifies

I do not know if Maple can do this without help

f := (n, x) -> Int(exp(-n*(x*cosh(t)+t)), t = 0 .. infinity);

# the previous relation is of the form

f__x(epsilon) = Int(g(t)*exp(phi__x(t)/epsilon), t = 0 .. infinity);

# with 
#   g(t)=1, 
#   phi__x(t)=-(x*cosh(t)+t)   (x is consider as a known parameter)
#   epsilon=1/n

# Rewritten this way one can apply the Laplace's Method:
#     g(t) and phi__x(t) are smooth functions
#     epsilon is a small positive parameter
# Doing this one must assume:
#     1/ that phi__x has a global maximum at some value t=c
#        (cf below)
#     2/ diff(phi__x(t), t$2) < 0 at t=c
#        The positiveness of x implies this condition is verified too 

phi__x := t -> -(x*cosh(t)+t) ;
c      := solve(diff(phi__x(t), t)=0, t);

diff(phi__x(t), t$2);
eval(%, t=c);         # negativeness of diff(phi__x(t), t$2) at t=c

# Laplace's Method:
# replace phi__x(t) by its taylor expansion up to order 2
# to get 

f__x(epsilon) = Int(exp(('phi__x(c)'+1/2*D[2]('phi__x(c)')*(t-'c')^2)/epsilon), t = 0 .. infinity);

h :=convert(taylor(phi__x(t), t=c, 3), polynom);

# verify this is identical to 

phi__x(c) +1/2*eval(diff(phi__x(t), t$2), t=c)*(t-c)^2;


# note that phi__x(c) doesn't depend on t, thus 

f__x(epsilon) = exp('phi__x(c)'/epsilon)*Int(exp(1/2*D[2]('phi__x(c)')*(t-'c')^2/epsilon), t = 0 .. infinity);

# The integral could easily be computed by hand (a ERF function)
# but I use Maple
#    I set A = D[2]('phi__x(c)') (which is negative)
#    and C=c to avoid premature evaluation and keep the expression simple.

K := int(exp(1/2*A*(t-C)^2)/epsilon, t=0..+infinity) assuming A < 0, epsilon > 0;

# Finally

J := simplify( 
       eval(K, [C=c, A=eval(diff(phi__x(t), t$2), t=c)]) 
     ) assuming x > 0;

# and thus the asymtotic expansion with regard to N:

Asympt_J := combine(asympt(eval(J, epsilon=1/N), N), exp);

In the plotsetup help page it is written :
          A list of colors cannot be provided to the plots[pointplot] or plots[pointplot3d] command.

Here is a workaround (Maple 2015)

List := [[[0, 0], 1], [[1, 2], 0], [[2, 3], 0], [[3, 1], 1]]:
C := n -> COLOR(RGB, op(`if`(List[n][2]=1, [1, 0, 0], [0, 0, 1])));

p := PLOT(seq( POINTS([List[i][1]], C(i)), i=1..nops(List) )):
# verify

plotsetup(ps, plotoutput=cat(currentdir(), "/Desktop/test"), plotoptions=`color=cmyk`);

Here is the content of the file test.eps (the extension is automatically added).

As Mapleprimes doesn't accept to load eps files I opened the eps file with a ps wiever and exported it to a png file... but believe me, test.eps contains colors.



I understand your dismay if you are not familiar with signal processing, DFT (Discrete Fourier Transform), and signal sampling.

You will find a few explanations in the attached file but I strongly advice you to read some book or course on the topic (I have a few that I used when I was a student but they are getting old and are in French).

Basically the DFT differs from the FFT on at least two points (I put aside considerations related to the sampling [=discretization] of the signal):

  • signals must be causal (because they are assumed to be real signals!)
  • the signal being observed on a finite period of time (or space), one must decide what this signal is outside of this observation window:
    • either it is considered as perioc
    • either it is considered as aperiodic

These considerations implie the FFT and the DFT geberally give different "raw" results.
But knowing the assumptions used by the DFT one can recover the result a FFT gives.

In your case, knowing these assumptions
(cf formulas used in DiscreteTransforms[FourierTransform] help page)
helps to explain the discrepancy betwen the exact and the numerical solutions, and also enable "adjusting" the numerical one.

Here is a solution

Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

eq__1 := W__mn^2*n__1+L__11*U__mn+L__12*V__mn+L__13*W__mn = C__1*`&Delta;T`:
eq__2 := W__mn^2*n__2+L__21*U__mn+L__22*V__mn+L__23*W__mn = 0:
A, b := GenerateMatrix([eq__1, eq__2], [U__mn, V__mn]);
sol := collect~(LinearSolve(A, b), W__mn):

# print solution
< [U__mn, V__mn] > = sol;

# print coefficients
[seq(a__||k, k=0..2)] =~ coeff~(sol[1], W__mn, [0, 1, 2]);
[seq(b__||k, k=0..2)] =~ coeff~(sol[1], W__mn, [0, 1, 2]);


I don't know what you expect from Maple  (what do you doubt of as the name of your file seems to mean)?
Given the abstract definitions of RV x and beta, you have absolutely no chance of finding (if that's what you were hoping for)  any closed form expression for q__ER nor Pi__ER.

Here is the farthest you can go with Maple (and probably with any other CAS or even by hand)



d := Distribution(
       PDF = (proc (t) options operator, arrow; f__D(t) end proc), 
       Quantile = (proc (q) options operator, arrow; int(f__D(t), t = -infinity .. q) end proc)

x := RandomVariable(d);

F := t -> Quantile(x, t);
F__c := t -> 1- F(t);

t -> Statistics:-Quantile(x, t);
t -> 1 - F(t);


d1 := Distribution(
        PDF = (t -> g__beta(t)), 
        Quantile = (q -> int(g__beta(t), t=-infinity..q))

beta := RandomVariable(d1);

G := t -> Quantile(beta, t);
G__c := t -> 1- G(t);


Pi__ER := beta__R -> piecewise(
            q/(G__c(beta__R)+m*G(beta__R))-x >= 0, 
            `and`(x-q/(G__c(beta__R)+m*G(beta__R)) > 0, q/G__c(beta__R)-x > 0), 
            x-q/G__c(beta__R) >= 0, 

# Differentiate first and take the mean after

dp := diff(Pi__ER(beta__R), q); 

# write the conditions under simpler forms

UV := [ op(1, rhs(op(1, dp))) = U, -op(1, rhs(op(-2, dp))) = V ];

dpUV := eval(dp, UV);
# simplify the writing of dpUV

dpUV := piecewise(_R <= U, op(2, dpUV), _R < V, op(4, dpUV), op(6, dpUV))


mean_dp := Mean(dpUV) assuming U < V:
mean_dp := eval(%, int=Int):
mean_dp := add(rhs~(map(Rule[`c*`], [op(mean_dp)])));

mean_dp := (-c+s)*(Int(f__D(_t)*_t, _t = -infinity .. U))+(p-c)*(Int(f__D(_t), _t = U .. V))+(m*p-m*s-c+s)*(Int(f__D(_t), _t = V .. infinity))

# let make "q" to appear (if you want it)

Krule  := op(1, denom(lhs(UV[2]))) = K;
UVrule := eval((rhs=lhs)~(UV), Krule);
mean_dp_with_q := eval(mean_dp, UVrule);

# It is obvious that q__ER defined by solve(mean_dp_with_q=0, q) cannot be computed 
# for f__D has no algebraic expression.
# So let's assume q__ER is known from any other mean

# Proceed as previously to compute Mean_p = Mean(Pi__ER):

pUV := eval(Pi__ER(beta__R), UV):
pUV := piecewise(_R <= U, op(2, pUV), _R < V, op(4, pUV), op(6, pUV)):

mean_p := Mean(pUV) assuming U < V:
mean_p := eval(%, int=Int):
mean_p := add(rhs~(map(Rule[`c*`], [op(mean_p)])));


mean_p__ER := eval(mean_p, q=q__ER);

I've just added a few "partially instanciated" examples at the end of the attached file to (try to) convince you that there is no hope in finding a formal expression  of q__ER and Pi__ER  until f__D and g__beta are abstrat



As long as you have not defined the PDF or CDF of D (that you must declare as local as D is a protected name) you will not be able to find the expression of q.

In the attached file you will see how far you can go FORMALLY with an abstract random variable D (that is a random variable for wich the PDF is just f__D(t) [I assume  D is continuous and f__D exists] ).
At the end you will find 4 examples where solve (..., q) doesn't fail.

local D:
d := Distribution(PDF = (t -> f__D(t)), Quantile = (z -> int(f__D(t), t=-infinity..q)));
D := RandomVariable(d);

TC := piecewise(q-D >= 0, c__o*(q-D), q-D < 0, c__u*(-q+D));

M := Mean(TC);
diff(%, q):
eval(%, op([2, -1], %)=1-op([1,-1], %));
IntegralTerm := op([1,-1], %);
eval(%%, IntegralTerm = J);
isolate(%, J);
expr := eval(%, J=IntegralTerm);

'Quantile(Q, q)' = Quantile(D, q);

# here is the rule to rewrite the lhs of expr as the value of the CDF of D at quantile q

MyRule := q -> subs(t=lhs(op(-1, lhs(expr))), Quantile(D, q)) = F__D(q):

Expr := eval(expr, MyRule(q));

solve(Expr, q);

# Example 1:

eval(Expr, F__D(q) = Quantile(Uniform(0, 1), q));
solve(%, q);
                          c__o + c__u

# Example 2:

eval(Expr, F__D(q) = Quantile(Exponential(lambda), q));
solve(%, q);
                    /          c__u        \    
                -exp|- --------------------| + 1
                    \  lambda (c__o + c__u)/    

# Example 3:

eval(Expr, F__D(q) = Quantile(BetaDistribution(a, b), q));
solve(%, q);

   /   c__u    \           /                        c__u    \
   |-----------|  hypergeom|[a, 1 - b], [1 + a], -----------|
   \c__o + c__u/           \                     c__o + c__u/
                          Beta(a, b) a                       

# Example 4:

eval(Expr, F__D(q) = Quantile(Normal(mu, sigma), q));
solve(%, q);

                /                            (1/2)\    
           1    |(c__o mu + c__u mu - c__u) 2     |   1
         - - erf|---------------------------------| + -
           2    \      2 sigma (c__o + c__u)      /   2

A few errors :

  • plot
    • the stntax is not correct: a range is needed for t
    • M is still undefined: you need to set it to a numerical value before plotting ; for instance
      plot(eval( dyr(t), M = 1), t=0..1);
  • dyr(5) cannot work because the expression of dyr(t) contains derivatives with respect to t: what would it mean to "derive with respect to 5"?
    eval(dyr(t), t=5)
              /               (1/4)                    (3/4)
            M \-0.6609839668 4      + 0.0005882438542 4     
               - 0.0001114078453 4     /
  • This is basically the sme problem for z(1) ; do this instead
    z := unapply(diff(dyr(t), t), t):
             /                (1/4)                    (3/4)
           M \-0.07155390824 4      + 0.0001061325461 4     
              - 0.00002412060546 4     /

As tomleslie opened the gate, here is a model which is infered from the observation of the data and which is aimed to respect some of their feautures:

  • Y is (seems to be) positive
  • Y is monotonically increasing (unless if Y[298]=Y[297]-1-

More of this a close examination of your data leads me to infer there are two adjacent X ranges:

  • for X fro 1 to (about) 220, Y evolves in some way
  • for X > 220 the evolution of Y seems to be the same as in the previous range, maybe just shifted and scaled in a propper way

A model that respect these requirements and observations is

th  := x -> a__1*(tanh((x-c__1)*b__1)+1) + a__2*(tanh((x-c__2)*b__2)+1):

If you discard the positivity constraint you can use this more general model

th := x -> a__1*(tanh((x-c__1)*b__1)+d__1) + a__2*(tanh((x-c__2)*b__2)+d__2)

Here is the result for this latter

Other solutions (better or worse) can be obtained by using similar models:

  • replace tanh by arctan
  • replace tanh by a sigmoid function

As a rule "before using a nonlinear fir, always look if the model is not linear under a proper transformation".
This will potentialy avoid numerical problems.
In your case the transformation is obvious:

y = A*exp(r*x) --> log(y) = log(A) + r*x


beta := LinearFit([1, x], X, log~(Y), x, output=parametervalues);

Model := x -> exp(beta.<1, x>):

    117.5619317177743 * exp(0.02875273320702703 * x)

  ScatterPlot(X, Y), 
  plot(Model(x), x=(min..max)(X), legend=evalf[4](expand(Model(x)))), 
  view=[default, 0..Model(X[-1])]

The fit is nevertheless rather poor.
So let's use this fit as a starting point and look if one can do better?
IMO Fit is not very efficient and I prefer doing the things "by hand" and to use
Optimization:-Minimize instead:

Prior := x -> A*exp(x*r):

# Objective function 

J := add((Y - Prior~(X))^~2):

opt := Optimization:-Minimize(J, initialpoint = {A=110, r=0.03});
  [3.26776e10  , [A = 10049.7774, r = 0.00735]]

# For information, this is the value of J for Model(x)
J := add((Y - Model~(X))^~2);

Post := eval(Prior(x), opt[2]);

  ScatterPlot(X, Y), 
  plot(Post, x=(min..max)(X), legend=evalf[4](expand(Model(x)))), 
  view=[default, 0..Model(X[-1])]

Even if Post is better than Model for it gives an objective function 100 times smaller, I'm not sure that we can't do even better by choosing different initial values for A and r...

For an unknown reason DrawGraph transforms each weight w as convert(evalf(w), string) and generates a TEXT structure for each weight which looks like this one

# Here is how the weight on the edge from not A to C is displayed
# Note that 1/7 has been transformed as "0.143"
op(-5, P)
    TEXT([0.7752048164, 0.2971539501], "0.143", ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))

whereas what is expected is

TEXT([0.7752048164, 0.2971539501], 1/7, ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))

Knowing this is quite easy to recover the initial form

# construct the graph WITHOUT weights displayed 

P__unweighted  := [ op(DrawGraph(G, style=tree, root=Omega, showweights=false)) ]:

# Extract the part of the plot that concerns the weights

EdgeWeights    := map(u -> if `not`(u in P__unweighted) then u end if, [ op(P) ]):

# Recode the weights as rationals

NewEdgeWeights := map( u -> subsop(2=convert(parse(op(2, u)), rational, 3), u), EdgeWeights):

# Reassemble the "weighted" structure

P__weighted    := P__unweighted[], NewEdgeWeights[]:

# Plot it
# as in my answer to your initial question I syill use Maple 2015, 
# so the rendering might be slightly different from yours

HighlightVertex(G, Vertices(G), white);
plots:-display(R(PLOT(P__weighted)), size=[600, 600]);

PS: I suspect the rational to float transformations are justified to avoid complicating the algorithm wich calculates vertex positions.
(the weight information "takes" on line instead of 2).
In my opinion a transformaton of the form 

r := 1/7:
rr := numer(r) %/ denom(r)

whould be better

PS: As I said you while answering your first question about probability graphs, you cann replace a vertex name like A by defining  A:=`#mo(A)` in order to print all the vertices with the same font.

I understand from your previous question that you are familiar with LaTeX?
In this case you could be find a few ideas here

Basically I generate the LaTeX code of a "tabular" structure directly in Maple.
The only prerequisite is yo know LaTeX.
I'm responsible of a computational code written in Maple and I have several test cases I need to run periodically to insure there is no code regression. I use to use the trick described in the attached file to generate automatically "execution reports" and insure traceability.

PS: I did not compile the tex cource to verify I didn' make mistakes, but I think I'm right


It's a consequence of the (rather comple) algorithm implicitplot uses (type showstat(`plots/implicitplot`)  in a worksheet and look at it).
In short, if numpoints if specified, implicitplot begins by computing a grid whose dimension will be nx by ny, where 

nx, ny := `$`(isqrt(numpoints+3)+1,2)

(for numpoints=1000 this gives a 33 by 33 grid).
This is the grid the curve is interpolated on.
The resolution intervenes later to determine the number of point where F(x, y) is to be evaluated.

For gridrefine=3 F(x,y) is calculated at 515 points, while 4119 points are used for gridrefine=6.
(Note these numbers are twice the numbers I obtained by examining the output of showstat(`plots/implicitplot`) ).

Now there is a third option named resolution, wich plays a role when explitely set to a value its default value is 0).
By setting resolution = k times 'the number of points where F(x,y) is computed' (k integer > 0) you can force implicitplot to really use all the points (and more) to get a ... better resolution of the curve.

In the attached file:

  • p3 is the plot obtained with gridrefine=3,
  • p6_0 is the plot obtained with gridrefine=6, and resolution=0 (default value),
  • p6_4 is the plot obtained with gridrefine=6, and resolution=4*4119 (or about).

A remark : 

M[i, j] := M[i, j] + something

# is useless, it's enough to write

M[i, j] := something

This is a shorter way to obtain the matrix M:
(maybe "expr" is less readable)


A := Matrix(2, 3, [[4, x, 1/4*x^2], [2, 3, x]]);
B := Matrix(3, 2, [[3, 5], [x, 2], [1/2*x^2, -x]]);

a := LinearAlgebra[ColumnDimension](A):
b := LinearAlgebra[ColumnDimension](B):

A := unapply(A, x):
B := unapply(B, x):

alpha := k -> x -> A(x)[k, ..]:
beta  := k -> x -> B(x)[.., k]:

expr := (i,j) -> alpha(1)(1/2)[i] * beta(1)(1/2)[j]
                 alpha(2)(1/2)[i] * beta(2)(1/2)[j]
                 int(alpha(1)(x+1/2)[i] * beta(2)(1/2*x)[j], x =-1/2..1/2);

M := Matrix(a$2, (i, j) -> expr(j, i));

The first example can be proved his way

CompleteSquare(a^2+2*a*b+b^2, a);
CompleteSquare(a^2+2*a*b+b^2, b);
                            (b + a) 
                            (b + a) 

Be patient for the second...

It works correctly, maybe you have made a false manipulation

Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

diff(y(x), x, x) = 1/y(x) - x*diff(y(x), x)/y(x)^2;
sol := dsolve({%, y(0)=1, D(y)(0)=0}, numeric);
plots:-odeplot(sol, [x, y(x)], x=0..1)

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