# set terminal png transparent nocrop enhanced font arial 8 size 420,320 
# set output 'prob2.8.png'
set format x "%2.0f"
set format y "%3.2f"
set key right top Right noreverse noinvert enhanced box linetype -1 linewidth 1.000 samplen 4 spacing 1 width 0 height 0 autotitles
set label 1 "mu" at 6.33333, 0.0170962, 0 left norotate back nopoint
set label 2 "sigma" at 9.31476, 0.0578888, 0 left norotate back nopoint
set arrow 1 from 5.83333, 0, 0 to 5.83333, 0.13038, 0  nohead back nofilled linetype 1 linewidth 1.000
set arrow 2 from 5.83333, 0.0578888, 0 to 8.81476, 0.0578888, 0  nohead back nofilled linetype 1 linewidth 1.000
set samples 200, 200
set xtics border mirror norotate 0.499000,2,17.0000
set ytics border mirror norotate 0.00000,0.0170962,0.170962
set title "negative binomial PDF using gamma approximation" 0.000000,0.000000  font ""
set xlabel "k, x ->" 0.000000,0.000000  font ""
set xrange [ -1.00000 : 18.0000 ] noreverse nowriteback
set ylabel "probability density ->" 0.000000,0.000000  font ""
set yrange [ 0.00000 : 0.170962 ] noreverse nowriteback
Binv(p,q)=exp(lgamma(p+q)-lgamma(p)-lgamma(q))
arcsin(x)=invpi/sqrt(r*r-x*x)
beta(x)=Binv(p,q)*x**(p-1.0)*(1.0-x)**(q-1.0)
bin_s(x)=n!/(n-int(x))!/int(x)!*p**int(x)*(1.0-p)**(n-int(x))
bin_l(x)=exp(lgamma(n+1)-lgamma(n-int(x)+1)-lgamma(int(x)+1)+int(x)*log(p)+(n-int(x))*log(1.0-p))
binom(x)=(n<20)?bin_s(x):bin_l(x)
cauchy(x)=b/(pi*(b*b+(x-a)**2))
chi(x)=exp((0.5*df1-1.0)*log(x)-0.5*x-lgamma(0.5*df1)-df1*0.5*log2)
erlang(x)=lambda**n/(n-1)!*x**(n-1)*exp(-lambda*x)
extreme(x)=alpha*(exp(-alpha*(x-u)-exp(-alpha*(x-u))))
f(x)=Binv(0.5*df1,0.5*df2)*(df1/df2)**(0.5*df1)*x**(0.5*df1-1.0)/(1.0+df1/df2*x)**(0.5*(df1+df2))
g(x)=exp(rho*log(lambda)+(rho-1.0)*log(x)-lgamma(rho)-lambda*x)
geometric(x)=exp(log(p)+int(x)*log(1.0-p))
halfnormal(x)=sqrt2invpi/sigma*exp(-0.5*(x/sigma)**2)
hypgeo(x)=(int(x)>mm||int(x)<mm+n-nn)?0:mm!/(mm-int(x))!/int(x)!*(nn-mm)!/(n-int(x))!/(nn-mm-n+int(x))!*(nn-n)!*n!/nn!
laplace(x)=0.5/b*exp(-abs(x-a)/b)
logistic(x)=lambda*exp(-lambda*(x-a))/(1.0+exp(-lambda*(x-a)))**2
lognormal(x)=invsqrt2pi/sigma/x*exp(-0.5*((log(x)-mu)/sigma)**2)
maxwell(x)=fourinvsqrtpi*a**3*x*x*exp(-a*a*x*x)
negbin(x)=exp(lgamma(r+int(x))-lgamma(r)-lgamma(int(x)+1)+r*log(p)+int(x)*log(1.0-p))
nexp(x)=lambda*exp(-lambda*x)
normal(x)=invsqrt2pi/sigma*exp(-0.5*((x-mu)/sigma)**2)
pareto(x)=x<a?0:b/x*(a/x)**b
poisson(x)=mu**int(x)/int(x)!*exp(-mu)
rayleigh(x)=lambda*2.0*x*exp(-lambda*x*x)
sine(x)=2.0/a*sin(n*pi*x/a)**2
t(x)=Binv(0.5*df1,0.5)/sqrt(df1)*(1.0+(x*x)/df1)**(-0.5*(df1+1.0))
triangular(x)=1.0/g-abs(x-m)/(g*g)
uniform(x)=1.0/(b-a)
weibull(x)=lambda*n*x**(n-1)*exp(-lambda*x**n)
carcsin(x)=0.5+invpi*asin(x/r)
cbeta(x)=ibeta(p,q,x)
cbinom(x)=ibeta(n-x,x+1.0,1.0-p)
ccauchy(x)=0.5+invpi*atan((x-a)/b)
cchi(x)=igamma(0.5*df1,0.5*x)
cerlang(x)=1.0-exp(-lambda*x)*(1.0+lambda*x+0.5*(lambda*x)**2)
cextreme(x)=exp(-exp(-alpha*(x-u)))
cf(x)=1.0-ibeta(0.5*df2,0.5*df1,df2/(df2+df1*x))
cgamma(x)=igamma(rho,x)
cgeometric(x)=(x<1)?geometric(0):geometric(x)+cgeometric(x-1)
chalfnormal(x)=erf(x/sigma/sqrt2)
chypgeo(x)=(x<1)?hypgeo(0):hypgeo(x)+chypgeo(x-1)
claplace(x)=(x<a)?0.5*exp((x-a)/b):1.0-0.5*exp(-(x-a)/b)
clogistic(x)=1.0/(1.0+exp(-lambda*(x-a)))
clognormal(x)=cnormal(log(x))
cnormal(x)=0.5+0.5*erf((x-mu)/sigma/sqrt2)
cmaxwell(x)=igamma(1.5,a*a*x*x)
cnegbin(x)=(x<1)?negbin(0):negbin(x)+cnegbin(x-1)
cnexp(x)=1.0-exp(-lambda*x)
cpareto(x)=x<a?0:1.0-(a/x)**b
cpoisson(x)=1.0-igamma(x+1.0,mu)
crayleigh(x)=1.0-exp(-lambda*x*x)
csine(x)=x/a-sin(n*twopi*x/a)/(n*twopi)
ct(x)=(x<0.0)?0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x)):1.0-0.5*ibeta(0.5*df1,0.5,df1/(df1+x*x))
ctriangular(x)=0.5+(x-m)/g-(x-m)*abs(x-m)/(2.0*g*g)
cuniform(x)=(x-a)/(b-a)
cweibull(x)=1.0-exp(-lambda*x**n)
fourinvsqrtpi = 2.25675833419103
invpi = 0.318309886183791
invsqrt2pi = 0.398942280401433
log2 = 0.693147180559945
sqrt2 = 1.4142135623731
sqrt2invpi = 0.797884560802865
twopi = 6.28318530717959
a = 1.0
alpha = 0.5
b = 2.0
df1 = 1
df2 = 1
g = 1.0
lambda = 0.6
m = 0.0
mm = 25
mu = 5.33333333333333
nn = 75
n = 10
p = 0.6
q = 0.5
r = 8
rho = 3.2
sigma = 2.98142396999972
u = 1.0
xmin = 0
xmax = 17
ymax = 0.170961981379082
xinc = 1
plot negbin(x), g(x - 0.5)