An easy but very limited solution: ?plottools/transform. It's not entirely turn-key: you have to specify the projection explicitly, independently of the orientation of the 3d plot. Say you want to project onto the plane spanned by <1,0,1> and <0,1,-2>. Then you first determine the images of the unit vectors: 

project := proc(u,v)
local uu, vv;
  uu := u / LinearAlgebra:-Norm(u, 2);
  vv := v - (v . u) * u;
  vv := vv / LinearAlgebra:-Norm(vv, 2);
  return w -> [w . uu, w . vv];
end proc:
pp := project(<1,0,1>, <0, 1, -2>);
             w -> [w . uu, w . vv]
pxyz := pp(<x, y, z>) assuming real;
             [1    (1/2)   1    (1/2)  2    (1/2)   1    (1/2)]
             [- x 2      + - z 2     , - x 5      + - y 5     ]
             [2            2           5            5         ]

Now you can use this to create a transform. Unfortunately, the 'shading = zhue' option has no equivalent in 2D plots; and worse, the corresponding element is not even taken out of the plot structure. We can do that ourselves with eval, though.

pt  := plottools:-transform(unapply(pxyz, [x, y, z])):
eval(pt(p), SHADING=(()->()));

Here p is the plot you defined above. Alas, since the shading is gone, you really can't see anymore what's going on...

Hope this inspires a better solution,

Erik Postma
Maplesoft. 

Hi Remus,

When I tried to compile simple.c just now to investigate, I initially had the same symptom as you did, but setting LD_LIBRARY_PATH to the $MAPLE/bin.$PLATFORM directory helped:

MAPLE=/home/epostma/maple15
export LD_LIBRARY_PATH=$MAPLE/bin.X86_64_LINUX

Note you need 'export' on the LD_LIBRARY_PATH definition (at least for sh-style shells).

Hope this helps,

Erik Postma
Maplesoft. 

Also off the cuff: I'm missing:

• the two nuggets in ?ArrayTools that are worth their weight in gold if you're after efficiency and reducing allocation: ?ArrayTools/Alias  and  ?ArrayTools/Copy (and to a lesser extent,  ?ArrayTools/BlockCopy);
• ?rtable_scanblock, a mini-programming language in itself;
• ?dsolve/numeric, in particular events, the 'known' option, and a related useful and little known fact about  ?evalhf described in the last bullet point of  ?evalhf/procedure ;
• ?Compiler/Compile .

You specify n -> infinity and p -> 0 and n * p -> lambda as valid on the same level, as a constraint on a process as it were - let n and p be such that eventually, n becomes arbitrarily large, p becomes arbitrarily close to zero, and n * p becomes arbitrarily close to lambda. So you specify a particular set of ways of approaching the point (infinity, 0) in (n, p)-space. I'm not sure if Maple can deal with such a general specification of a limiting process. However, if such a limit exists, then you can first let n * p -> lambda, keep (say) p as a function of n and lambda, and then let n go to infinity. This immediately gives you the last equation: evaluates to .

The real difficulty would be in getting this to fail in cases where it truly does depend on the family of ways of approaching this point...

Hope this helps,

Erik Postma
Maplesoft.

Hi Adam,

if you use Arrays (capital A) instead of arrays (lowercase a), then it's easier to massage the data into the right format. arrays are a relic of long ago that we keep around only because there might be code that people are using that still relies on it, but for all practical purposes, Arrays are superior. See their help pages:  ?Array and ?array(deprecated). One consequence (and actually a reason that lists and Arrays are in most cases much more intuitive than arrays) is that you'll need to quote the axis labels you supply to pointplot; otherwise Maple tries to compute ln(p) and ln(q) for you, which you don't want.

In this case, it might be even easier to just use plain lists. Just remove the call to array. I also used tilde operators to define lnx and lny. Then you end up with the following code (cleaned up one typo: use of 'p' in definition of lny):

