Thanks, Joe. I thought there would be a way to do it in a direct way, by summing c(j) over j from 0 to l-1. What's wrong with the sum here? 

Thanks, Joe. I thought there would be a way to do it in a direct way, by summing c(j) over j from 0 to l-1. What's wrong with the sum here? 

Sorry for the subs question. I was trying to compute the probabilities for the values that relative error M (in approximating A with F) can take on. The minimum of such values is zero, the maximum is 2 and there is a finite number of them in between. I failed to interpret this 

f := `assuming`([simplify(Probability(M = x))], [x >= 0]); pr := simplify(subs(p = 1/3, %)); plot(pr, x = -1 .. 3)

so I asked Maple whether M could be described by a means of a generalised density, using the Dirac delta.
Maple returned an expression involving the Dirac delta function,  but I failed to find out how to make use of it or whether it is valid. 

Perhaps I should generalise my question: what would be the standard approach to describe the (absolute) relative error M probabilistically? How to plot this description? 

Sorry for the subs question. I was trying to compute the probabilities for the values that relative error M (in approximating A with F) can take on. The minimum of such values is zero, the maximum is 2 and there is a finite number of them in between. I failed to interpret this 

f := `assuming`([simplify(Probability(M = x))], [x >= 0]); pr := simplify(subs(p = 1/3, %)); plot(pr, x = -1 .. 3)

so I asked Maple whether M could be described by a means of a generalised density, using the Dirac delta.
Maple returned an expression involving the Dirac delta function,  but I failed to find out how to make use of it or whether it is valid. 

Perhaps I should generalise my question: what would be the standard approach to describe the (absolute) relative error M probabilistically? How to plot this description? 

I think this conversion should be done automatically by Maple, perhaps in future versions. If the series converges, I'd like to know it direct from Maple, not from other sources.

I think this conversion should be done automatically by Maple, perhaps in future versions. If the series converges, I'd like to know it direct from Maple, not from other sources.

yes, that the correct result, but what I want is a closed form expression involving h in place of numbers X and Y in  Xalpha +Yalpha^2. Your proc does what mine proc labelled A does, but what I want is an expression in alpha and h only, without sums. The triple sum expression is a minimal example that causes problems.I don't know why.

The expression I'd like maple to provide for me can be obtained manually, by looking up the coefficient sequences at Sloan's OEIS, but I think there must be a way to do it automatically with maple. Is there such a way?   

yes, that the correct result, but what I want is a closed form expression involving h in place of numbers X and Y in  Xalpha +Yalpha^2. Your proc does what mine proc labelled A does, but what I want is an expression in alpha and h only, without sums. The triple sum expression is a minimal example that causes problems.I don't know why.

The expression I'd like maple to provide for me can be obtained manually, by looking up the coefficient sequences at Sloan's OEIS, but I think there must be a way to do it automatically with maple. Is there such a way?   

That's the kind of proof I wamted. For some reason, simplify(eq,trig) with 13.00 considers cos(theta)*(4*cos(theta)^2-3) simpler than cos(3*theta) and report accordingly.   Any general strategy that helps  not overlook useful simplifications?

That's the kind of proof I wamted. For some reason, simplify(eq,trig) with 13.00 considers cos(theta)*(4*cos(theta)^2-3) simpler than cos(3*theta) and report accordingly.   Any general strategy that helps  not overlook useful simplifications?

Yes, a simple typo. Thanks. Should read:
3) one can show, rather generally, that inequality S<sigma^2/mu^2 is equivalent to m1>m2.
 

Yes, a simple typo. Thanks. Should read:
3) one can show, rather generally, that inequality S<sigma^2/mu^2 is equivalent to m1>m2.
 

Recently, we had to think of creating a forecasting wiki.  We didn't consider  the wikibook  line at all. Eventually we chose the mediawiki engine (of so many others) and http://www.siteground.com/wiki-hosting.htm (of so many others). It would be excellent if we have allmaple.info or the like wiki website in the visible future. It is best to separate from wikipedia. NB: It is not an ad or propaganda, I just did some deep research for our own needs.

Thanks for the discussion, Joe. The relavant maple help page does not seem to distinguish between the Pascal and Polya cases of NBD, giving the condition of x>0 without specifying the set of numbers to which it applies. Even with an integer-valued x, it fails to work properly for p<>0.5; try [seq(print(tmp(i, .7)), i = 1 .. 4)], for example. Can the problem be cured locally, without stats package?
 

Thanks to all - very interesting comments and solutions. I've re-coded my original example taking into account your advice - F_4 is clearly best and my F_5 similar. It seems seq and explicit loop structures are similarly efficient for this example? Is there anything else in the code for F_5 or F_4 that slows it down? For example, I read it somewhere that floor, round, ceil and the like are very slow and must be avoided. Thanks. The current performance is as below: F[5] := proc (n) local S, H, s; S := Statistics:-Sample(Statistics:-RandomVariable((':-Poisson')(1))); s := proc (x) options operator, arrow; floor(convert(S(x), `+`)) end proc; H := proc (g) options operator, arrow; map(proc (t) options operator, arrow; `if`(t = 0, 0, s(g)) end proc, g) end proc; [seq(H(floor(S(1)[1])), i = 1 .. n)] end proc testt := proc (n::posint) #F[1]:=FFFF (original);F[2]:=MyF (I'm not sure I understand)[];F[3]:=FF1(scott);F[4]:=FFFF[2] #(acer[]);F[5]:=new: local t, T, i; for i to 5 do t := time(); F[i](n); T[i] := time()-t end do; [seq(F[i] = T[i], i = 1 .. 5)] end proc: > testt(100); testt(100); testt(100); testt(100); testt(100); print(`output redirected...`); # input placeholder [F[1] = 0.547, F[2] = 0.109, F[3] = 0.188, F[4] = 0.015, F[5] = 0.016] [F[1] = 0.219, F[2] = 0.062, F[3] = 0.188, F[4] = 0.016, F[5] = 0.015] [F[1] = 0.188, F[2] = 0.078, F[3] = 0.156, F[4] = 0.016, F[5] = 0.016] [F[1] = 0.203, F[2] = 0.062, F[3] = 0.188, F[4] = 0.016, F[5] = 0.015] [F[1] = 0.219, F[2] = 0.062, F[3] = 0.172, F[4] = 0.016, F[5] = 0.015] > testt(2000); testt(2000); testt(2000); print(`output redirected...`); # input placeholder [F[1] = 4.329, F[2] = 1.781, F[3] = 4.500, F[4] = 0.406, F[5] = 0.344] [F[1] = 4.625, F[2] = 1.781, F[3] = 4.703, F[4] = 0.360, F[5] = 0.437] [F[1] = 4.734, F[2] = 2.016, F[3] = 4.937, F[4] = 0.438, F[5] = 0.375]