why does it take forever to solve a function

im trying to compute value of "pd" for various values of "lnd"..... but it takes too much time without any results and when i terminate i get this message :

Operation terminated by user during sym/symsum (line 65)

In Untitled5 (line 12)
pd=(symsum(nchoosek(M,i)*(pd1)^(i)*(1-pd1)^(M-i), i, s, M));

equations i'm trying to implement are attached

M=8;
N=32;
K=14;
c=2;
s=1;
T=18.623390280430364437245197384232;

      lnd=30*rand;
  pd1=2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T/(1+lnd)))/(gamma(N+1+T/(1+lnd)))-...
      (symsum(nchoosek(N/2,i)*(gamma(K+i)*gamma(N-K-i+1+T/(1+lnd)))/(gamma(N+1+T/(1+lnd))), i, K, (N/2))));
  pd=(symsum(nchoosek(M,i)*(pd1)^(i)*(1-pd1)^(M-i), i, s, M));
  plot(lnd,pd,'b.')

Answer 1

[Walter Roberson]

symsum is often the wrong function to call when you want to calculate a definite sum (that is, the upper and lower bound are numeric.) symsum() works by trying to calculate the general formula for the indefinite summation, and then substituting in the upper and lower bounds. Calculating the general formula can take a long time.

Instead, when you have numeric upper and lower bounds, you should typically generate all of the individual entries; vectorize if possible and otherwise use a for loop or arrayfun(). Then with those in hand, sum() the result.

I corrected some errors in your calculations of pd1

M = 8;
N = 32;
K = 14;
c = 2;
s = 1;
T = 18.623390280430364437245197384232;

nlambda = 50;
lambdavals = linspace(0,30,nlambda);

pd = zeros(1, nlambda);

for L = 1 : nlambda

    lambda = lambdavals(L);

    T1lambda = T./(1+lambda);                 %constant sub-expression factored out
    pd1temp = arrayfun(@(i) nchoosek(N/2,i) * gamma(K+i) * gamma(N-K-i+1+T1lambda), K : N/2, 'Uniform', 0); %constant denominator removed from terms
    pd1sum = sum( [pd1temp{:}] ) ./ gamma(N+1+T1lambda);      %this is the first symsum
    pd1 = 2 * K * nchoosek(N/2,K) * gamma(K) * gamma(N/2-K+1+T1lambda) / gamma(N/2+1+T1lambda) - pd1sum;

    pdtemp = arrayfun(@(i) nchoosek(M,i)*(pd1)^(i)*(1-pd1)^(M-i), s:M, 'Uniform' ,0);
    pdsum = sum( [pdtemp{:}] );                       %this is the second symsum

    pd(L) = pdsum;
  end

plot(lambdavals, pd)

You will find that this will take less than a second to compute.

[Fatma Abdullah]

thank you....but please would you tell me if i want to simulate this function "pd" by monte carlo accept reject method,,,i tried the following but the same problem rises and im not quite sure i could shape it the way you did

if true
 M=8;
N=32;
K=14;
c=2;
s=1;
T=18.623390280430364437245197384232;
lnd=0:5:30;
pd1=(2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T./(1+lnd)))./(gamma(N+1+T./(1+lnd)))-...
(symsum(nchoosek(N/2,i)*(gamma(K+i)*gamma(N-K-i+1+T./(1+lnd)))./(gamma(N+1+T./(1+lnd))), i, K, (N/2)))));
%lnd=0:5:30;
pd=(symsum(nchoosek(M,i)*(pd1).^(i).* (1-(pd1)).^(M-i), i, s, M));
O=10^7;
Mx=1.2*max(pd);
lndi=10*rand(1,O);
pj=10*rand(1,O);
pd1i=2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T./(1+lndi)))./(gamma(N+1+T./(1+lndi)))-...
  (symsum(nchoosek(N/2,i)*(gamma(K+i)*gamma(N-K-i+1+T./(1+lndi)))./(gamma(N+1+T./(1+lndi))), i, K, (N/2))));
f_lndi=(symsum(nchoosek(M,i)*(pd1i).^(i)* (1-(pd1i)).^(M-i), i, s, M));
comp=pj*Mx;
result=f_lndi>=comp;
index=find(result);
lnd_new=lndi(1,index);
[pd1,lnd1]=hist(lnd_new,100);
dlnd=(lnd1(2)-lnd1(1));
h=(pd1/length(lnd_new))/dlnd;
plot(lnd1,h,'b.');
end

[Walter Roberson]

Do not use symsum() for adding definite symbolic terms. Calculate the individual terms and sum() them.

pd1=(2*K*nchoosek(N/2,K)*((gamma(K+i)*gamma(N-K+i+1+T./(1+lnd)))./(gamma(N+1+T./(1+lnd)))-...

If pd1 is related to the image that you posted, then you have coded it incorrectly. I detected multiple errors in your coding and fixed them in my version. In particular, your gamma(N references should be gamma(N/2 references in the first line of your pd1 code.

[Fatma Abdullah]

thank you so much,,,i will try to do that and fix the code