epostma

1394 Reputation

17 Badges

13 years, 134 days
Maplesoft

Social Networks and Content at Maplesoft.com

Maple Application Center
I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are answers submitted by epostma

Hi GaoCG,

You'll need to provide initial values for your parameters so that Maple has a starting point for iterating towards a locally optimal solution. You can do this with the initialvalues option explained on the Statistics/Fit help page. Also, there seems to be an extra parenthesis in the call to Fit; I assume that's the one just before ', X, Y, sigma', right?

It's not so easy to find values of the parameters that lead to numbers that are anywhere near reasonable. Maybe you know the actual value of the parameters for some other instance of the model? Since there's a factor of sigma^2-(650*10^6)^2 in the denominator, where sigma^2 is on the order of a few hundred, this is a big negative number; since we need to find values that are positive, we need to compensate for that. One option would be to set A to a big negative number.

Another thing to remember is that there might be bad cancellation going on: because you're subtracting such a big number from sigma^2, if you want to see the effect of sigma on the result, you'll need sufficiently many digits to represent it accurately. It might be useful to increase Digits to a higher number, like 25.

Hope this helps,

Erik Postma
Maplesoft.

Hi hossayni,

What does currentdir(); print in these two worksheet? Maple will search for .m files in that directory, and if it doesn't match you either need to specify the path more fully, or change the current directory. See ?currentdir.

Hope this helps,

Erik Postma
Maplesoft.

Hi strawfields93,

The issue is in your use of 2n - if you use 2*n (in all three places where it occurs) it should work. (It did for me, at least.)

HTH,

Erik Postma
Maplesoft.

Hi upperhilsdale,

Wondering what your input is given as, before you turn it into splines - is it just an unordered set of data points, or of weighted data points? If it's either of those, then the easiest solution is to use the Statistics package.

# Generate x-values and weights:
X1 := Statistics:-Sample(Normal(5, 2), 10^4):
W1 := Statistics:-Sample(ChiSquare(3), 10^4):
X := cos~(X1) + X1:
W := abs~(sin~(X1 + W1) + W1):

Statistics:-PointPlot(W, 'xcoords' = X);
Statistics:-Median(X);
Statistics:-Median(X, 'weights' = W);

Median gives you the 50% quantile; if you'd want to explore other values you can use the Quantile command.

