acer

32490 Reputation

29 Badges

20 years, 7 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@mmcdara In my Maple 2019.1 the computation of Variance(Z(X+Y)) completes, slowly, under these assumptions on sigma__x and sigma__y.

It also sometimes completes without assumptions, if interrupted and re-computed. That's a sign that int is struggling. An important difficulty is that in nested/iterated form int(int(int(...))) the inner integrals are not aware that the outer variables of integration may be real, or lie within a specific range, etc.

But using the so-called collapsed form int(...,[...]) it computes reasonably quickly under those same assumptions.

I will submit another Software Change Request, to the effect that Statistics should probably use that collapsed form for all its multiple intergals -- whether active or inert.

The same goes for VectorCalculus:-int and Student:-MultivariateCalculus:-MultiInt. (I have a growing collection of simple examples where it makes the difference between success versus failure, or a quick versus a slow computation.)

So, for now, value(CollapseNested(Combine( some_inert_nested_integral ))) can be a useful alternative.

restart:

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

with(Statistics): with(IntegrationTools):

X := RandomVariable(Normal(mu__x, sigma__x)):
Y := RandomVariable(Normal(mu__y, sigma__y)):

U := X+Y:
Z := q -> cos(q):
v__z := Variance(Z(X)):
V := subs({mu__x=Mean(U), sigma__x=StandardDeviation(U)}, v__z);

-(1/2)*(2*cos(mu__x+mu__y)^2*exp(sigma__x^2+sigma__y^2)-2*cos(mu__x+mu__y)^2-exp(2*sigma__x^2+2*sigma__y^2)+1)*exp(-2*sigma__x^2-2*sigma__y^2)

bar := CollapseNested(Combine(Variance(Z(X+Y),inert))):

raw1 := CodeTools:-Usage( value(bar) ) assuming sigma__x>0, sigma__y>0;

memory used=349.05MiB, alloc change=76.00MiB, cpu time=3.69s, real time=3.41s, gc time=478.88ms

(1/4)*(-exp(sigma__x^2+sigma__y^2)-2*exp((2*I)*mu__x+sigma__x^2+sigma__y^2+(2*I)*mu__y)+2*exp(2*sigma__x^2+(2*I)*mu__x+(2*I)*mu__y+2*sigma__y^2)+exp((4*I)*(mu__x+mu__y))-exp((4*I)*mu__x+(4*I)*mu__y+sigma__x^2+sigma__y^2)+1)*exp(-2*sigma__x^2-(2*I)*mu__x-(2*I)*mu__y-2*sigma__y^2)

S1 := simplify(evalc(simplify(raw1)));
is(S1=V);

(1/2)*(-cos(2*mu__x+2*mu__y)-1)*exp(-sigma__x^2-sigma__y^2)+(1/2)*exp(-2*sigma__x^2-2*sigma__y^2)*cos(2*mu__x+2*mu__y)+1/2

true

raw3 := CodeTools:-Usage( Variance(Z(X+Y)) ) assuming sigma__x>0, sigma__y>0;

memory used=5.68GiB, alloc change=345.59MiB, cpu time=55.16s, real time=50.86s, gc time=7.13s

-(1/4)*(exp(sigma__x^2+sigma__y^2)+2*exp((2*I)*mu__x+sigma__x^2+sigma__y^2+(2*I)*mu__y)-2*exp(2*sigma__x^2+(2*I)*mu__x+(2*I)*mu__y+2*sigma__y^2)-exp((4*I)*(mu__x+mu__y))+exp((4*I)*mu__x+(4*I)*mu__y+sigma__x^2+sigma__y^2)-1)*exp(-2*sigma__x^2-(2*I)*mu__x-(2*I)*mu__y-2*sigma__y^2)

S3 := simplify(evalc(simplify(raw3)));
is(S3=V);

(1/2)*(-cos(2*mu__x+2*mu__y)-1)*exp(-sigma__x^2-sigma__y^2)+(1/2)*exp(-2*sigma__x^2-2*sigma__y^2)*cos(2*mu__x+2*mu__y)+1/2

true

 

Download stats_example.mw

@qilianshan In your input you can type Ctrl-L to insert such an equation label reference.

@wsningrat I was striving to provide you with code whose approach was straightforward and where the plotting methodology seemed natural.

