Recursive function for replacing multiple for loops

Hi I want to implement a recursive function which could replace the following code in Matlab:

p=0.2;
n=10;
a=5;
p1=0;

for i = 0:1:(n-a)

      for j = 0:1:(n-i-a)

            for k = 0:(n-i-j-a)

               for l = 0:(n-i-j-k-a)

                   for m = 0:(n-i-j-k-l-a)

                       p1=p1+(p*(1-p)^i)*(p*(1-p)^j)*(p*(1-p)^k)*(p*(1-p)^l)*(p*(1-p)^m);

                   end

                end

            end

       end

end

I could have used the above code if had a=5 or 10. But in my case, value of n is constant like n=100 and value of a can be up to 100, i.e, n>=a, which makes it difficult to change the number of for loops on each value of a. I will be thankful if someone helps me in implementing such a recursive function which could replace the above for loops.

[David Goodmanson]

Hello Ameer, What determines the number of for loops that are done? Is that determined by 'a'? Since i,j,k etc. can be zero, there could be a very large number of for loops.

[David Goodmanson]

I could have asked this question better. What determines the number of independent indices i,j,k, ... that you have? You could have stopped at, say, just two variables and gotten an answer.

Answer 1

[David Goodmanson]

Hello Ameer, suppose you have d for loops, which means d independent variables i,j,k, ... then your expression reduces to p^d sum (1-p)^(i+j+k+ ...) and by computing the number of times a given sum occurs you can get to the following:

N = n-a;
d = 3     %  this is the number of independent i,i,k, ... variables
          %  which is the number of for loops
p2 = 0;
for q = 0:N
  p2 = p2 + (factorial(q+d-1)/(factorial(q)*factorial(d-1)))*(1-p)^q;
end
p2 = p^d*p2

which appears to agree with your calculation.

Answer 2

[Jan]

Why do you want a recursive function? It would be easier to solve this using an index vector:

p  = 0.2;
n  = 10;
a  = 5;
nv = 5;
v  = zeros(1, nv);
p1 = 0;
ready = false;
while ~ready
  p1 = p1 + prod(p * (1-p) .^ v);
  % Update the index vector:
  ready = true;       % Assume that the WHILE loop is ready
  for k = nv:-1:1
    v(k) = v(k) + 1;
    if v(k) <= n - a - sum(v(1:k-1))
      % This is your "0:(n-a -i-j-k-l-...)" criterion
      ready = false;  % No, WHILE loop is not ready now
      break;          % v(k) increased successfully, leave "for k" loop
    end
    v(k) = 0;  % v(k) reached the limit, reset it and iterate v(k-1)
  end
end

The idea is to use one index vector v = [i, j, k, l, m, ...] instead of a bunch of FOR loops. The inner loop "for k" performs the update of this vector and this is not restricted to a certain number of loops.

The power operation is very expensive. If speed matters, calculate the powers once before the loop:

d = (1 - p) .^ (0:n-a);
...
while ~ready
  p1 = p1 + prod(p * d(v+1));

[Ameer Ahmed]

Hi Jan!

How do we choose the value of nv? For example, for p=0.9, n=8, and a=1, what value should be assigned to nv? Should it be equal to the valus of 'a'?

Whenever, main program calls the function of this code, it takes too much time in this function. Can we speed up execution of this code?

Answer 3

[Ameer Ahmed]

Hi David & Jan! Thanks for the response. @David: Number of for loops are determined by a. @Jan: Recursive function is not required. I only need a feasible solution which could replace multiple for loops which are determined by a.

Answer 4

[Ameer Ahmed]

Thanks Jan and David for solving my problem. Great job done by both of you. Many many thanks.