My workaround is similar to Tom's, and additionally it'll continue to work if and when Maplesoft changes the parameters, regardless of how they're changed (even if it's more complicated than just switching their order):

MyGammaD:= proc(k, theta)
uses St= Statistics;
local a,b, GD:= 'GammaDistribution'(a,b);
   eval(GD, solve({St:-Mean(GD) = k*theta, St:-Variance(GD) = k*theta^2}, {a,b}))
end proc:

I'm not claiming to know whether Maplesoft has ever changed the parameters. The above works regardless of changes.

Even more bothersome to me than this thing with GammaDistribution is that there doesn't seem to be standardization across mathematical software and literature as to whether the second parameter of Normal is sigma or sigma^2.

I don't know what you're trying to achieve with the unevaluation quotes, but they have no effect whatsoever when used on a procedure parameter. Other than that, it looks like you're trying to match expressions that are invariant under op. Any conditional can be made a type by using type satisfies. So, the type is satifies(x-> x = op(x)). This won't error out if op(x) is multi-term, so there's no need for square brackets.

With tools like unapply and ->, there's no practical limit to the level of "meta differential operators" that you can create, each capable of being iterated with @@. In this example, I define an operator Newton which takes an "operator template" (f,c) and returns the differential operator that constructs the Newton's method iteration operator for solving f(x) = c. This final operator can be iterated with @@. Each level of these operators is only one line of code.

Newton:= proc(f,c) local x; subs(_c= c, _c-> unapply(x - f(x,_c)/D[1](f)(x), x)) end proc:
Sqrt:= Newton(((x,c)-> x^2-c), c): 
(Sqrt(2)@@4)(1.5); #1.5 is arbitrary starting point.
                          1.414213562
%^2;
                          1.999999999

The output produced by Newton itself is very ugly to look at. But the result of Sqrt(2) is what you'd expect.

Regardless of whether its number of elements is even or odd, the median^[1] of a nonempty^[2] sorted list X equals

add(X[[floor,ceil]((nops(X)+1)/2)])/2

That idea is used in the following, where I also address an issue that you'll encounter when your input has symbolic real constants that you may not want to clobber with evalf in the output:

#Note that the default sort order doesn't respect the usual order
#of real numbers when the numbers are symbolic:
sin~([$1..4]): (sort, evalf)(%);
     [sin(1), sin(2), sin(3), sin(4)], [0.8414709848, 0.9092974268, 0.1411200081, -0.7568024953]

#We can get the desired order by sorting with evalf:
(x::list)-> ((S,n)-> (S, add(S[[floor,ceil]((n+1)/2)])/2, add(S)/n))(sort(evalf(x)), nops(x)):
   