If your original data truly represents 'points on a curve', then it's almost certainly best to do what acer suggests, but just for fun you could possibly generate random data points according to the curve and then take the median of those points - however, it's quite tricky to get Statistics:-Sample to work with a piecewise-defined PDF (I didn't manage to make it work in the short time that I tried).

HTH,

Erik Postma
Maplesoft.

Hi PrimZor,

You ask how to transform the result into a "usable" result. What would you like to do with your results? For example if you want to find out the values of x3-1 at the solution(s) of your system of equations, you can write:

a1 := solve({x^2=1, x>0}, x);
eval(x^3-1, a1);

or

a2 := solve({x^2=1, x>-5}, x);
eval~(x^3-1, [a2]);

Unfortunately, you'll need to know beforehand whether there will be 0, 1, or more solutions with this syntax. It's a bit more predictable if you enclose the variable(s) you want to solve for in list brackets:

solve({x^2=1, x>0}, [x]);

returns a list of solutions (just the one in this case), each being a list of equations (also just one in this case) defining the solution. So this way you can always do

a3 := solve({x^2=1, x>0}, [x]);
eval~(x^3-1, a3);
a4 := solve({x^2=1, x>-5}, [x]);
eval~(x^3-1, a4);

and get a list of the values at all the solutions of your equation.

Hope this helps,

Erik Postma
Maplesoft.

Hi Mac Dude,

The answer really is: it depends. In general using Vectors and Matrices is a safe choice, so I'd stay with what you're doing now if you don't have a pressing reason to switch. I would say that for most applications Vectors and Matrices are probably the best choices. If you would find particular performance problems with these two, you can always reconsider.

By the way, lists are not a superset of Vectors; they're different animals altogether. All of the types you've posted are distinct. (There is a common term for Vectors and Matrices (and Arrays), though, which is rtable, for rectangular table.)

Hope this helps,

Erik Postma
Maplesoft

On the two machines where I tested it, this seems to be about 10% faster:

listhist2 := proc(L)
    [seq([i, numboccur(i, L)], i = {op(L)})];
end proc;

Erik.

Hi mitko,

I don't know if you're interested in getting the execution of this single line of code sped up; it's unlikely to help you much at this point, but I saw an opportunity for improvement and figured I'd let you, and maybe more importantly the rest of the forum, know.

There are in fact two options for speedup that I see here. The first is to use ArrayTools:-Alias instead of convert(..., Vector). In particular, if your Matrix is called M, then the following alternatives will - often - give you the same result, or at least a result that looks the same:

v1 := convert(M, 'Vector');
v2 := ArrayTools:-Alias(M, [numelems(M)]);

There are two differences between v1 and v2. The first is the cause of the word 'often' above. Maple can store Matrix data in a number of ways: influenced by its indexing function (a.k.a. shape; see ?rtable_indexfcn), and its storage and order options (see ?Matrix and ?rtable_options). The default is no indexing function, rectangular storage, and Fortran order; if that's the case, then v1 and v2 will be the same Vector.

The other difference is the cause for ArrayTools:-Alias to be faster - especially if your initial Matrix is large. convert(..., Vector) creates a new Vector and copies the Matrix data into it. ArrayTools:-Alias, on the other hand, simply creates a new reference to the same block of memory, treating it as a Vector now. (Read the warnings on the package: this can go horribly wrong if you strip the attributes off the resulting Vector, and let the original Matrix be destroyed while your new reference is still around.) If you are constructing temporary objects, to be used for the construction of a different object, as you are doing here, then ArrayTools:-Alias is an excellent choice.

This causes a speedup of a factor of 20 for a 3x3 Matrix to almost 6000 for a 300x300 Matrix. (I like to test this sort of stuff using the new features of the iterations option of CodeTools:-Usage in Maple 16, where you can say 'iterations = 2*Unit(s)' to get it to choose the number of iterations so that it runs for about 2 seconds.)

The other option for speedup is using the kernel function `<|>`, which you normally wouldn't call explicitly but which is the functionality behind the Matrix/Vector shortcut syntax <a | b | c>. (It does this together with its sibling, `<,>`.) In particular, if vs is a sequence of column Vectors, then

M1 := ArrayTools:-Concatenate(2, vs);
M2 := `<|>`(vs);

will return the same Matrices - but in the test I did (with 10 Vectors of 10 elements each), `<|>` was about a hundred times faster. This is because `<|>` is pure kernel code whereas ArrayTools:-Concatenate has a substantial amount of library code.

Erik Postma
Maplesoft.

Hi mah00,

The case where sigma is degenerate (the case in which you are) is a special case that you'll need to handle specially. The PDF for this case is not well-defined. Let that sink in: there is no PDF for alpha - it's impossible. Wikipedia has a nice paragraph about it: http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case.

So you need to consider what the alternatives are. You don't write why you want to know the PDF, but I think what I would recommend in most cases is to find an n x 12 transformation matrix A and the vector of means m such that alpha = A . X + m, where X is an n-dimensional vector of independent standard Gaussians (i.e. mean 0 and standard deviation 1) and A has rank n. You can then do whichever analysis you need to do using X.

You can find this as follows. First, find the vector of means: this is simply the vector containing the ExpectedValue of each alpha_i. Subtract this from alpha to obtain alpha'. Now alpha' has mean 0 in every entry.

Now compute your matrix of covariances with respect to alpha'. Let's call this matrix sigma'. Like sigma, it will be a 12 x 12 symmetric matrix of incomplete rank. That means that it has some singular values that are zero. Compute these with

u, s := LinearAlgebra:-SingularValues(sigma_prime, 'output = [U, S]');

Check which entries of s are zero; for example, maybe this is nrs 11 and 12. That means the actual information is in entries 1 .. 10. We have now found that n = 10. Now set

A := u[.., 1 .. 10] . Matrix(s[1 .. 10], 'shape = diagonal') . u[1 .. 10]^+;

This will be your matrix A. You have alpha = A . X + m; use that instead of alpha.

Let us know if you need any more help.

Erik Postma
Maplesoft.

Hi Kitonum and Xilyte,

I saw that Kitonum found a great solution, but you can see that Maple gets stuck in a situation where the local search methods do not find a great solution. I played with it for a bit and adapted it to the following solution, where we first fit the one exponent coefficient (assuming that the other is twice as small) and then, fixing one, fit the second exponent:

restart: with(plots): with(Statistics):

X:=<0, 7.17, 18.11, 34.34, 57.95, 91.84, 139.98, 207.94, 303.48, 437.57, 625.87, 890.96, 1265.91, 1800>:
Y:=<.44, .43, .42, .41, .40, .39, .38, .37, .36, .35, .34, .33, .32, .31>:

expr := a+b1*exp(-t/c1)+b2*exp(-t/c2);

fit:=Fit(eval(expr, c2 = 2*c1), X, Y, t, parameterranges=[a=0..1, c1=1..100000], output=[parametervalues, leastsquaresfunction]);
fitplot1 := plot(fit[2], t=0..1800, thickness=2):
graph1 := pointplot(X, Y, symbol = BOX, symbolsize=15, axes = BOXED):
display(graph1, fitplot1);

fit2:=Fit(eval(expr, c1 = 191.),X, Y, t, parameterranges=[a=0..1, c2=1..100000]);
fitplot2 := plot(fit2, t=0..1800, thickness=2, color=red):
graph1 := pointplot(X, Y, symbol = BOX, symbolsize=15, axes = BOXED):
fitplot1 := display(fitplot1, transparency=0.5):
display(graph1, fitplot1, fitplot2);

Hope this helps,

Erik Postma
Maplesoft.

In the Appendix at the end of the Programming Guide (available on the web at http://www.maplesoft.com/support/help/Maple/view.aspx?path=ProgrammingGuide/Appendix) there is a lot of information about tables. In particular, you can see that a Maple-language table is a TABLE dag, the contents of which are in a HASHTAB dag which consists of HASH dags. So these guys are hash tables. That means they're pretty efficient - the relevant operations are amortized O(1) in principle.

If you have specific questions about the implementation, let us know.

Hope this helps,

Erik Postma
Maplesoft.

I think you'd have to write this yourself, but it wouldn't be too much work. The equivalent of the ratweight command would simply write the weights into a table (global, or even better, (optionally) user-supplied) and you'd write a command called something like evaluateWithRatWtLvl that takes an expression and an integer (and optionally the weight table, or otherwise using the global weight table), expands the expression, and throws away terms with too high weight.

If you'd like, you could write a package that overloads certain operators, like + and ^, to do an automatic call to evaluateWithRatWtLvl after performing their computation. This will only apply to code you type at the top level, though, not to internal computations done in Maple library code.

Hope this helps,

Erik Postma
Maplesoft.

The parameters in your expression to be fitted, Amplitude, seem to be x, t, n, and S (after, as Axel pointed out, you substitute I for i). However, in the call to NonlinearFit, you give parameter ranges for a wholly different set of parameters - De, S, and tau. That doesn't work - you should either rewrite Amplitude to depend (only) on De, S, and tau, or supply ranges for (only) x, t, n, and S.

Hope this helps,

Erik Postma
Maplesoft.

Hi alix,

This type of thing is usually easier in Maple with expressions than with functions. If you have an expression

expr := x + 3;

then you can get the operands x and 3 using ?op:

op(expr)

will return x, 3. The function combining these to get expr back can be obtained with op(0, expr) ; it will return `+` in this case.

Things get a little complicated if there are data structures such as ?tables or ?Vectors involved, or indeed procedures such as in your example.

Hope this helps,

Erik Postma
Maplesoft.

You can't really define w as both a bivariate function w(ξ,s) and a univariate function w(ξ). I'd suggest using notation like

w2 := (ξ,s) -> N(s) . w1(ξ);

where w2 is the bivariate and w1 the univariate function. Note that I wrote the function definition in the non-ambiguous arrow notation; writing it as 

w2(ξ,s) := N(s) . w1(ξ);

may work, depending on your settings, but the former always does. Note also the use of the period (.) for matrix and vector multiplication.

If you show us a little more of what you're trying to achieve, I bet we'll find out that it's even better to not define w2 itself, but just a definition equation

w2def := w2 = (ξ,s) -> N(s) . w1(ξ);

and then use eval(something, w2def) in some well-chosen spots. This way, you can control when Maple applies this definition of w2; otherwise, it is always applied immediately when Maple encounters the function w2 and you will never get a result that involves w2.

As to your other question: this is not always easy.  We would need to see a bit more of the problem you're trying to solve in order to answer it. By the way, it may or may not work to use the same symbol s for the index of x and its argument; if you substitute a value like 2 or 5 for s, you would get x2(2) and x5(5), not xs(2) or xs(5). It might be better to write this as xs(s) - Maple will never break up a single "word" like that. (There are other solutions involving things called "inert identifiers" - but I'd advise to do that only if it's really worth the added complexity.)

Hope this helps,

Erik Postma
Maplesoft.

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