with(plots);
q := [892., 1012., 1060., 987., 680., 739., 809., 1275., 946., 874., 720., 1096.];
p := [1.23, 1.15, 1.10, 1.20, 1.35, 1.25, 1.28, .99, 1.22, 1.25, 1.30, 1.05];
lnx := ln~(p);
lny := ln~(q);
pointplot({seq([ln(p[t]), ln(q[t])], t = 1 .. 12)}, symbol = circle, axesfont = [TIMES, 14], labels = '[x = ln(p), y = ln(q)]');

Now for the fit. It's almost always a good idea to look at the help page when you want to use a command, so that you can see what the valid ways of using the command are (the 'calling sequences'). Here's the one you'll need: ?CurveFitting/LeastSquares. One of the calling sequences uses the x- and y-data in separate lists (or Arrays or Matrices), and that's what we have, so it would be easiest to use that.

yfit := CurveFitting:-LeastSquares(lnx, lny, xfit);

This will define yfit as something like 7.2 - 1.9 * xfit. In order to obtain a plot of this, it's easy enough to use the  ?plot command. In order to combine this together with the result of the pointplot command, we need to give that plot a name. Let's do that first and then combine the plots with the plots:-display command. We can now move the plot options for the axis font and labels to the plots:-display command.

plot1 := pointplot({seq([ln(p[t]), ln(q[t])], t = 1 .. 12)}, symbol = circle):
plot2 := plot(yfit, xfit = min(lnx) .. max(lnx)):
plots:-display(plot1, plot2, axesfont = [TIMES, 14], labels = '[x = ln(p), y = ln(q)]'); 

In general, if you expect you'll ever want to display a variable name such as p or q, then you save yourself a bit of trouble if you don't ever assign to that variable - use a different variable name to hold the actual data. If you would replace p by, say, pdata in the code above, you wouldn't need the awkward quotes in the plot options - you could just use ln(p) without Maple thinking you somehow want the logarithm of a data structure. The same holds for q, of course.

Hope this helps,

Erik Postma
Maplesoft. 

The easiest solution is to use Statistics:-Sample instead. You get a Vector by default, but you can obtain Arrays by supplying a range or a list of ranges:

with(Statistics):
a1 := Sample(Uniform(3, 5),  25); # 25-element Vector
a2 := Sample(Uniform(3, 5), 1 .. 25); # Array(1 .. 25)
a3 := Sample(Uniform(3, 5), [1 .. 25]); # Array(1 .. 25)
a4 := Sample(Uniform(3, 5), [1 .. 6, 1 .. 4]); # Array(1 .. 6, 1 .. 4)
# If you want to do this many times for small rtables:
p := Sample(Uniform(3, 5)):
a5 := p([1..6, 1..4]); # Equivalents to the other calling sequences above are also supported.

(Minor quibble: ArrayTools:-RandomArray(..., distribution=normal) doesn't sample from the "normal distribution on the unit interval" (that's impossible) but from the standard normal distribution.)

This of course gives you access to all other distributions as well - Sample(Normal(-2, 4), 1 .. 25), Sample(StudentT(2), 10), etc.

Hope this helps,

Erik Postma
Maplesoft. 

Hi bernie,

Your best bet is probably using the parameterranges option. It's very important for Statistics[NonlinearFit] to have a good idea of the ranges of the parameters.

You don't give function and model, so I can't verify that this is enough for your case. But for example:

xvalues := Vector([seq(10^5 .. 10^6, 10^4)]):
yvalues := Vector([seq(x*ln(1.618e6 - x), x = 10^5 .. 10^6, 10^4)]):
with(Statistics):

# This will give the error you discuss, because Maple tries values of a less than 10^6:
NonlinearFit(x * ln(a - x), xvalues, yvalues, x, output = [parametervalues, residuals]);

# This will give a successful fit:
NonlinearFit(x * ln(a - x), xvalues, yvalues, x, output = [parametervalues, residuals],
    parameterranges = [a = 10^6+1 .. 10^7]);