%(sin~([$1..4]));
     [-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714810

#We can do the same thing plus not clobber the input
#by using sort(..., output= permutation):
(x::list)-> 
   ((P,n)-> (x[P], add(x[P[[floor,ceil]((n+1)/2)]])/2, add(x)/n))
      (sort(evalf(x), 'output'= 'permutation'), nops(x))
:   
%(sin~([1,2,3,4]));
                                      1          1         
    [sin(4), sin(3), sin(1), sin(2)], - sin(3) + - sin(1), 
                                      2          2         

      1          1          1          1       
      - sin(1) + - sin(2) + - sin(3) + - sin(4)
      4          4          4          4 
      
evalf([%]);
     [[-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714811]

^[1]I mean median in the sense that it's taught in elementary math, not the more-technical sense discussed by @sand15 .

^[2]And do we really want to give serious consideration to defining the median of an empty list?

The vertical axis's tickmarks are positioned at -1/2, -3/2, ..., -(2*n-1)/2. Thus, this works:

restart:
n:= 10:
RM:= LinearAlgebra:-RandomMatrix(n):
NewTickMarks:= [($1..n)+~1/2 =~ A||(1..n)]:
Statistics:-HeatMap(RM, tickmarks= [NewTickMarks, (1-lhs=rhs)~(NewTickMarks)]);

You have 4 real-valued dimensions: x1, x2, x3, and f1. You're only allowed 3. If you set one of those 4 to a constant value, then you could plot the remaining 3 as a level surface. Then you could use time as a 4th dimension and make an animation. Sometimes I use gradations of color to represent extra dimensions.

On the left side of the assignment in the double loop, you have uppercase P. I think that you meant lowercase p.

As you read through this explanation, if you get to any point that you don't understand, please ask about it. You need to understand each step before going to the next step.

First you should strive to understand s(x). Using some plots if you need to, convince yourself that:

1. s is periodic with period 1,
2. on its fundamental period 0 <= x <= 1, s is equivalent to piecewise(x <= 1/2, x, x <= 1,1-x).

Let me know if you need help with the plot command.

It follows that s'(x) is also periodic with period 1, and is equivalent on its fundamental period to 

piecewise(x=0, undefined, x < 1/2, 1, x = 1/2, undefined, x < 1, -1, x = 1, undefined). 

Note that s'(x) is always an integer when it's defined. We'll use that fact later. Also, s'(x) is defined for all real x except when x = p/2 for some integer p.

Now define S_M(x) = s(M*x)/M, for any positive integer M. It follows from basic facts about periodic functions (no calculus required) that S_M(x) has fundamental period 1/M. From the chain rule, we have that S_M'(x) =  s'(M*x) wherever it's defined. In particular, it's either 1 or -1 wherever it's defined. And it's defined for all real x except when x=p/(2*M) for some integer p.

Proposition 1: If you have two periodic functions f(x) and g(x) with fundamental periods P and Q and P/Q is an integer, then the period of h(x) = f(x) + g(x) is P. 

[Edit: The red characters have been deleted from the next line.]

Now define b_N(x) = sum(s(2^n*x)/2^n, n= 1..N). Its period is 1/2^N (by that proposition) and it's differentiable for all real x except when x = p/2^(N+1) for some integer p. When it's defined, b_N'(x) is a sum of exactly N values, each being either 1 or -1. So, it's an integer between -N and N, inclusive.

Your original problem was b₄'(9.05), which equals b₄'(.05). The critical points are p/32 for any integer p. The value is S₂'(.05) + S₄'(.05) + S₈'(.05) + S₁₆'(.05) = 1 + 1 + 1 + (-1) = 3 because 1/32 < .05 < 2/32.

Note that we can easily compute b_N'(x) for any positive integer N and any real x by looking at the first N digits after the decimal point in the binary (base-2) expansion of frac(x). You just add 1 for every 0 and subtract 1 for every 1.

It can be created easily like this:

n:= 4:
M:= Matrix(n, scan= band[1,0], [[1$(n-1)], [$1..n]]);

So, the band[1,0] says that there is 1 subdiagonal and 0 superdiagonals, and the initializer is given by those diagonals, in that order.

To optimize the storage, use

M:= Matrix(
   n, 
   shape= triangular[lower], storage= band[1,0],
   scan= band[1,0], [[1$(n-1)], [$1..n]]
);

So, when you use scan= band[...], then you can make your initializer logically follow your diagonal structure, rather than using some complicated indexing scheme that goes by rows and columns.

Yes, it is easy to build a matrix by its submatrices. I usually use the angle-bracket constructors for that. If you give an example, I'll create it.
 

I implemented the L.D.L^T algorithm for hardware-float matrices for you. It's not as fast as NAG code, but it's not bad. I think that you'll find it useful. Let me know.

Would someone who knows more about compiled Maple than I do suggest some improvements for the timing of this?
 

Download LDLT.mw

The seq command does preserve the order. However, you took the output of seq and made it into set by using {...}. Sets do not preserve order, which is consistent with their mathematical definition. If you simply remove those {...}, then you'll get what you expect.

You can assign a default value to a parameter, and thereby make it optional, by directly assigning to it in the proc line:

sim:= proc(x, n, o:= exactly) 

You can do this and also put a type restriction on the parameter like this:

sim:= proc(x, n, o::identical(exactly,minimum,maximum):= exactly)

There's a huge number of different parameter types, and those help pages that you referenced could use a lot more examples.

Why not just use @@ instead of ^ in the first place? 

MyOp:= ((A+B)@@2)(x); #your 1st question

subs(A= f, B= g, MyOp); #your 2nd question

Yes, the function is NumberTheory:-PrimeCounting. If your Maple is too old to have that, identical functionality is provided by numtheory[pi].

I suppose that your F is a PDF. If it is reasonably smooth^[1], then this can be done by Statistics:-Sample. See the third example on its help page, which deals with a custom distribution. I suspect that the efficiency of this example could've been improved by including the basepoints and range options. I'll look into that later.

^[1]Twice differentiable with finitely many inflection points on its support.