expectation calculation in matlab

weetabixharry said:

Depending on your requirements, either sum() or trapz() are often fine to approximate integrals. (This won't generally allow you to calculate an expectation... just approximate it, assuming ergodicity and stationarity of any random signals).

What are you trying to calculate the expectation of? M? If so, you will need multiple realisations of M. It's easiest to store these in a 3D array, so M,:,1), M,:,2), ... , are the different realisations.

- - - Updated - - -

Oh wait... are you trying to calculate the expectation of a discrete random variable, X? Does P denote probability in your equation?

For a discrete random variable, you don't need to calculate an integral... just the sum in your equation. Let's start with a scalar random variable, Y, before considering the matrix X. In Matlab, we can just write:

Code:

M = 0;
for i = 1 : N
[INDENT]M = M + P(i)*Y(i);[/INDENT]
end

where N is the number of possible discrete outcomes.
Or, we can do exactly the same thing much more compactly using column vectors:

Code:

M = P.'*Y;

As an example, let's consider a fair die. When we roll dice, there are six possible values, each with probability 1/6. So, we have:

Code:

Y = [1;2;3;4;5;6];
P = (1/6)*ones(6,1);
M = P.'*Y

and Matlab will tell you the correct answer that the expected value of Y is 3.5.

Now, we just need to generalise this so X can be a matrix. As I said above, the easiest way to store X is so that X,:,1), X,:,2), ..., are the different realisations. Then, the first way to calculate the expectation is almost the same as before:

Code:

M = zeros(2,2);
for i = 1 : N
[INDENT]M = M + P(i)*X(:,:,i);[/INDENT]
end

If you want to do this very quickly and efficiently for large numbers of realisations, please let me know. This will be a little trickier, as Matlab does not optimise so neatly for higher-dimensional arrays. We could come up with a solution using something like this or this, so it would not be too difficult.

Click to expand...

Y = zeros(2,2);
for k=1:N

....

X(2*k-1:2*k,1:2) = ....;
Y = Y + P(i)*X(2*k-1:2*k,1:2);

....

end

    0         0
    0    1.0000
    0.5000   -0.5000
   -0.5000    0.5000
    0.5000    0.5000
    0.5000    0.5000
    1.0000         0
         0         0
    0.5000    0.5000
    0.5000    0.5000
    0.5000    0.5000
    0.5000    0.5000
    0.5000    0.5000
    0.5000    0.5000

weetabixharry said:

the easiest way to store X is so that X,:,1), X,:,2), ..., are the different realisations

Click to expand...

X = zeros(2,2,N);
Y = zeros(2,2);
for k=1:N

....

X(:,:,k) = ....;
Y = Y + P(k)*X(:,:,k);

....

end

X = zeros(2,2*N);
Y = zeros(2,2);
...

weetabixharry said:

In my first post I said that:

So, let me clarify that. In Matlab, we want to create a 3D array of size (2 x 2 x N) to store the N different realisations of X:

Code:

X = zeros(2,2,N);
Y = zeros(2,2);
for k=1:N

....

X(:,:,k) = ....;
Y = Y + P(k)*X(:,:,k);

....

end

You will find that this gives a slightly more convenient output in Matlab.

There was nothing technically wrong with your approach and it is actually really useful to understand how to do that type of "sliced" indexing (e.g. 2*k-1:2*k). The only comment I would have is that you should always pre-allocate your arrays to the correct size:

Code:

X = zeros(2,2*N);
Y = zeros(2,2);
...

Otherwise, your code will execute slowly as memory is reallocated every time the array grows. (And also I think your P(i) should be P(k)).

Click to expand...

for k=1:N
X(k) = [(2x2matrix)]
end

for k=1:N
X(2*k-1:2*k,1:2) = [(2x2matrix)]
end