As it happens it was already reasonably quick. But for your interest here it is, but faster, and (I hope) still reasonably straightforward as far as the plotting goes. I also include Tom's approach, which is also pretty fast, and you could easily go with either (I leave it to you to decide which suits your style).   captive_breeding_ac2.mw

@chandrashekhar 

Your claims appear to be untrue. Please point out what is incorrect, specifically, with the worksheet I now attach.

Perhaps you should cut and paste what was in your own post (rather than work with any prior worksheet), and test your claims.

Please check that the expressions appearing below are the same as what you wrote in your original Question.

restart;

#
# I cut and pasted the original expression.
#

firstexpression := 4*beta*c*lambda^2*mu^2+2*beta*c*lambda*mu^3
 +2*beta*c*lambda*mu^2*nu+beta*c*mu^3*nu-4*c*lambda^2*mu^2*sigma
 -2*c*lambda*mu^3*sigma-2*c*lambda*mu^2*nu*sigma-c*mu^3*nu*sigma
 +4*beta*lambda^3*sigma+4*beta*lambda^2*mu*sigma+2*beta*lambda^2*nu*sigma
 +2*beta*lambda*mu^2*sigma+2*beta*lambda*mu*nu*sigma+beta*mu^3*sigma
 +beta*mu^2*nu*sigma+4*lambda^2*mu^2*sigma+2*lambda*mu^3*sigma
 +2*lambda*mu^2*nu*sigma+mu^3*nu*sigma;

4*beta*c*lambda^2*mu^2+2*beta*c*lambda*mu^3+2*beta*c*lambda*mu^2*nu+beta*c*mu^3*nu-4*c*lambda^2*mu^2*sigma-2*c*lambda*mu^3*sigma-2*c*lambda*mu^2*nu*sigma-c*mu^3*nu*sigma+4*beta*lambda^3*sigma+4*beta*lambda^2*mu*sigma+2*beta*lambda^2*nu*sigma+2*beta*lambda*mu^2*sigma+2*beta*lambda*mu*nu*sigma+beta*mu^3*sigma+beta*mu^2*nu*sigma+4*lambda^2*mu^2*sigma+2*lambda*mu^3*sigma+2*lambda*mu^2*nu*sigma+mu^3*nu*sigma

#
# I cut and pasted the expression which
# chandrashekhar claims was simplified by hand.
#

secondexpression := beta*mu^3+(2*lambda+nu)*(beta*mu^2
 +2*beta*lambda*mu+2*beta*lambda^2+(1-c)*(2*lambda+mu)*mu^2);

beta*mu^3+(2*lambda+nu)*(beta*mu^2+2*beta*lambda*mu+2*beta*lambda^2+(1-c)*(2*lambda+mu)*mu^2)

#
# chandrashekhar claims this will produce 1.
#
# It does not do so.
#

simplify(firstexpression/secondexpression,size);

(((-2*c*beta+2*sigma*(c-1))*lambda+(-c*nu-sigma)*beta+nu*sigma*(c-1))*mu^3-2*(2*lambda+nu)*((c*beta-sigma*(c-1))*lambda+(1/2)*sigma*beta)*mu^2-2*beta*sigma*lambda*(2*lambda+nu)*mu-2*beta*sigma*lambda^2*(2*lambda+nu))/(((2*c-2)*lambda-beta+(c-1)*nu)*mu^3+2*(2*lambda+nu)*((c-1)*lambda-(1/2)*beta)*mu^2-2*beta*lambda*(2*lambda+nu)*mu-2*beta*lambda^2*(2*lambda+nu))

#
# chandrashekhar claims this will produce 0.
#
# It does not do so.
#

simplify(firstexpression-secondexpression,size);

((2*c*beta-2*(sigma-1)*(c-1))*lambda+(c*nu+sigma-1)*beta-nu*(sigma-1)*(c-1))*mu^3+2*((c*beta-(sigma-1)*(c-1))*lambda+(1/2)*beta*(sigma-1))*(2*lambda+nu)*mu^2+2*beta*lambda*(sigma-1)*(2*lambda+nu)*mu+2*beta*lambda^2*(sigma-1)*(2*lambda+nu)

 

Download claims2.mw

