@vv Ah, I understand now. Yes, your way is better than mine. Perhaps it'd be worthwhile to issue an error if the integral is significantly different from an integer.

@Carl Love Okay, here's a general procedure for symbolic contour integration (for simple closed contours, counterclockwise, blah, blah). For each point a returned by singular (other than infinity), I compute the numeric contour integral of 1/(z-a) around the contour. If the absolute value of this is significantly^[*1] larger than 0, then a is used.

^[*1]"Significantly" defaults to 10.^(2-Digits), and can be changed by the keyword parameter eps.

ContourInt:= proc(
   f::algebraic, 
   C::(name= algebraic), 
   trange::name= range(realcons), 
   {eps::positive:= 10.^(2-Digits)}
)
local z:= lhs(C), t:= lhs(trange), r:= rhs(C);
   add(
      map2(
         residue,
         f,
         select(
            Z-> abs(evalf(Int(diff(r,t)/(r-eval(z,Z)), trange))) > eps, 
            remove(z-> type(rhs(z), infinity), op~({singular(f,z)}))
         )
      )
   )*(2*Pi*I)
end proc:

Usage for your two cases:

ContourInt(f(sqrt(2),z), z= 2*exp(I*t), t= 0..2*Pi);
ContourInt(f(2,z), z= cos(t)+2*I*sin(t), t= 0..2*Pi);
 

I don't know whether this is more or less complicated or robust than computing the winding number, as suggested by VV above.

@Joe Riel I don't think that Maple V has 2D Input.

@Josolumoh The normalization process that I described won't work for d < 0 because the process requires B >= 0 for all x in the "support" and B > 0 for at least some x. I put "support" in quotes because the technical definition of it is "those x such that B > 0". But if there are some x in the "support" such that B = 0, we can easily work around that, but we can't work around B < 0.

I think that you'll get better search results if you change "cloud" to "cluster".

Do you need all the data to be below the hypothenuse (this seems relatively easy), or just most of it below?

@Josolumoh Yes, I have confirmed that for d = 0 the sum over the support x= 0..infinity is exactly 1 for any nonnegative integer r and any b > 0. Do you wish to proceed with this, using d = 0? 

To proceed for d > 0, the function B would need to be "normalized" for each b by dividing by the sum over the support for that r, d, and b. I'm a bit uncomfortable doing this because I don't understand the significance of the probability function. But it would be a valid probability mass function after the normalization. Perhaps you were thinking of the Beta negative binomial distribution (Wikipedia), also known as the inverse Markov-Polya distribution or the generalized Waring distribution?

@emendes Sorry for putting a mysterious number, 49, into the code. That was my bad. It came from your 50, minus 1. In the code below, I've improved that by renaming your 50 as nspan.

The following code uses a single procedure NestList for hardware floats, software floats, and exact arithmetic:

NestList:= proc(f, x, n::nonnegint, nmlzr:= (), {digits::posint:= trunc(evalhf(Digits))})
uses LT= ListTools;
local 
   k, R:= rtable(0..n, [x]),
   Iterate:= proc(f, R::rtable, n::nonnegint, nmlzr:= (x-> x))
   local k;
      for k to n do R[k]:= nmlzr(f(R[k-1])) od;
      R
   end proc,
   Types:=                [complex(float), radalgnum,    anything],
   nmlzrs:= table(Types=~ [x-> x,          evala@Normal, simplify]),
   NMLZR:= `if`(nmlzr = (), nmlzrs[LT:-SelectFirst((t,x)-> x::t, Types, x)], nmlzr)
;
   Digits:= digits;
   if x::complex(float) and Digits <= trunc(evalhf(Digits)) then
      R:= hfarray(0..n, [x]); 
      try return evalhf(Iterate(f, R, n, x-> x)) catch: print("evalhf failed") end try 
   fi;
   Iterate(f, R, n, NMLZR)
end proc
:
DoPlot:= (aaa::rtable)->
   plot(
      Re~(<<[$-nstep-1..nspan-1]> | aaa^%T >), 
      style= pointline, symbol= solidcircle, symbolsize= 4, 
      thickness= 0, view= [default, 0..1], _rest
   )
:
(nstep, nspan):= (12, 50);
logistic:= x-> 4*x*(1-x):
fp:= [solve((logistic@@5)(x) = x, x)]:
cfp:= allvalues(fp[5]):
p:= nstep+nspan:
plots:-display(<
   DoPlot(NestList(logistic, evalhf(Re(cfp)), p), title= "Using evalhf") |
   DoPlot(
      NestList(logistic, evalf[30](Re(cfp)), p, digits= 30), 
      title= "Using 30 Digits"
   ) |
   DoPlot(NestList(logistic, cfp, p), title= "Using exact arithmetic")
>);

@Josolumoh 

I have some ways to help you, but the function B that you propose as your pmf (or ProbabilityFunction) doesn't sum to 1 for any finite b. Rather, it looks like the limit of its sum is 1 as b approaches infinity, as shown in this plot:

B:= (r,d)-> b-> x-> 
   binomial(x+r-1, x)
   / (1+d*x)
   * (b/2/(1 + d*x + b/2))^r
   * ((1+d*x)/(1+d*x+b/2))^x
;
plot(b-> evalf(Sum(B(3,2.5)(b)(x), x= 0..infinity)), 0..40);

Once you correct this pmf, I'll be happy to help you generate samples.

It's difficult to help you because I can only see your output, not the input that you used to generate it.

There are several easy ways in Maple to map an indexing operation over a set or list of indices.

@Christopher2222 An easier way is to just use Statistics:-Fit instead of Statistics:-LinearFit:

Statistics:-Fit(3*x+a, pts1, pts2, x);

Note that x must be the independent variable and a the unknown parameter. Thus x must be the last argument, not a.

@acer Thank you, and thank you for teaching me about Tabulate. I do really like the way that the default prettyprinting for rtables (which is obviously inherited by DataFrames) sets the column widths and row heights. I don't relish the idea of needing to set the weights. Is there a way (perhaps something in Typesetting) to access the widths that would be used for an rtable? or just the width of an ordinary algebraic expression? Then those could be used (perhaps with some scaling) to set the weights.

It isn't only the lambda that isn't typeset.

@Carl Love Second hint (for "by hand" computation, assisted by Maple): As upsilon --> infinity (either positive or negative), the simple rational function inside the square root, (upsilon - 4)/upsilon, approaches 1. So, just replace it thus:

f:= ...your original expression...:
assume(t > 0);
f1:= simplify(subs(epsilon= 1/upsilon, f));
f2:= simplify(subs((upsilon-4)/upsilon= 1, f1);

@bsoudmand I get manifestly real expressions (with Re and Im resolved). My guess is that you have an earlier version of Maple such that evalc won't automatically "map over" my fancily labelled expressions. Try this:

simplify(evalc~([Re,Im](Ec))) assuming positive;

Kitonum's Answer shows that evalc is not necessary in this case (a rational function is an easy case). However, in general I think it will be necessary.

@acer I was wrong. I was misled by some spurious timings that I can't reproduce.

Instead of cot~(...), you could use map[evalhf, inplace](cot, ...). It would be much faster, but your point is still valid.