# In general, because NonlinearFit uses Optimization which performs local optimization
# only (not global optimization), it is even better to additionally give an initial guess
# for your parameter values:
NonlinearFit(x * ln(a - x), xvalues, yvalues, x, output = [parametervalues, residuals],
    parameterranges = [a = 10^6+1 .. 10^7], initialvalues = [a = 2e6]);

Hope this helps,

Erik Postma
Maplesoft.

Hi makh,

If you're trying to perform a "normal" factorization, then you don't want to use the calling sequence Factor(polynomial) mod m . That is for factoring a polynomial modulo a prime number; e.g. Factor(x^2 + 2) mod 3 will return (x + 2)(x + 1) because the difference between x^2 + 2 and (x + 1)(x + 2) is divisible by 3 for integer x.

For regular factorization, not modulo a prime, you'll want to use the  ?factor command (with a lowercase f ); e.g. factor(x^5 - 1) returns a factorization into two factors.

Hope this helps,

Erik Postma
Maplesoft. 

Hi john2,

I had missed the interval [0, 10] that you mentioned. Yes, in such an interval it should be possible to define a kurtosis. If the function you're investigating is nonnegative on the interval (and not 0 everywhere), you can scale it to become a probability distribution (meaning you make sure that its integral over the interval is 1), and then you're home and dry. If the function is negative on some subintervals, then unless the integral of the function over the interval is 0, you can do the same and just ignore the nonnegativity requirement for probability distributions.

So let's take your example (which is nonnegative on [0, 10]). First we scale to have integral 1.

expr0 := x^3 + 2*x^2 + 3*x;
interval := 0 .. 10:
c := int(expr0, x = interval);
expr := expr0 / c;

Now we have a probability distribution. We can continue by hand or use the Statistics package. Let's first see what continuing by hand would look like (this is probably what you'd want to do if your expression is negative anywhere in the interval). We compute the mean of this distribution, then the variance, and finally the kurtosis.

mean := int(x * expr, x = interval);
variance := int((x - mean)^2 * expr, x = interval);
kurtosis := int((x - mean)^4 * expr, x = interval) / variance^2;

This corresponds to the definition you have used - there is an alternative definition of kurtosis for which you would need to subtract 3 from the value you'd get from the computation above (and from your definition in the comment).

Using the Statistics package, you'd define expr as before, and then continue as:

with(Statistics):
dist := Distribution(PDF = unapply(piecewise(x < lhs(interval), 0, x < rhs(interval), expr, 0), x));
kurtosis2 := Kurtosis(dist);

In either case you get the same result, of course (this function is somewhat leptokurtic - it has relatively fat tails). If you want a numerical approximation instead of the fraction you get as the exact result, you can either use the evalf function by typing evalf(kurtosis) or evalf(kurtosis2), or if using the Statistics package, type

kurtosis3 := Kurtosis(dist, numeric);

For all of these cases, if you have a function with a finite integral over all the reals, you can use the interval comprising all the real numbers by specifying interval := -infinity .. infinity. This would work, for example, if you started with something like

expr0 := (x^3 + 2*x^2 + 3*x) * exp(-x^2);

Hope this helps,

Erik Postma
Maplesoft.

W.l.o.g. we may assume that D <> 0, and then we can solve the equation defining the hyperspace for x4; rescaling we get

.

Thus we get a subset of the 3d cube -1 <= xi <= 1 (1 <= i <= 3) defined by

 ;

in other words, we take the subset between two arbitrary parallel planes. This reduces the question from an algebraic to a geometric one.

So we get at most two faces that are not subfaces of the original cube, plus at most six faces from the original cube, for a maximum of 8. That means the dodecahedron and icosahedron are out.

If we'd want to get the octahedron, we'd need all eight of these faces; in particular, all faces of the original cube, which are pairwise parallel and otherwise perpendicular. There are no perpendicular faces in the octahedron, so it's also out.

We can clearly obtain the cube itself by moving the two planes far away from the origin (that is, make e, a, b, c all much smaller than 1; that is, make D much larger than A, B, C, E).