I have already asked you to upload a worksheet that contained the expressions as you expected them. You post comment after comment without any such worksheet. It's not clear that you are really interested in help us to help you.

For what it might be worth, here are three forms of the first expression (as I copied it) that are compact. 

(mu^2*(beta-sigma)*(mu+2*lambda)*c+(2*beta*lambda^2
+2*mu*(mu+beta)*lambda+mu^2*beta)*sigma)*(nu+2*lambda)
+sigma*(nu+beta+2*lambda)*mu^3;

(mu^2*(beta-sigma)*(mu+2*lambda)*c+2*(beta*(lambda^2+1/2*mu^2)
+mu*(mu+beta)*lambda)*sigma)*(nu+2*lambda)
+sigma*(nu+beta+2*lambda)*mu^3;

mu^2*(beta-sigma)*(nu+2*lambda)*(mu+2*lambda)*c
+(nu+2*lambda)*(2*beta*lambda^2+2*mu*(mu+beta)*lambda
+mu^2*beta)*sigma+sigma*(nu+beta+2*lambda)*mu^3;

You should specify what metric you want to use to measure size of the expression. Maple's length command is not best. The commands MmaTranslator:-Mma:-LeafCount or `codegen/cost` may be better.

Does anyone else get that the OP's original and his manually simplified expression are actually equivalent?

And Christian's reformulation?

Perhaps I'm missing something in my copy&paste.

I find it difficult to take this seriously if you deliberately avoid providing a complete example that reproduces all the problem's characteristic issues.

@mmcdara Why do you consider that non-numeric indexing won't work?

And inside a procedure, you only need to declare the base-name theta as local.

restart:
vars := seq(theta[i], i=1..4);
Array(1..4, i -> diff(w(vars), theta[i]));

vars := theta[1], theta[2], theta[3], theta[4]

Array(%id = 18446884382838275238)

restart:
MyIndices := [a, b, c, d];
vars := seq(theta[n], n in MyIndices);
Array(1..4, i -> diff(w(vars), vars[i]));
 

MyIndices := [a, b, c, d]

vars := theta[a], theta[b], theta[c], theta[d]

Array(%id = 18446884382838268022)

restart;

W := w(theta[a],theta[b],theta[c],theta[d]);

w(theta[a], theta[b], theta[c], theta[d])

diff(W, theta[a]);

diff(w(theta[a], theta[b], theta[c], theta[d]), theta[a])

seq(diff(W,theta[v]), v=[a,b,c,d]);

diff(w(theta[a], theta[b], theta[c], theta[d]), theta[a]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[b]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[c]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[d])

vars := seq(theta[n], n in [a, b, c, d]);
seq(diff(W,r), r=vars);

theta[a], theta[b], theta[c], theta[d]

diff(w(theta[a], theta[b], theta[c], theta[d]), theta[a]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[b]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[c]), diff(w(theta[a], theta[b], theta[c], theta[d]), theta[d])

 

Download diff_indexed_name.mw

 

@q41b3 The name P was already used, so some other name(s) are needed to which to assign the constructed matrices. You could also use a new base-name such as, say, Pmat.

You don't necessarily need to use concatenation to obtain the new names to which to assign the matrices. You could also using indexing on that base-name. Eg, Pmat[0]:=P and then in the loop Pmat[n] := evalm(Pmat[n-1] &* P[0])

And, no, you don't necessarily need to create/declare a structure in advance, to hold those Pmat[n] matrices. If you don't then Maple will implicitly create Pmat as a table (into which you can index). But if you prefer you could also create an Pmat:=Array(1..10) in advance.  The table structure is handy because you don't need to assign Pmat:=table([]) in advance (though you could) and its size is not fixed.

Using concatenation to construct the names isn't very elegant, in my opinion.

If you put the code into a procedure then using Pmat as a table or Array is convenient because you only have to declare the one local Pmat.

Here is the code using the more modern (captilized) Matrix instead of the deprecated (lowercase) matrix. You can notice that evalm is not required for arithmetic, and it prints directly. And matrix-matrix multiplcation is now done using the noncommutative dot (ie, . ).

I have also done this using 1D Maple notation (plaintext input) rather than marked-up 2D Input. I find it easier to edit, and unlike 2D Input it is WYSIWYG. It can be set as a global preference for new Worksheets under Tools->Options from the main menubar.

restart;

