# Nested double sum does not converge

I am quite a beginner, and have the following problem. I have a sum that comes from Poisson distribution and for each index there corresponds a binomial-type sum. I try to find its root for the probability that is in the binomial-type sum (as a function of other variables/parameters) so that I could plot it.

But I get a warning that says that Hypergeometric PFQ does not exist and the parameters are not consistent, and that the sum does not converge. I should appreciate any help.

Clear[c, d, a, p]

stochprob[d_, a_, c_] :=

p /. FindRoot[
Sum[Exp[-a]*(a^h)/(h!) Sum[
Binomial[h, k] p^k * (1 - p)^(h - k) 1/(k + 3)^2, {k, 0,
h}], {h, 0, Infinity}] == (1 - d)*c, {p, 1/3}]


I am not sure does the code below produce correct results, but it does provide some numbers reasonably fast. (I guess not because of the negative p, but may be this code is useful anyway...)

I made several changes of the code in the question:

• using NSum instead of Sum;

• calling NSum (or Sum) once with specification of two ranges;

• defining a function, BiSum, that takes numerical arguments only.

Note that I set the method for NSum to something that would produce results faster. It is probably a good idea to experiment with using different methods and/or using different precision goals.

Clear[BiSum]
BiSum[a_?NumericQ, c_?NumericQ, d_?NumericQ, p_?NumericQ] :=
NSum[Exp[-a]*(a^h)/(h!) Binomial[h, k] p^
k*(1 - p)^(h - k) 1/(k + 3)^2, {h, 0, Infinity}, {k, 0, h},
Method -> {"WynnEpsilon", "ExtraTerms" -> 20}, NSumTerms -> 200]

BiSum[1, 1, 1, 1/3]

(* 0.0962654 *)

Clear[c, d, a, p]
stochprob[d_, a_, c_] :=
p /. FindRoot[BiSum[a, c, d, p] == (1 - d)*c, {p, 1/3},
MaxIterations -> 20]

stochprob[0.1, 1, 1]

(* -3.98788 *)