mmcdara

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.

@MaPal93 

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

@dharr 

  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:
V=~0;
eval(Eqp,%);

Thanks for all your comments

@MaPal93 

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

matrix_inverse_mmcdara.mw

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?

 

restart:

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

(1)

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

 

 

Download DyadicProduct.mw

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

 

@Dkunb 

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

restart:

with(plots):

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

 

51

 

51

(1)

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

(2)

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

(3)

C  := 0.2:
LC := NULL:
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,
              contourplot(
                add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [nodes]),
                a=alpha[i]..alpha[i+1],
                b=beta[j]..beta[j+1],
                contours=[C]
              )
    end if:
  end do:
end do:

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

 

 


 

Download LevelCurve_2D.mw

@acer 

Thznk you acer

@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 
    (M__P+C)*sin(t)
    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

@acer 

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

@acer 
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:
 

restart
with(DocumentTools):
with(DocumentTools[Layout]):
with(StringTools):
MV := parse(Trim(StringSplit(Split(interface(version), ",")[2], "Maple")[2])):

P := NULL:
for i to 22 do
  Threads:-Sleep(0.2):
  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 ) )
  else
    InsertContent( Worksheet( Group( Input( Tlayout ) ) ) )
  end if:
end do:

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

@vs140580 
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 https://en.wikipedia.org/wiki/Binomial_distribution 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 https://mathoverflow.net/questions/323845/can-the-sum-of-identically-distributed-dependent-bernoulli-trials-be-binomially (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));
use
[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:

  • https://projecteuclid.org/journals/annals-of-probability/volume-3/issue-3/Poisson-Approximation-for-Dependent-Trials/10.1214/aop/1176996359.full
  • https://www.d.umn.edu/~yqi/mydownload/12letters.pdf

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

@acer 

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

@acer 

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