## dharr

Dr. David Harrington

## 4791 Reputation

19 years, 34 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## randomize...

@radaar My code was not Metropolis and is not much use without the whole project. But basically I used

Statistics:-Sample(Statistics:-RandomVariable(Normal(0,sdev)),nfreqs) in a loop (I wanted Normal distribution)

Note that if you rerun a worksheet you get the same numbers, but you can use randomize() to get different numbers every time you run the worksheet. The actual sequence is very long.

## zero determinant...

@Magma You say "and also the determinant of A should be non-zero" but I'm not clear whether than means it follows automatically or it has to be another test.

I would have thought if you have an irreducible polynomial of degree n which doesn't have zero as a root (such as the one you gave), then the determinant is automatically going to be nonzero? So you only have to test if the characteristic poiynomial is equal to the required polynomial, and then you are done (or could do the determinant test only when the polynomials match to improve efficiency).

## Monte Carlo...

@radaar I recently used the Statistics method for a Monte Carlo method of estimating parameter errors in nonlinear regression and it seemed to work well.

But I don't understand - your timing tests suggested it was faster than rand() and now you are saying rand is faster? Any of them can be used for Metropolis, surely.

## not random?...

Why don't you think the Statistics package method is random? This is what I would use.

## log(L) simplification...

Nice worksheet!

To simplify the log of the product to a sum you need to tell Maple N::posint (not just integer). The following works:

LL := expand(log(L)) assuming positive;

I'm assuming in the first part you didn't want to use the Statistics package, but to avoid typing errors for the PDF, I might have done

Statistics:-PDF(Normal(mu,sigma),x);

## like totally positive matrices...

@Magma A final suggestion might be to look into the literature on totally positive matrices (all minors nonnegative) or strictly totally positive matrices (all minors positive), where there are improvements on how many determinants you need to check - for example you only need check n^2 minors for totally positive matrices, not all of them. These originate in relationships between various determinants, which might help for the MDS case though it seems the MDS case will likely be more difficult. On the other hand, all strictly totally positive matrices are MDS, so if you just want to find some MDS matrices, this might be helpful.

A good starting point is J. Pen~a, SeMA Journal, 62 (2013) 61.

## Determinant methods...

@Magma My point was mainly that doing the inverse is slower than the adjoint, and you can keep everything in integer arithmetic. And you're right that mul is probably not the best way to check all zeroes. But thinking about it, adjoint has to do the same determinants anyway, so I think the inverse idea really is not any better. I tried to do it recursively (often efficient in Maple), but it took forever.

But note that Determinant has the possibility of using many different algorithms, so you could optimize the routine to some extent. And see ?LinearAlgebra,Generic,Determinant for working in different fields.

@Magma If Determinant(A) is 0, then it's not MDS and you are done, but otherwise just calculate the Adjoint rather than the inverse. Then for an integer matrix the Adjoint is another integer matrix. Multiplying all the entries and checking if the answer is zero is also integer-only operations. So for 16x16 this is certainly feasible. Surprisingly, mul(Adjoint(A)) is much faster than working out the Adjoint and then using mul in a seperate step [edit: no, the display just slows things down], and takes about 0.02 s on my laptop.

Then you can work down to the smaller matrices.

## oops...

@Thomas Richard that's embarrasing; its years since I made that mistake in my own worksheets. Perhaps pi should be prettyprinted in grey and Pi in blue, since it is a common error.

## shouldn't but does...

@vv taylor can do this, but I agree it shouldn't. I guess it just passes it to asympt.

 > restart;kernelopts(version);

 > taylor( (2 * pi * n) ^ (1/(2*n)), n = infinity, 2);

 > asympt( (2 * pi * n) ^ (1/(2*n)), n, 2);

 >

## with plot...

@eslamelidy plot can be done as on the attached (adjust imax to suit).

final4.mw

## table and integrals...

@eslamelidy If you use a variable b[2], then this is part of a table called b, so you cannot then use b as a variable name. I fixed this and showed how to get the table outputs you want. I'm assuming you want all these quantities at z=500.

However, you have not properly converted the integrals. The integrand of M__2 had only the integration variable t (and not z) so I used eval to substitute t for z to make the differential equation (this was just to avoid making editing errors and save me time, you could just enter the z's by hand). I then had a differential equation for M__2. But then this gets multiplied by z in creating M, which is a combination of M__2 etc. But there should be no differential equation for M. You have somehow combined into a single differential equation for M, but this has z and t together, and when you substitute t for z you just get zero, which you can see in the output of dsol3: M(z) and Mc(z) are always zero.

So you only convert M__2 etc that were integrals to differential equations, not M or Mc or N or Nc.

## yes...

@eslamelidy Yes, in my method the integrals (M[1] etc) are calculated at the same time as the (other) differential equations, which makes it efficient. (You can compare these with the integrals you had after if you like.)

## Partially done...

@eslamelidy I've started to change the worksheet, and think that if you do the rest it should work, Just to repeat, you don't need the integrals if you have the corresponding ODEs, for example f(z) = Int(sin(t),t=0..z) is the solution to diff(f(z),z)=sin(z) with initial condition f(0)=0.

final2.mw

## Simpler formulation...

@eslamelidy You have done a lot of cutting and pasting to get the "known" functions H(z) etc, which presumably are related simply to M[1] etc, but this is not clear in the worksheet you have given, so it is not obvious to me how to implement the details of the idea I am suggesting. You need to clearly specify how H(z) etc were constructed from the simpler functions.