how to implement an infinite series?

im trying to implement the function in the attached image but i get that error "

Attempt to reference field of non-structure array.

at the step where (v) is being calculated. i dont know what it means and i dont know how to fix it

if true
clc;
close all;
N=32;
K=14;
c=0.8;
p=1;
pc=p*gamma(1+2/c);
T=399.01923099362968691492508329876;
s=T/(pc*(1+lambda);
syms i j f
nlambda = 50;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
for L = 1 : nlambda
  lambda = lambdavals(L);
 for i=K:N/2
     w=nchoosek(N/2,i) %first series over i
     z=arrayfun(@(j)nchoosek(K+i-1,j)*(-1)^(j),0 :K+i-1, 'uniform',0)%second series over j
     v=symsum((-1)^(f)*(N+j-K-i+1)^(f)*gamma((f+1)*c/2)/(factorial(f)*(p*(T/(pc*(1+lambda))))^((f+1)*c/2)),i ,0, Inf)%third series over f
      pd1temp=w*z*v
       end
     pd1 = c * K * nchoosek(N/2,K) * pd1temp;   %the whole function
     end
end

[Jan]

No useful comments, outcommented code, no auto-indentation, much code not related to the question, incomplete error message which does not mention, which line is failing: This makes it hard to post an answer. Please invest some time to narrow down the question. Post the relevant part of the code only.

[Fatma Abdullah]

any ideas????!!!

[Fatma Abdullah]

HELP !

[Walter Roberson]

Are there any constraints on the values of c, k, or N ? For example is N constrained to be an even positive integer? Is k constrained to be a non-negative integer less than or equal to N/2 ? Is c constrained to be between 0 and 1?

[Fatma Abdullah]

YES, N is constrained to be positive from 0 to 32.....k also constrained to be positive from 0 to 14 ....c could be any positive value from 0 to 2 including fractions

[Walter Roberson]

So N could be non-integral or could be odd-valued, so N/2 is not necessarily an integer? Should

 (N/2)
 ( k )

be understood as the generalization

gamma(N/2+1)/(gamma(k+1)*gamma(N/2-k+1))

instead of as

factorial(N/2) / (factorial(k) * factorial(N/2-k))

?

And to confirm, k could be (for example) pi, since 0 < pi < 14?

I notice in your code that you set k = 14. You say that "positive from 0 to 14" and since 0 is not positive, that would suggest that your bounds are "exclusive" rather than "inclusive", the interval (0,14) instead of [0,14] -- but if so then k could not be 14. Please clarify.

Is N required to be an integer? Is N required to be an even integer, which is to say an integer that is a multiple of 2? Is k required to be an integer? Can k be 0? Can c be 0? Can c be 2 ?

[Fatma Abdullah]

N is always integer positive even number,so yes N is an integer that is a multiple of 2, k is required to be an integer positive between value>0 and N/2 and could be even or odd so the limits is ]0,N/2] , c is limited between value>0 and infinity including fractions and is required to be positive so the limits of c is ]0,inf[

[Fatma Abdullah]

was that helpful ?!!,,, um wiling to clear any further questions

[Walter Roberson]

I have not been able to find any closed form solution for the infinite sum. It looks to me as if numeric approximation might work well, but numeric approximation is not the solution to an infinite summation.

[Fatma Abdullah]

okay....what can i do then?....can you suggest me any solutions i can try?

[Walter Roberson]

The denominator of your summation includes rho and s to the power of (f+1)*c/2 with f growing to 0. You could imagine that c might be small, but since it is positive, no matter how small it is, (f+1)*c/2 is going to grow indefinitely large as f approaches infinity. So we have a denominator that includes (rho*s) to an indefinitely large power as f approaches infinity.

Now, if rho*s is greater than 1, then as the power (f+1)*c/2 grows indefinitely large, so would (rho*s)^((f+1)*c/2) . That would lead you to an indefinitely large positive denominator. If the numerator of the fraction is growing more slowly then the implication is that the fraction as a whole would converge to 0, allowing for the possibility of a finite summation.

But if rho*s is less than 1, then as the power (f+1)*c/2 grows indefinitely large, then (rho*s)^((f+1)*c/2) would grow indefinitely small. If the numerator of the fraction is not growing smaller at a faster rate, then this would lead to an indefinitely small denominator which would explode, leading to divergence of the sum.

Can (rho*s) be less than 1? If you work through your equations, rho*s = 399.01923099362968691492508329876 / (gamma(2/c + 1)*(lambda + 1)) . Your c is permitted to range up to infinity, so 2/c can decrease indefinitely towards, so gamma(2/c+1) can head towards gamma(1) = 1. So what about lambda? We do not have an expressed limit on lambda so far; looking at the way lambdavals is constructed, it would appear that lambda would always be at least 1 and could be at least 1000, possibly more, so that gives you 399/1000 which is less than 1, so Yes, rho*s can indeed be less than 1. gamma(2/c+1) would have to be less than about 1/4 to prevent that, and that combination cannot happen for non-negative c.

Maybe we had better have another look at the numerator. The N+ etc looks non-negative, so that to the power f is going to grow indefinitely large. gamma((f+1)*c/2) is going to grow indefinitely large. With those two growing indefinitely large, if the denominator grows indefinitely small, we have divergence.

So back to the denominator then, where I have neglected the factorial(f) . Perhaps that could be shown to grow more quickly than (rho*s)^((f+1)*c/2) falls.... ah yes, it appears that for constant rho*s that is indeed the case, that factorial(f) grows faster than (rho*s)^((f+1)*c/2) decreases. So maybe we do not need to worry? But then there is the question of whether the gamma in the numerator cancels that out... I can see circumstances in which it can...

At the moment it is looking to me as if the sum can be divergent, but it needs more examination.

[Walter Roberson]

... Yes, with a few test values I can establish that the sum is divergent for c >= 2 and larger lambda > 1, and I can establish that the terms of the sum converge towards 0 for c < 2 and larger lambda. At the moment I am not certain of the behaviour for lambda = 1 exactly.

[David Goodmanson]

Hello fatma and Walter,

I took a look at just the inner, infinite series. For c = 0, the series diverges simply because gamma(0) = inf. Aside from that, and with the help of the large argument expression for the gamma function, for large n the ratio between successive terms is

|A_n+1 / A_n| = B (c/2)^(c/2) n^(c/2 -1)
where
B = (N+j-k-i+1)/(rho s)^(c/2).        [modified again, second time]

The ratio test says the series is

divergent for      c > 2
convergent for 0 < c < 2
for c = 2:
  divergent for    |B| > 1
  convergent for   |B| < 1
  undetermined for |B| = 1

Hopefully the numerical results are in accord with this. How quickly the series converges is of course another matter.

[Walter Roberson]

Thank you, David.

[David Goodmanson]

I corrected the comment above to change the expression for B to include ^(c/2), and the convergent-divergent conclusions are the same.

[David Goodmanson]

I wanted to mention that you can replace the infinite sum by a (numerical) integral from 0 to inf. In some cases the sum needs a huge number of terms before the terms start to get smaller, and I think the integral may work over a wider range of parameter values than the sum does.

If you use the identity

gamma(a) = Integral {0,infinity} x^(a-1) e^(-x) dx

on the gamma function in the numerator and sum the series, then add and subtract a known integral, you can arrive at

 c2 = c/2;
 B = (N+k-j-i+1)/(rho s)^(c/2)

 result = (1/(rho*s)^c2)*( (1/(B*c2)) + integral(@(x)b(B,c2,x),0, inf) )

function y = b(B,c2,x)
  y = (exp(-x)-1).*x.^(c2-1).*exp(-B*x.^c2);
  y(x==0) = 0;
end

qwrqwerth a

[Fatma Abdullah]

i think i do understand your contributions; but unfortunately i still haven't been able to understand how to numerically apply this function on matlab

[Walter Roberson]

You had indicated that c can be 0 to infinity. We calculate that if c > 2 then the summation is infinite and so does not need to be carried out. We calculate that if c == 2 exactly then under some circumstances the summation is infinite and so does not need to be carried out. Therefor there is no point doing the summation unless c < 2, or unless c == 2 exactly and a particular relationship holds on some of the other variables.

[Walter Roberson]

David, you wrote,

B = (N+k-j-i+1)/(rho*s)^(c/2). [modified]

for c = 2: divergent for |B| > 1

That test looks at any one B value, but the B values change because there are two other levels of summation happening that involve changing some of those values. I just wanted to check some things along those lines:

In the (N+k-j-i+1) component, N and k are fixed in all levels of the summation, but i and j can increase relative to their starting values; an increase of i or j would give a smaller numerator, which would lead to a smaller B value (provided the terms are non-negative, which should be the case given the constraints.) Therefore we should be able to test the first combination for B, and if it gives a value < 1 then it is not possible for any of the other possibilities to end up with B > 1.

The starting combination is i = k and j = 0, which would lead to (N+k-0-k+1) = N+1 . So looking overall, if c == 2 and (N+1)/(rho* s)^(c/2) > 1 then there would be divergence. Taking into account that this is the c==2 case, that could be rewritten as N > (rho*s) -1 would give divergence.

Let's see... for c = 2, s works out as T/(rho*(1+lambda)) so rho*s works out as T/(1+lambda), so N > T/(1+lambda) - 1 . The right hand side is least for larger lambda . We are not given a hard constraint on lambda but according to the code it can be up to 10^(30/10) = 1000, and T is about 400. For simplicity let lambda = 999 so 1+lambda = 1000 and approximate T as 400, so the constraint becomes N > 1/4 - 1 so Yes, there would be a problem . The largest supported lambda would become T/(N+1)-1 which would be about 12 for N = 32 and the given T.

On the other hand, the code is evaluating the summation for multiple lambda, so for the given N and T ranges, this mostly tells us that we can stop going fairly early under this circumstance of c = 2 exactly.

[Walter Roberson]

The user sent me a copy of the paper, which is "Performance analysis of ordered-statistic greatest of-constant false alarm rate with binary integration for M-sweeps", by X.W. Meng. Glancing at it, two things stand out soon:

1. Yes, c=2 is a case of interest; and
2. they talk about doing "binary integration", which is not something I am familiar with:

"Binary integration is a non-coherent approach offering most of the benefits of integration with a simpler and less expensive receiver implementation [7]. This technique is based on detecting a specified number (S) of the available pulses (M). If at least S of the pulses are detected, a ‘hit’ is declared. This is sometimes referred to as coincidence or S-out-of-M detection."

I have not read deeply enough into the paper to understand how that affects the situation.

[David Goodmanson]

Hi Walter, thanks for your continued interest in this. Unfortunately I seem to be losing it trying to write the expression for B, which should not be that hard. Per the numerator in the original posting, it's actually

B = (N+j-k-i+1)/(rho s)^(c/2).        [modified again]
not
B = (N+k-j-i+1)/(rho s)^(c/2).        [modified]

so i and j come in with opposite sign and you can't make the same conclusion that you have. I apologize if this costs you some time.

When c=2 the gamma functions go away and you get a straightforward geometric series that converges for |B|<1. A worse problem is when c < 2 and |B| is only moderately large. This code shows log10 of the individual terms:

n = 0:63000;
c = 1.8;
B = 3;
y = .434*(n*log(B) + gammaln((n+1)*(c/2)) - gammaln(n+1));
plot(n,y)

and it tops out at around n = 23,000 where the value of the term itself is roughly 1e992, eventually dropping to 1e-20 at around n = 62,600. Not a formula for success, so |B| had better be small. I can't vouch for its accuracy but the integral method does give a normal looking answer here, .2847.

[Fatma Abdullah]

Walter....binary integration is a method used in the paper and its a definition "binary integration" made to explain a certain method of combining the back scattered signals from a target in the radar receiver , its not a special type of mathematical integration as it might first come to your mind......

[Fatma Abdullah]

I've forgot to mention that rho would actually be cancelled as described in the attached img......deeply sorry for not mentioning that earlier

I've tried to code differently as:

if true
  clc;
close all;
M=8;
N=32;
K=14;
c=1.2;
pc=p*gamma(1+2/c);
s1=1;
T=54.199575;
syms f
nlambda = 30;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
pd = zeros(1, nlambda);
for L = 1 : nlambda
lambda = lambdavals(L);
Gsum=0;
S=T*p/(p*gamma(1+2/c)*(1+lambda))
pc=p*gamma(1+2/c);
for i=K:N/2
   i;
   w=nchoosek(N/2,i); %first series over i
   z=0;
   sumjf1=0;
   for j1=0:K+i-1
     z=(nchoosek(K+i-1,j1)*((-1)^j1));%second series over j1
   v1=subs((((-1)^(f))*((N+j1-K-i+1)^(f))*gamma((f+1)*c/2)/(factorial(f)*((T/(gamma(1+2/c)*(1+lambda)))^((f+1)*c/2)))),f,0,inf);%infinite series over f
   v1sum=sum(v1);
    pd1temp=v1sum*z
    sumjf1=sumjf1+pd1temp;
   end
   sumjf=sumjf1*w;
 Gsum=Gsum+sumjf ;
end
   pd1 =1*c* K * nchoosek(N/2,K) * Gsum    %the whole function
    pdtemp = arrayfun(@(i) nchoosek(M,i)*(pd1)^(i)*(1-pd1)^(M-i), s1:M, 'Uniform' ,0);%pd1 is implemented in this function of pd to draw pd against lambda(db)
  pdsum = sum( [pdtemp{:}] ); 
              pd(L) = pdsum;
end
plot(lambdavals2, pd);hold on
axis([0 30 0 1])

end

[Fatma Abdullah]

When i draw this equation it was totally different from what i expected( almast a constant straight line with pd=0 for all values of lambda )what i expected is in the attached img (the curve I've pointed)

[Walter Roberson]

Your original code had s=T/(pc*(1+lambda); which is missing a ')' and so is ambiguous. The image you just posted has T divided by rho_c all multiplied by 1+lambda but you appear to have implemented as T divided by (rho_c times 1+lambda)

[Fatma Abdullah]

its T/[pc*(1+lambda)] and its the same in the image just a different way of writing ....i might have missed the parenthesis in the code

[Fatma Abdullah]

and the value of T is not always 399 it could be any value >0 to thousands .....its calculated according to other equations mentioned in the paper

Answer 1

[Walter Roberson]

You have

i=13;

which assigns a specific numeric value to i . Then, a few lines later, you have

syms pfa1 i T j f

which overwrites i with a symbolic value, equivalent to doing

i = sym('i');

So now i is symbolic with no value other than itself, and has completely lost track of the fact it had previously been assigned 13.

Then later in your code you assume that i is still 13.

[Fatma Abdullah]

no i need (i) as a symbolic value actually, i assigned i=13 just to get rid of that warning of considering it complex,its stupid i know, and my code in the question wasn't well written but i've edited it now

[John D'Errico]

But you still have i defined as symbolic AND as a loop variable.

[Fatma Abdullah]

okay ....i've edited the code to be like that and it takes literally forever,,,is there any way to fix it to run faster and get me results

if true
  M=8;
N=32;
K=14;
c=0.8;
SS=1:8;
i=13;
p=1;
pc=p*gamma(1+2/c);
T=399.01923099362968691492508329876
nlambda = 50;
lambdavals2 = linspace(0,30,nlambda);
lambdavals = 10.^(lambdavals2/10);
pd = zeros(1, nlambda);
stop_count = nlambda;
for L = 1 : nlambda
  lambda = lambdavals(L);
  pd1temp=0;
  s=T./pc.*(1+lambda)
 for i=K:N/2
     w=nchoosek(N/2,i); %first series over i
     for j=0:K+i-1
        z= nchoosek(K+i-1,j)*(-1)^(j);
    % z=arrayfun(@(j)nchoosek(K+i-1,j)*(-1)^(j),0 :K+i-1, 'uniform',0)%second series over j
    % v=arrayfun(@(f)
     v=symsum((-1)^(f)*(N+j-K-i+1)^(f)*gamma((f+1)*c/2)/(factorial(f)*(p*s)^((f+1)*c/2)),f ,0, Inf);%third series over f
      pd1temp1=w*z*v;
      pd1temp=pd1temp1+pd1temp;
     end
            pd1 = c * K * nchoosek(N/2,K) * pd1temp
     plot(pd1,lambdavals2);hold on
 end
end

end

[Torsten]

It does not make sense to use "symsum" if no symbolic variables are defined.

Best wishes

Torsten.

[Fatma Abdullah]

it doesnt matter even if i put the line

if true
  syms f

end

still nothing changes,,,it takes forever, nothing comes out, and eventually matlab just doesnt respond and shuts down !!!

[Walter Roberson]

You asked for an infinite sum; that can potentially take a long time and a lot of memory unless convergence can be established.

That is why I have been asking the questions about the constraints: because knowing the constraints can make it much easier to establish whether something converges or not.

[Fatma Abdullah]

N is always integer positive even number,so yes N is an integer that is a multiple of 2, k is required to be an integer positive between value>0 and N/2 and could be even or odd so the limits is ]0,N/2] , c is limited between value>0 and infinity including fractions and is required to be positive so the limits of c is ]0,inf[