4019 Reputation

17 Badges

6 years, 177 days

MaplePrimes Activity

These are replies submitted by mmcdara

@Carl Love 

My answer was based on my understanding, of what I thought the OP wanted.
Regarding the dualaxisplot, I totally agree with you: in this case, it looks more like comparing two graphs and a simple display seems more appropriate.


I have thought it would be obvious, sorry:

plots:-display(`<|>`(TwoMat(A_jk, A_ki), TwoMat(A_jk, A_ki, scale = true)));

# Or if you just want unscaled matrices:
TwoMat(A_jk, A_ki)

Maybe you are puzzled by the plot if A_ki ???
In this case look to this matrix or do 

plots:-matrixplot(A_ki, heights = histogram)

to understand what happens


  1. I need to think again to that.
    "minimizing only finds one of the solution": right, but I believe that there is (maybe I'm wrong) a single set of values for the lambda(s) and an infinity for the mu(s); something as if the solution was a variety of dimension 1 in a space of dimension 8 (???).
    Observation: starting from different initial points didn't give me different results for the lambda(s), only for the mu(s).
  2. Observation: It's clear that 6 equations depend only on the lambda(s), but applyng fsolve on this sub system gives no solution whatever the value of Digits or the initial point (for those I used).
  3. By the way V=~0 is not a solution:

Thanks for all your comments


Optimization:-NLPSolve is faster and more stable than fsolve:

Splitting the system Eqp in two sub-systems reveals that its solution verifies 

{mu__jk = -mu__ki, mu__ki = mu__ki}


The dyadic product can be realized in a simple way without using LinearAlgebra:-KroneckerProduct:

N := 3: 
A := Vector(N, symbol=a): 
B := Vector(N, symbol=b): 
M1 := LinearAlgebra:-KroneckerProduct(A, B^+): 

# alternative way

M2 := Matrix(N$2, (i, j) -> A[i]*B[j]):

Which method is the fastest?



N := 3:
A := Vector(N, symbol=a):
B := Vector(N, symbol=b):
LinearAlgebra:-KroneckerProduct(A, B^+):
M := Matrix(N$2, (i, j) -> A[i]*B[j]):

%%, %

Matrix(3, 3, {(1, 1) = a[1]*b[1], (1, 2) = a[1]*b[2], (1, 3) = a[1]*b[3], (2, 1) = a[2]*b[1], (2, 2) = a[2]*b[2], (2, 3) = a[2]*b[3], (3, 1) = a[3]*b[1], (3, 2) = a[3]*b[2], (3, 3) = a[3]*b[3]}), Matrix(3, 3, {(1, 1) = a[1]*b[1], (1, 2) = a[1]*b[2], (1, 3) = a[1]*b[3], (2, 1) = a[2]*b[1], (2, 2) = a[2]*b[2], (2, 3) = a[2]*b[3], (3, 1) = a[3]*b[1], (3, 2) = a[3]*b[2], (3, 3) = a[3]*b[3]})


N := 1000:
A := Vector(N, symbol=a):
B := Vector(N, symbol=b):
CodeTools:-Usage( LinearAlgebra:-KroneckerProduct(A, B^+) ):
CodeTools:-Usage( Matrix(N$2, (i, j) -> A[i]*B[j]) ):

memory used=91.63MiB, alloc change=339.64MiB, cpu time=3.36s, real time=2.85s, gc time=1.65s

memory used=45.79MiB, alloc change=7.63MiB, cpu time=876.00ms, real time=615.00ms, gc time=356.94ms




As it sometimes happens, the built-in procedure is neither the fastest nor the most parsimonious in terms of memory allocation.



A faster algorithm to plot the level curve f(a, b)=C.
(same toy problem than in my previous answer).



hA := 0.2:
hB := 0.2:

alpha := Array([seq(0..10, hA)]):
beta  := Array([seq(0..10, hB)]);
PA := numelems(alpha);
PB := numelems(beta);

beta := Vector(4, {(1) = ` 1 .. 51 `*Array, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})






NN := Matrix(PA, PB, (i, j) -> cos(alpha[i]/3+beta[j]/3)*sin(alpha[i]*beta[j]/9))

NN := Vector(4, {(1) = ` 51 x 51 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})


phi := (a, b) -> [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)]

proc (a, b) options operator, arrow; [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)] end proc


C  := 0.2:
for i from 1 to PA-1 do
  for j from 1 to PB-1 do
    nodes  := NN[i, j], NN[i, j+1], NN[i+1, j+1], NN[i+1, j]:
    bounds := (min..max)(nodes):
    if verify(C, bounds, interval) then
       loc := add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [nodes]):
       s   := solve(loc, b):
       LC  := LC,
                add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [nodes]),
    end if:
  end do:
end do:

display(LC, view=[0..10, 0..10])






Thznk you acer


I take this limit as a clever and practical alternative that, for the time being, enables me to go further.
Nevertheless I would really understand why dsolve does the things differently depending on the way the system is writte:

  1. Why do sys_2 and sys_4 provide different results?
    I thought that sys_2 was automatically rewritten ender the form sys_4?
  2. Why do sys_2 and sys_5 give different expressions?
    Is it that having the term 
    in sys_5 precludes early simplifications?

What is disturbing is that dsolve knows how to build the the exact solution of very complex equations and, as in this case, seems unable to put the final touch


Your 2015 version will help me while telecommuting (right now) and I will use your 2022 as soon as I get back to the office.
Thank again

If I'm not mistaken I think that I asked a question about a progress bar, maybe three years ago, and that you gave a very satisfying answer.
Strangely I'm not able to put a finger on this question and I don't remember where I saved your original answer.
Nevertheless I keep using it under different forms, in particular:

MV := parse(Trim(StringSplit(Split(interface(version), ",")[2], "Maple")[2])):

P := NULL:
for i to 22 do
  P := P, plot([[i, i^.5]], style = point):
  Tlayout := Table(Column(), widthmode=percentage, width=40, alignment=left, interior=none ,
    Row( Cell( Textfield( InlinePlot(plots:-display(P, title=typeset(Iteration=i))) ) ) )

  if MV > 2018 then
    InsertContent(Worksheet( Tlayout ) )
    InsertContent( Worksheet( Group( Input( Tlayout ) ) ) )
  end if:
end do:

I think that using DocumentTools[Layout] enables more customization than Tabulate.

Sorry for the misunderstanding
And thanks for your clarification

A few remarks:

  • Your problem is not completely defined:
    • Let Sn the event "individual xn gets sick"; this event has a probability p to occur.
      But are the events Sn and Sm independent?
      More generally, let x={x1, ..., xN} a set of N such individuals and Prob( S1)=...= Prob( SN)=p: are the events S1, ..., Smutually independent?
    • The probability to get sick generally depends on the period of time on considers.
      You make no mention of that.
  • In order to model correctly your real problem it(s necessary to firstly answer the two the questions above.


A probabilistic model, but not the only one, could be this one:

  • Assuming the probability to get thick doesn't depend on any time period, the event Sn is the outcome of Bernoulli random variable Ber(p).
  • Assuming mutual independence, the event "K individuals out of N get sick" is the outcome of a Binomial random variable Bin(N, p).
    (see line 1).
    Remark: mutual independence is a necessary and sufficient condition for the sum of N Bernoulli random variable of identical parameter p be a Binomial random variable, for a counter example see here (for instance)).

Given this model Statistics:-Sample(Binomial(100, .3), 1))  returns a realisation of Bin(N=100, p=0.3), that is a value of K (a float in fact, thus the truncation used by Carl Love).


PS instead of using 
[seq](trunc(x), x= Statistics:-Sample(Binomial(100, .3), 9));
[seq](Student[Statistics]:-Sample(BinomialRandomVariable(100, .3), 9));
(Student[Statistics]:-Sample(...) returns an integer as Statistics:-Sample(...) returns a float)


The case of "perfect" dependence:

  • One still assumes the probability to get thick doesn't depend on any time period.
  • If everybody gets thick instantaneously as soon as a single gets thick (perfect dependence) then the event "K individuals out of N get sick" has probability 0 if K < N and 1 if K=N.

For non-dependent Bernoulli random variables you can see here:


@acer @tomleslie

Thanks to both of you.
I've been on vacation since yesterday and will be back at the office next Monday, so I can't give you, right now,  the exact version of Maple I used.
All I can say is that it was Maple 2020. ???? on Windows 10.
I will come back to this thread on Monday with more information.


Very disturbing last figure.
Stacking histograms is quite unusual and gives a very wrong idea of the empirical distributions.


Thanks for your reply

@Shameera : the answer of your deleted answer is avaliable here  233822-How--To-Solve-This-Equation-By-Using

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