This leaves the tetrahedron. I believe we can use the same argument as for the octahedron - it has only faces that make an angle of about 71 degrees with each other, so at most one of the parallel planes can lead to a face of the tetrahedron and at most one of the original faces of the cube. That leaves too few.

I'd be glad to be proven wrong in any of my arguments, because the result is a little sad: only the trivial case of a cube seems to be achievable.

Hope this helps,

Erik Postma
Maplesoft. 

It looks like neither beta1[1] nor alpha(t) are assigned a value. (Maybe you meant beta1(1) and alpha1(t)?) Also, the variable s doesn't occur in anything other than your plot command, so Maple can't really vary it. And finally, if I change beta1[1] to beta1(1) and set alpha := alpha1, I get a complex value for indrukking (by asking for evalf(indrukking)). So there are a few things to fix...

As a general guideline, I wouldn't set Digits to less than 10 if your computations aren't taking very long - if you just want to display fewer digits (instead of computing fewer digits), I'd recommend using something like interface(displayprecision = 4), which carries the full precision of the setting of Digits all around but displays only 4 digits.

Hope this helps,

Erik Postma
Maplesoft.

Hi Steweb,

I don't think I understand - what would you expect f(x) to be if not int(x^2, x = 0 .. 1) = 1/3? You define f as the function that maps an argument a to int(a * x, x = 0 .. 1), and if you substitute x for a you get the result above.

Maybe you wanted Maple to assume that a is independent of x, compute the integral, and only then substitute the argument you give? There are a few ways to do that. The most correct way would be simply to use something for a that is actually independent of x, like x1 or y. (Or a, for that matter.) Another option would be to have Maple evaluate the integral before you assign the result, as a procedure, to f; you could do that with f := unapply(int(a * x, x = 0 .. 1), a), or by letting f be not a procedure but just an expression - as f := int(a * x, x = 0 .. 1) and then using eval(f, a = something) where you used f(something) before.

Hope this helps,

Erik Postma
Maplesoft.

Hi Per,

The fitting routines in Statistics all work in the floating point domain only. That means that the result you get will be accurate in a floating point sense - the inputs will be regarded as being specified up to a certain number of digits, and the result will be accurate up to a certain number of digits as well.

As with all other floating point computations, you can set the required accuracy by changing the Digits environment variable. If you want to have a reasonable sense of confidence that all n digits that you see in the answer are correct (that is, that this n-digit number provides a closer fit than any other n-digit number - not that it is an n-digit approximation to a perfect fit: such a perfect fit may not exist!), you should set Digits to a higher value initially (maybe n + 3 or so, potentially depending on the number of data points you have), then round the result to n digits. You could do that as follows:

oldDigits := Digits:
Digits := Digits + max(3, ilog10(numelems(X)));
answer := PowerFit(X, Y, x);
Digits := oldDigits:
answer := evalindets(answer, 'float', x -> evalf(SFloat(x)));

That last line is complicated by the split personality of floating point computations in Maple. Maple supports floating point numbers with an arbitrary number of digits, which are implemented in software, and hardware-precision floating point numbers which are more efficient but provide limited accuracy (15 decimal digits on all platforms, currently). It tries to use hardware floats in many situations (including in the Statistics package whenever the setting of Digits allows it), but to round a number to a given accuracy, we have to force it to be a software floating point number.

Another meaning of your question might be that you're looking for a symbolic answer, not a floating point one. In that case, first of all you would definitely not want to have any floating point numbers in your input. Furthermore, the Statistics package can't help you for this case. Thankfully, it's fairly straightforward to implement the required linear algebra yourself:

exactPowerFit := proc(X :: Vector[column],
                      Y :: Vector[column],
                      x :: name,
                      $)