N := 2: A := -N: B := N:

q := 0.3: p := 0.5: sa := 0.9: sb := 0.1: r := 1 - p - q:

dimP := 2*N + 1:

P := Matrix(dimP, dimP):

P[1, 1] := sa: P[1, 2] := 1 - sa:
P[dimP, dimP] := sb: P[dimP, dimP - 1] := 1 - sb:

for i from 2 to dimP - 1 do
    P[i, i - 1] := q;
    P[i, i] := r;
    P[i, i + 1] := p;
end do:

Pmat[0] := P;

Matrix(5, 5, {(1, 1) = .9, (1, 2) = .1, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = .3, (2, 2) = .2, (2, 3) = .5, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = .3, (3, 3) = .2, (3, 4) = .5, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = .3, (4, 4) = .2, (4, 5) = .5, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = .9, (5, 5) = .1})

for n to 10 do
    Pmat[n] := Pmat[n-1] . Pmat[0];
end do:

Pmat[10];

Matrix(5, 5, {(1, 1) = .5448343619400003, (1, 2) = .11532848416000005, (1, 3) = .12176249215000005, (1, 4) = .14498696925000004, (1, 5) = 0.7308769250000002e-1, (2, 1) = .3459854524800001, (2, 2) = .10282244927000003, (2, 3) = .15926588350000004, (2, 4) = .2516929085000001, (2, 5) = .14023330625000005, (3, 1) = .21917248587000002, (3, 2) = 0.9555953010000001e-1, (3, 3) = .18078069908, (3, 4) = .3246936532, (3, 5) = .17979363175, (4, 1) = .15658592679, (4, 2) = 0.9060944705999999e-1, (4, 3) = .19481619192, (4, 4) = .35339349113, (4, 5) = .2045949431, (5, 1) = .14208247422000003, (5, 2) = 0.9087118245e-1, (5, 3) = .19417712229, (5, 4) = .36827089758000003, (5, 5) = .20459832346000004})

Pmat[0]^11;

Matrix(5, 5, {(1, 1) = .5448343619400002, (1, 2) = .11532848416000004, (1, 3) = .12176249215000003, (1, 4) = .14498696925000004, (1, 5) = 0.7308769250000002e-1, (2, 1) = .3459854524800001, (2, 2) = .10282244927000002, (2, 3) = .1592658835, (2, 4) = .25169290850000003, (2, 5) = .14023330625000002, (3, 1) = .21917248587000002, (3, 2) = 0.955595301e-1, (3, 3) = .18078069908, (3, 4) = .3246936532, (3, 5) = .17979363175000002, (4, 1) = .15658592679000002, (4, 2) = 0.9060944706e-1, (4, 3) = .19481619192000005, (4, 4) = .35339349113000007, (4, 5) = .20459494310000007, (5, 1) = .14208247422000003, (5, 2) = 0.9087118245000002e-1, (5, 3) = .19417712229000003, (5, 4) = .3682708975800001, (5, 5) = .20459832346000004})

 

Download TransitionMatix_ac_Matrix.mw

@q41b3 

When you type else if you are doing two things with two otherwise unrelated key words:
1) You are starting an else clause in the current if..end if (of which there can only be one)
2) You are then starting a new, deeper level if..end if .

Your original went wrong because it started several nested/deeper if statements that lacked their corresponding end if bits. But even with those included there would still be another problem of trying to put multiple else clauses on the outermost if..end if. See the last attempted example below.

Perhaps the following will make it clear. The use of white-space and indentation is there to try and make it easier to understand.

restart;

A := proc(x)
  if x>4 then
    4;
  elif x>3 then
    3;
  elif x>2 then
    2;
  elif x>1 then
    1;
  end if;
end proc:

A(4.5), A(3.5), A(2.5), A(1.5);

4, 3, 2, 1

# This is harder to read and manage. Each "if" requires
# its own "end if". And the nesting is done because
# any "if..end if" can have only a single "else" clause.
#
B := proc(x)
  if x>4 then
    4;
  else if x>3 then
         3;
       else if x>2 then
              2;
            else if x>1 then
                   1;
                 end if;
            end if;
       end if;
  end if;
end proc:

B(4.5), B(3.5), B(2.5), B(1.5);

4, 3, 2, 1

