A somewhat more sophisticated approach is to note that if you translate one of the vertices of the tetrahedron to the origin of a Euclidean space, then there is a linear transformation of the unit simplex to the tetrahedron.  Consequently, the volume of the tetrahedron equals the volume of the unit simplex (1/6) multiplied by the determinate of the linear transformation.  The matrix of the linear transformation is easily found, its columns are the three vectors that correspond to the edges extending from the origin of the translated tetrahedron. 

So, if the vertices are the vectors V1, V2, V3, V4, then

M := < V2-V1 | V3-V1 | V4-V1>: # matrix of linear transformation
vol := LinearAlgebra:-Determinant(M)/6:

You may have to take the absolute value if the transformation doesn't preserve orientation

Note that the approximation is only good for small alpha.

The usual way to do that is with zip,

 zip((x,y)->[x[],y], a,b);
                  [[a1, a2, b1], [a3, a4, b2], [a5, a6, b3]]

An interesting alternative is

use T = ListTools:-Transpose in 
   T([T(a)[],b])
end use;
                  [[a1, a2, b1], [a3, a4, b2], [a5, a6, b3]]

Consider the specialization

X1 := RandomVariable(Geometric(1/2)):
Sample(X1, 10);
                                  [4., 0., 1., 0., 1., 0., 1., 0., 0., 2.]

Note that 0 is a valid sample.  ln(0) is undefined, so Mean(ln(X1)) is also undefined.

I had modified your code to correct what may be the issue, but that led to the problem of g5, which I forget to confirm existed in the original, it doesn't, at least not identically.

Here's what I noticed.  When using assumption it is generally necessary to apply them before creating expressions that use the variable, otherwise those expressions get the "unassumed" variable which is different from the variable with assumptions.  To avoid that, I moved your assumptions to the start of the script.  When that is done, I get

indets([{mu = Mean(Y5), sigma^2 = Variance(Y5)},{c,b}]);

                                   {b, c, c, mu, sigma}

Note that c appears twice.  Those correspond to different variables, albeit with the same apparent name. I haven't looked closer, but believe the immediate problem is coming from the use of

assume(GAMMA((1/2)*c) > 0, 0 < sqrt(2^c*GAMMA((3/2)*c)*GAMMA((1/2)*c))*sigma/(2^c*GAMMA((3/2)*c))):

Insert all the values into a Matrix, say, M, then create the spreadsheet with ExcelTools:-Export(M, "mydata.xls");

Don't use square brackets as parentheses, they do not act that way in Maple.  Also, use Pi instead of pi.

That depends on what you mean be flattening.  Depending on the problem, you might be able to do

map(op, L);

or (newer Maple required):

op~(L);

Look at g5.  The mean and variance are not independent, there is no solution (for b and c).

Compilation isn't possible here with the calls to solve.  You might try substituting ?fsolve for solve, which could make the calls faster.  If a and b are simple expressions of x, you might be able to solve everything symbolically then create a procedure that returns the specific solutions given parameter values.

There are two problems with your input.  First, Maple uses square-brackets to delineate lists, not as parentheses.  Second, the differential equation is second order, there is no reason to specify the initial second derivative (dsolve complains). The following works,

f := ((diff(y(x), x))*(diff(y(x), x, x))/sqrt(y(x)^2+diff(y(x), x))+diff(y(x), x, x)) = x^6:
integ := dsolve({f, y(0) = 2, D(y)(0) = 1.5}, numeric, range = 0 .. 10):
integ(1);
                                           d
         [x = 1., y(x) = 3.51232025340609, -- y(x) = 1.59936902917771]
                                           dx

That is a linear operator on the tangent space.  I don't see how you plan to assign a sign to it.  You can use the ?DifferentialGeometry package to compute it,

with(DifferentialGeometry):
DGsetup([x,y,z],M):
f := x^2 + x*y + sin(z);
                                  2
                            f := x  + x y + sin(z)

ExteriorDerivative(f);
                        (2 x + y) dx + x dy + cos(z) dz

This particular example can be handled by declaring and using a discrete-time variable (which only change at events).  Here oldy is declared as a discrete-time float and its value at t=3 is saved.  Then at t=5 it is transferred to y:

restart;
eq := { diff(y(t),t) = y(t) };
ic := { y(0) = 2, oldy(0)=0 };

discvars := [ oldy(t) :: float ]:

ttrip := 5:

evts := [ NULL
          , [ t = ttrip-2,  oldy(t) = y(t) ]
          , [ t = ttrip,    y(t) = oldy(t) ]
        ]:

integ := dsolve( eq union ic, numeric, 'events'=evts, 'discrete_variables' = discvars):

plots:-odeplot(integ, [t,y(t)], 0..6);

You can use the ?gfun package, ?rsolve, and ?sum to handle this:

L := [1, 5, 9, 13, 17]:
# Convert the list to a recurrence equation
R := gfun:-listtorec(L,u(n));
        R := [{-u(n + 2) + 2 u(n + 1) - u(n), u(0) = 1, u(1) = 5}, ogf]

# Use rsolve to solve the recurrence.
u := rsolve(R[1],u);
                                 u := 4 n + 1

# Note that u(0)=1, not u(1)=1.  Let's shift that so v(n)=u(n-1)
v := eval(u, n=n-1);
                                 v := 4 n - 3

# Use sum to compute a symbolic expression for the n-th partial sum; factor it.
V := factor(sum(v, n=1..n));
                               V := n (2 n - 1)

# Plug-in a few values
eval(V, n=1);
                                       1

eval(V, n=100);
                                     19900

# check
add(v, n=1..100);
                                     19900

For finite n you could do

n := 30:
y := add(x[i]^2,i=1..n);
solve({seq(diff(y,x[i]), i=1..n)}, {seq(x[i],i=1..n)});

Thinks get more interesting if you want to take the derivative of an indefinite sum, say

n := 'n':
y := Sum(x[i]^2, i=1..n):
diff(y, x[k]);

You would really like this to return something like

charfcn[1..n](k)*2*x[k]

See ?charfcn for details.  I don't know that there is a package for doing this; maybe in the Application Center.