local
    A :: Vector[column],
    M :: Matrix,
    lnX :: Vector[column],
    lnXsum :: algebraic,
    lnY :: Vector[column],
    n :: nonnegint,
    v;

    n := numelems(X);
    if numelems(Y) <> n then
        error "dependent and independent data should have equally many data points, "
        "received: %1 and %2, respectively", numelems(Y), n;
    elif n < 2 then
        error "need at least two data points, received: %1", n;
    elif hastype([X, Y], 'float') then
        WARNING("this procedure is not meant to be used with floating point values");
    end if;

    # Precomputation.
    lnX := ln~(X);
    lnY := ln~(Y);
    lnXsum := add(v, v in lnX);

    M := <n, lnXsum | lnXsum, add(v^2, v in lnX)>;
    A := <add(v, v in lnY), add(lnX[v] * lnY[v], v = 1 .. n)>;

    A := LinearAlgebra:-LinearSolve(M, A);

    return exp(A[1]) * x^A[2];
end proc;

Now you can do

X := <seq(i, i = 1 .. 6)>;
Y := <seq(i^2, i = 1 .. 6)>;
exactPowerFit(X, Y, x);

and get a symbolic result. Or more interestingly,

Y := <seq(i^2 + 1, i = 1 .. 6)>;

and find a large expression for the coefficient and exponent that give the closest fit in a least-squares sense (in the log domain). It's a little too big to print here comfortably.

Hope this helps,

Erik Postma
Maplesoft.

Hi Herclau,

Yes, this can be solved using the Statistics package, but it's not fully automated. It works as follows. First load the package and define a random variable with the appropriate distribution.

with(Statistics):
 X := RandomVariable(Poisson(lambda)):

Your function P(x, lambda) is now available as ProbabilityFunction(X, x).

ProbabilityFunction(X, x);

returns . The sum of probabilities for x = .. a is available as CDF(X, a). For example, you can obtain the first expression you had in your post as CDF(X, 10) and the second expression as 1 - CDF(X, 9).

If I understand correctly, you want to find the value of lambda for which the probability that X > 10 is 0.05, and the value for which the probability that X CDF, but it's the wrong inverse: it finds the value of a resulting in a certain value given a fixed value of lambda, whereas we're looking for a value for lambda given a fixed value of a. (This inverse map is called Quantile.) So we will still need to use fsolve.

alpha := 0.05;
inputs := ;
LuValues := fsolve~(CDF~(X, inputs) -~ alpha, lambda, 0 .. infinity);
LiValues := fsolve~(1 -~ CDF~(X, inputs) -~ alpha, lambda, 0 .. infinity);
PointPlot([LiValues, LuValues], xcoords = convert(inputs, list));

Hope this helps,

Erik Postma
Maplesoft.

PS: Upon rereading, I realized I should subtract 1 from the inputs for computing the LiValues, as follows:

LiValues := fsolve~(1 -~ CDF~(X, inputs -~ 1) -~ alpha, lambda, 0 .. infinity);

For input = 0 the problem becomes solving 0.95 = 0, which of course has no solution. This throws off PointPlot. To fix this:

map[inplace](x -> `if`(nargs = 0, undefined, x), LiValues);

Then call PointPlot as before.

In addition to Joe's and maple fan's ideas, try to make sure that you set the datatype option of your Vectors to float if the entries are indeed floats and you have Digits set to no more than 15. Good diagnostic information is available if you set infolevel[LinearAlgebra] to a sufficiently high number (10 should be more than enough). If I try with two datatype=float vectors after setting infolevel, I get this diagnostic information:

DotProduct:   "calling external function"
DotProduct:   "NAG"   hw_f06eaf

That's probably the fastest that Maple can do - unless your vectors have some special structure (e.g. they're sparse).

Hope this helps,

Erik Postma
Maplesoft.

PS: I would expect that CUDA is not terribly effective for dot products, since the amount of data to be transferred is in the same order of magnitude as the number of elementary operations to be done, and CUDA, or generally GPGPU computing, typically works best when the number of operations is an order of magnitude larger than the number of data transfers - such as matrix multiplication (with near-square matrices).