Question: How to define an abstract operator 'MathematicalExpectation'?

Hi,
I would like to present you a recurring problem of mine.

Context:

It is very common, in Modeling and Simulation activities, to account for Uncertainties.

In order to set these ideas down, consider a computational code F (typically a code that solves a set of PDEs in space (M) and time (t) ) and its response Y = F(M,t) (here Y is a short for Y(M,t)).

This response Y usually depends also on some set P of parameters (each of them considered as a scalar quantity) .
A more convenient way to note Y  is Y(M, t | P) where the "|" character is used here to express that Y is considered as a function of M and t for each given value of P.
Generally one does not consider Y as a whole but more often one focuses on some quantity of interest (Q) derived from Y through applying it some operator G (for example the operator max(Y) over some space and time domain).

Applying G generally makes  Q to appear as a function of P alone.
In Uncertainty Quantification activities, a major concern is to understand how uncertainties about P modify the values of Q ?
The classical framework is to consider P as a (multi-dimensionnal) random variable. 
One of the most common problems is then to assess some basic characteristics of Q, where Q is considered as a function of P (a quick and notional notation is Q = H(P) = (G°F)(X, t | P)).

The simpler and faster method to do this is based on a Taylor expansion of H (provided some conditions hold) arround some particular point p* (p* could be the mean of the multi-dimensional distribution of P)
One writes, for every value p of P :

Q = H(p) = H(p*) + (H')^t (p-p*) + 1/2 (p-p*)^t H" (p-p*) ...

where
  H'  is the gradient vector of H according to P at point P = p*
  H" is the hessian matrix of H according to P at point P = p*
  (p-p*) is the vector of differences, assumed to be "small"

Let E the usual "mathematical expectation" operator.
Let us assume  p* denotes the mean of P.
Then, applying E to  the previous formula gives :

E(Q) = E(H(p*)) + E(...) + ....

Assuming some conditions hold, the first order mathematical expectation E(Q) of Q is simply E(H(p*)) = H(p*) = H(E(P))

Some little algebra gives the first order approximation of the variance V(Q) = E(Q^2)-(E(Q))^2 of Q : 

V(Q) = (H')^t V(P) H'  where V(P) is the variance matrix of P

-----------------------------------------------------------------------------------------------------------------------

Problem

I would like to define the operator E so that I could derive automatically approximations of the first 4 moments of Q, up to any desired order.
In particular, order 2 is often necessary as soon as (G°F) is highly non linear regarding P ; and estimations of the 3rd and 4th moments is of great help to determine how much symetric or flat is the distribution of Q.

The idea is to define an operator E with suitable properties and to apply it to a multivariate taylor expansion of Q^n where n is any positive integer

I tried to do this by my own (look to the supplied .mw file) but I do not have sufficient skill in Maple to complete the job.

Could someone help me ?

Even if I am not qualified in saying this, I believe that this type of approximation of the different moments of Q could be included in a future release of Maple ?


Thanks in advance 

Download MathematicalExpectation.mw