Recursion to compute convolution

Hello,

consider the sum of i.i.d. random variables with :

Is it possible to implement this recursively?

If it is, can anybody tell me how I can do that or give at least a hint?

[Image Analyst]

If you want convolution, simply use the built-inconv() function. Don't do it manually - the hard way.

[Alex Schmdit]

Thank you. But how do I use it here? conv() takes only to input arguments?

[Image Analyst]

I don't know what that formula is, especially the index for P which is S(n-1)==(k - j), but I'm just going by what you say - that you want convolution. There is a built in function for that, that takes several arguments. First one is the main signal. Second argument is the filter weights of the scanning/sliding filter window. Third argument is whether you want the full convolution, only the valid elements, or an output that is the same size as the input signal. Look up the documentation and online tutorials about what convolution is. Find one that has diagrams that show the window as it slides along -- I'm sure there are lots of convolution tutorials out there.

[Alex Schmdit]

I am not sure we are talking about the same thing. If I want to compute the probability of a random variable and know its distribution, I can easy compute it. The distribution of sum of random variable is computed using convolution. For example take . Therefore:

The second factor would be . How do I use conv() here?

Answer 1

[Bruno Luong]

Here is a hint, you must replace XXX and YYY with something (after all it's your homework)

function P = Psum(n, k, alpha)
if k < n
    error('not valid parameter');
end
P1 = @(j) 1/j^alpha - 1/(j+1)^alpha;
if n == 1
    P = XXX;
else
    P = 0;
    for j=1:k-n+1
        P = P + YYY *P1(j);
    end
end
end

[Alex Schmdit]

Thank you for answer. So in the case n==1, XXX=PSum(k-j)^P1(j).

In the other case I get: YYY=PSum(k-j). Is that a possible way which works in Matlab?

[Bruno Luong]

Wrong in both.

You should get (I think)

>> Psum(7,10,2)

ans =

    0.1250

[Alex Schmdit]

Why do you set P=0 in the else part? I has to be YYY=PSum(n-1,k-j,alpha)

[Alex Schmdit]

XXX=P1(j)

[Bruno Luong]

XXX=P1(k)

You are now good!

[Alex Schmdit]

Thank youi very much. Can you explain me why you resest P=0 in the else part?

[Bruno Luong]

To make a P = P + ... in the for-loop returns the right sum at the end.

[Alex Schmdit]

Last question. Why do you put in P1(k) and not j? j would not work. But what is the logic behind this?

[Bruno Luong]

For n = 1, the code must return P(S1 = k) that is P(X = k). You don't sum on j (or you sum on single j=k if you insist)

[Alex Schmdit]

Thank you very much:) Have a nice day.

[Alex Schmdit]

If I use this code to compute different convolutions in an other program, I get the following error:

Out of memory. The likely cause is an infinite recursion within the program.

Why do I get an infinte recursion?

[Bruno Luong]

Do you call with n<=0? or not integer?

[Alex Schmdit]

Yes it happens. Should I add: if(n<=0) return;?

[Bruno Luong]

function P = Psum(n, k, alpha)
if k < n
    error('not valid parameter k');
end
P1 = @(j) 1/j^alpha - 1/(j+1)^alpha;
if n == 0
    % just do it, I have no desire to explain it
    P = double(k == 0);
else
    P = 0;
    for j=1:k-n+1
        P = P + Psum(n-1, k-j, alpha)*P1(j);
    end
end
end

[Alex Schmdit]

Sorry this leads to a negative proability

[Alex Schmdit]

Sorry I get the same.

[Bruno Luong]

I give here a non-recursive solution using CONV suggested by Paul

%
function P = Psum_nonrecurse(n, k, alpha)
% n scalar,  k can be vector
if (n < 1) || any(k < n)
    error('not valid parameter');
end
P1 = @(j) 1./j.^alpha - 1./(j+1).^alpha;
M = P1(1:max(k)-n+1);
P = M;
for i=1:n-1
    P = conv(P, M);
end
P = P(k-n+1);
end

[Alex Schmdit]

Thank you very much for your help. Can you explain to me what the for loop does?

[Bruno Luong]

When you have two independent integer-discrete random variables X and Y, with repectively PDF A and B, the pdf of random variable X+Y is the convolution of A and B:

conv(A,B)

With that in mind, the for-loop computes then sequentially the pdf of

S(i+1) = S(i) + X(i+1).

If you only need to calculate P(S(n)=k) for a scalar k, this methods computes a lot of things and throw away a big part, as oppose to the recursive method. However if you need for a bunch of k, use CONV is a better approach. This function can be called with a vector k (but same n), not the recursive method that can only deal with scalar k.

[Alex Schmdit]

Sorry for the late answer. I think I get your point. Thank you very much. I appreciate your help:)

Answer 2

[Image Analyst]

If dist is your distribution, then the distribution of the sum would be

distSum = conv(dist, dist, 'full');

No recursion needed.

[Alex Schmdit]

What should I fill in for dist?

[Image Analyst]

That is your probability distribution function for your random variable. Just one, not the distribution of the sum of the two but just one alone. For example for a uniform distribution it would be all ones, like

n = 500; % Whatever.
dist = ones(n, 1) / n;
subplot(2, 1, 1);
plot(dist, 'b-', 'LineWidth', 2);
title('PDF (Histogram) of One Variable', 'FontSize', 15);
ylim([0, 0.0025]);
grid on;
distSum = conv(dist, dist, 'full');
subplot(2, 1, 2);
plot(distSum, 'b-', 'LineWidth', 2);
grid on;
title('PDF (Histogram) of Sum of Two Variables', 'FontSize', 15);

where n is the resolution you want to digitize your continuous PDF at.

[Image Analyst]

Looks like you've accepted Bruno's answer, so I guess you got it all solved now.

Answer 3

[Paul]

Let the vector M be the distribution of integer valued X. For example, if X takes on values 1-6 uniformly, then

M = ones(1,6)/6;

Now you can compute the distribution of the sum of the random sample of X as:

Mn = M;
for ii = 2:n
    Mn = conv(Mn,M);
end

This isn't recursive, so I'm not sure if that's the implementation you want. Also, it may not be efficient depending on the values of numel(M) and n because Mn is growing each time through the loop. But I think it gives you the result you seek. You'll have to be careful the indexing of Mn if X can take on values less than 1.

[Paul]

Assuming X can only take on values j, j >= 1, then M(j) = P(X == j). If j can take on values < 1, then you'll have to do some reindexing because Matlab uses only 1-based indexing.