# This is not valid syntax, since any "if..end if" can
# have only a single "else" clause. The error message
# is indicating that the appearance of the second "else"
# is unexpected.
#
proc(x)
  if x>4 then
    4;
  else if x>3 then
         3;
       end if;
  else if x>2 then
         2;
       end if;
  else if x>1 then
         1;
       end if;
  end if;
end proc:

Error, reserved word `else` unexpected

Download elif.mw

@tomleslie Isn't that what my answer already did above (except for the highly trivial aspect of making T a parameter of an additional layer of encapsulating procedure)?

It looks like what the OP may be after, judging by his picture.

Given the disparity between the OP's original odeplot call and his posted picture, I'd suggest that understandability trumps maximal computational efficiency on a first run. (Things like maximal efficiency, adaptive plotting vs fixed mesh, ease of specifying which curves are included, etc, are additional considerations.)

At present these is no direct option of the Tabulate command to control the font size for either text or math.

If I recall correctly the font used for numeric entries (data) should match the worksheet's "2D Input" style, while the row and column headers should match the worksheet's "Text" style. You can see the main menubar's Format->Styles to check your worksheet's defaults for those.

When you simply print the DataFrame in the usual way then what you see is displayed in the "2D Output" style.

Can I ask, which look too big for you: the row/column headers or the data entries? What Operating System and Maple version are you using?

Are you interested in a solution which would allow you to customize the Tabulate entries individually cell-by-cell, or would you be ok with a solution that only let you put a blanket adjustment over all data entries and another for all row/column headers?

Would you be interested in a procedure that get's what you need, but only handles DataFrames and not other structures, including those with math and plots in them? (It's slightly easier to write one for your use now, then it is for me to come up with an edit to Tabulate itself). Or perhaps you wouldn't be interested in any customized alternative right now, because it's be awkward to re-use it?

 

@mugwort Perhaps you might consider uploading your actual problems, in worksheets, when you first post future Questions.

@Carl Love Yes, a dedicated tridiagonal solver can beat a general band solver, as you've nicely shown. And even more so if the band solver is not set up for convenient re-use of pre-factorization.

I have written solvers (and re-usable prefactorization based solvers) for both band and tridiagonal myself; as you've shown it is reasonably straightforward. The more specialized solver outperforms the more general, and the prefactorized usage does better still.

So, if one can quickly write a numeric solver that handles a specific Matrix shape better, why doesn't Maple cover it? That is a sensible and interesting question, I think. I think that -- leaving aside wider issues of priorities -- the answer can include this: Writing the solver can sometines be the least involved part of the software development, with coordinated integration into the product requiring the larger part of time resources.

@Rouben Rostamian  

[edit] I realize that your worksheet may be primarily expository, and it is very nice btw. I add these comments below merely as followup to Carl's remarks. I would not be surprised is overhead of creation of Matrix A, etc, may dominate at this stage of the work. And that's fine.[/edit]

Maple has a band (float) solver that is fast here.  At sizes m=n=50 you may not be able to distinguish the performance differences. But as the problem size increases the band solver's performance superiority comes through.

At large enough Matrix size n the scheme of precomputing the full rectangular LU decomposition leads over precomputing the inverse. Worst of all, naturally, is doing the full rectangular LU+back/forward solving for each rhs V.

(Maple doesn't offer precomputing the P,LU for a band storage, in part because it needs a structure larger than input AT and things get a bit awkward functionality-wise. And it's a performance bug that MatrixInverse doesn't utilize the band solver.)

But there is an additional risk to multiplying by the inverse, in that it may not be as numerically robust.

In summary, calling LinearSolve with the float[8] band storage Matrix A is likely best for performance and robustness (unless you write your own dedicated, optimized solver for this tridiagonal case, of course).

You can utilize infolevel[LinearAlgebra] to see what's being called externally. I didn't look at the scheme for performance considerations/optimizations (inplace stuff, etc).

You can play with the N=n and maxiter=m in this attachment. (The NAG function names in the userinfo message are actually their LAPACK equivalents in the Intel MKL these days.)

band_LU.mw

@tomleslie It was not possible for many years, and it was introduced to 1D before 2D, and Section 6.9 of the Programming Guide still says the body must be a single expression.

First 216 217 218 219 220 221 222 Last Page 218 of 594