--- /dev/null
+#
+# $Id: prob.dem,v 1.9 2006/06/14 03:24:09 sfeam Exp $
+#
+# Demo Statistical Functions version 2.3
+#
+# Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl
+
+# History:
+# -- --- 1992 Jos van der Woude: 1st version
+# 06 Jun 2006 Dan Sebald: Added some variety and plotting techniques for
+# better visual effect. More tutorial in nature.
+
+print " Statistical Library Demo, version 2.3"
+print "\n Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl"
+print "\n\n\n\n\n\n\n"
+print "NOTE: contains 54 plots and consequently takes a lot of time to run"
+print " Press Ctrl-C to exit right now"
+pause -1 " Press Return to start demo ..."
+
+load "stat.inc"
+
+eps = 1.0e-10 # Supposed to be float resolution (nice if were defined internally)
+
+## Gamma function
+xmin = -5.5
+xmax = 5
+ymin = -10
+ymax = 10
+unset key
+set xzeroaxis
+gsampfunc(t,n) = t<0?0.5*1/(-t+1.0)**n:1.0-0.5*1/(t+1.0)**n
+set parametric
+set trange [-1:1]
+set sample 200
+set xrange [xmin : xmax]
+set yrange [ymin : ymax]
+set xlabel "x"
+set ylabel "gamma(x)"
+set arrow 1 from 0,ymin to 0,ymax nohead lt 0
+set arrow 2 from -1,ymin to -1,ymax nohead lt 0
+set arrow 3 from -2,ymin to -2,ymax nohead lt 0
+set arrow 4 from -3,ymin to -3,ymax nohead lt 0
+set arrow 5 from -4,ymin to -4,ymax nohead lt 0
+set arrow 6 from -5,ymin to -5,ymax nohead lt 0
+set title "gamma function, very useful function for probability"
+plot gsampfunc(5*t,5)-6, gamma(gsampfunc(5*t,5)-6) lt 1, \
+ gsampfunc(5*t,5)-5, gamma(gsampfunc(5*t,5)-5) lt 1, \
+ gsampfunc(5*t,4)-4, gamma(gsampfunc(5*t,4)-4) lt 1, \
+ gsampfunc(5*t,3)-3, gamma(gsampfunc(5*t,3)-3) lt 1, \
+ gsampfunc(5*t,2)-2, gamma(gsampfunc(5*t,2)-2) lt 1, \
+ gsampfunc(5*t,1)-1, gamma(gsampfunc(5*t,1)-1) lt 1, \
+ 5*gsampfunc(5*t,2), gamma(5*gsampfunc(5*t,2)) lt 1
+pause -1 "Hit return to continue"
+ymin = ymin/2
+ymax = ymax/2
+set yrange [ymin:ymax]
+set ylabel "lgamma(x)"
+set arrow 1 from 0,ymin to 0,ymax nohead lt 0
+set arrow 2 from -1,ymin to -1,ymax nohead lt 0
+set arrow 3 from -2,ymin to -2,ymax nohead lt 0
+set arrow 4 from -3,ymin to -3,ymax nohead lt 0
+set arrow 5 from -4,ymin to -4,ymax nohead lt 0
+set arrow 6 from -5,ymin to -5,ymax nohead lt 0
+set title "log gamma function, similarly very useful function"
+plot gsampfunc(5*t,5)-6, lgamma(gsampfunc(5*t,5)-6) lt 1, \
+ gsampfunc(5*t,5)-5, lgamma(gsampfunc(5*t,5)-5) lt 1, \
+ gsampfunc(5*t,4)-4, lgamma(gsampfunc(5*t,4)-4) lt 1, \
+ gsampfunc(5*t,3)-3, lgamma(gsampfunc(5*t,3)-3) lt 1, \
+ gsampfunc(5*t,3)-2, lgamma(gsampfunc(5*t,3)-2) lt 1, \
+ gsampfunc(5*t,3)-1, lgamma(gsampfunc(5*t,3)-1) lt 1, \
+ 5*gsampfunc(5*t,3), lgamma(5*gsampfunc(5*t,3)) lt 1
+pause -1 "Hit return to continue"
+reset
+
+# Arcsinus PDF and CDF
+r = 2.0
+mu = 0.0
+sigma = r / sqrt2
+xmin = -(r+1)
+xmax = r+1
+unset key
+set zeroaxis
+set xrange [xmin : xmax]
+set yrange [-0.2 : 1.2]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 50*xmax+1
+set title "arcsin PDF with r = 2.0"
+plot arcsin(x, r)
+pause -1 "Hit return to continue"
+set title "arcsin CDF with r = 2.0"
+set yrange [-0.2 : 1.2]
+plot carcsin(x, r)
+pause -1 "Hit return to continue"
+
+# Beta PDF and CDF
+p = 5.0; q = 3.0
+mu = p / (p + q)
+sigma = sqrt(p**q) / ((p + q ) * sqrt(p + q + 1.0))
+xmin = 0.0
+xmax = 1.0
+#Mode of beta PDF used
+ymax = (p < 1.0 || q < 1.0) ? 2.0 : 1.4 * beta((p - 1.0)/(p + q - 2.0), p, q)
+set key right box
+set zeroaxis
+set xrange [xmin : xmax]
+set yrange [0 : ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 200
+set title "beta PDF"
+plot beta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \
+ beta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \
+ beta(x, 0.5, 2.5) title "p = 0.5, q = 2.5"
+pause -1 "Hit return to continue"
+set yrange [0:1.1]
+set title "incomplete beta CDF"
+set key left box
+plot cbeta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \
+ cbeta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \
+ cbeta(x, 0.5, 2.5) title "p = 0.5, q = 2.5"
+pause -1 "Hit return to continue"
+
+# Binomial PDF and CDF
+n = 25; p = 0.15
+mu = n * p
+sigma = sqrt(n * p * (1.0 - p))
+xmin = int(mu - 4.0 * sigma)
+xmin = xmin < -2 ? -2 : xmin
+xmax = int(mu + 4.0 * sigma)
+xmax = xmax < n+2 ? n+2 : xmax
+ymax = 1.1 * binom(int((n+1)*p), n, p) #Mode of binomial PDF used
+unset key
+unset zeroaxis
+set xrange [xmin : xmax]
+set yrange [0 : ymax]
+set xlabel "k ->"
+set ylabel "probability density ->"
+set ytics 0, ymax / 10, ymax
+set format x "%2.0f"
+set format y "%3.2f"
+set sample (xmax - xmin) + 1
+set title "binomial PDF with n = 25, p = 0.15"
+plot binom(x, n, p) with impulses
+pause -1 "Hit return to continue"
+set ytics autofreq
+set xzeroaxis
+set title "binomial CDF with n = 25, p = 0.15"
+set yrange [-0.1 : 1.1]
+set ytics 0, 0.1, 1.0
+plot cbinom(x, n, p) with steps
+pause -1 "Hit return to continue"
+
+# Cauchy PDF and CDF
+a = 0.0; b = 2.0
+#cauchy PDF has no moments
+xmin = a - 5.0 * b
+xmax = a + 5.0 * b
+ymax = 1.1 * cauchy(a, a, b) #Mode of cauchy PDF used
+set key left box
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100
+set title "cauchy PDF"
+a=0
+b=2
+plot [xmin:xmax] [0:ymax] cauchy(x, 0, 2) title "a = 0, b = 2", \
+ cauchy(x, 0, 4) title "a = 0, b = 4"
+pause -1 "Hit return to continue"
+set title "cauchy CDF"
+plot [xmin:xmax] [0:1.0] ccauchy(x, 0, 2) title "a = 0, b = 2", \
+ ccauchy(x, 0, 4) title "a = 0, b = 4"
+pause -1 "Hit return to continue"
+
+# Chi-square PDF and CDF
+k = 4.0
+mu = k
+sigma = sqrt(2.0 * k)
+xmin = mu - 4.0 * sigma
+xmin = xmin < 0 ? 0 : xmin
+xmax = int(mu + 4.0 * sigma)
+k = 2.0
+ymax = (k > 2.0 ? 1.1*chisq(k - 2.0, k) : 0.5) #Mode of chi PDF used
+set key right box
+set zeroaxis
+set xrange [xmin+eps : xmax] #Discontinuity at zero for k < 2
+set yrange [0:ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100
+set title "chi-square PDF"
+set key right box
+set samples 15*20+1
+keystr(k) = sprintf("k = %d", k)
+plot k = 1, x==0?1/0:chisq(x, k) title keystr(k), \
+ k = 2, x==0?1/0:chisq(x, k) title keystr(k), \
+ k = 3, chisq(x, k) title keystr(k), \
+ k = 4, chisq(x, k) title keystr(k), \
+ k = 5, chisq(x, k) title keystr(k), \
+ k = 6, chisq(x, k) title keystr(k), \
+ k = 7, chisq(x, k) title keystr(k), \
+ k = 8, chisq(x, k) title keystr(k)
+pause -1 "Hit return to continue"
+set yrange [0:1.1]
+set key bottom right box
+set title "chi-square CDF"
+plot k = 1, cchisq(x, k) title keystr(k), \
+ k = 2, cchisq(x, k) title keystr(k), \
+ k = 3, cchisq(x, k) title keystr(k), \
+ k = 4, cchisq(x, k) title keystr(k), \
+ k = 5, cchisq(x, k) title keystr(k), \
+ k = 6, cchisq(x, k) title keystr(k), \
+ k = 7, cchisq(x, k) title keystr(k), \
+ k = 8, cchisq(x, k) title keystr(k)
+pause -1 "Hit return to continue"
+
+# Erlang PDF and CDF
+lambda = 1.0; n = 2.0
+mu = n / lambda
+sigma = sqrt(n) / lambda
+xmax = int(mu + 5.0 * sigma)
+n = 1.0
+ymax = n < 2.0 ? 1.0 : 1.1 * erlang((n - 1.0) / lambda, n, lambda) #Mode of erlang PDF used
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set title "erlang PDF"
+set key top right box
+l1 = 1.0; l2 = 0.5
+set arrow 1 from 2,0.8 to 0.33,erlang(0.33,1,l1)
+set arrow 2 from 2,0.8 to 0.33,erlang(0.33,1,l2)
+set label 1 "n = 1, exponential r.v." at 2.1,0.8 left
+keystr(n,lambda) = sprintf("lambda = %0.1f, n = %d", lambda, n)
+plot [0:xmax] [0:ymax] n = 1, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \
+ n = 1, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda), \
+ n = 2, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \
+ n = 2, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda)
+pause -1 "Hit return to continue"
+unset label 1
+unset arrow 1; unset arrow 2
+set title "erlang CDF"
+set key bottom right box
+plot [0:xmax] [0:1.1] n = 1, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \
+ n = 1, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda), \
+ n = 2, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \
+ n = 2, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda)
+pause -1 "Hit return to continue"
+
+# Thanks to mrb2j@kelvin.seas.Virginia.EDU for telling us about this.
+# Extreme (Gumbel extreme value) PDF and CDF
+alpha = 1.0; u = 0.0
+mu = u + (0.577215665/alpha) # Euler's constant
+sigma = pi/(sqrt(6.0)*alpha)
+xmin = mu - 6.0 * sigma
+xmax = mu + 6.0 * sigma
+ymax = 1.1 * extreme(u, u, alpha) #Mode of extreme PDF used
+ymax = int(10*ymax)/10.0
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100
+set title "extreme PDF"
+set key top left box
+plot [xmin:xmax] [0:ymax] extreme(x, 1.0, 0.5) title "alpha = 0.5, u = 1.0", \
+ extreme(x, 0.0, 1.0) title "alpha = 1.0, u = 0.0"
+pause -1 "Hit return to continue"
+set title "extreme CDF"
+plot [xmin:xmax] [0:1.1] cextreme(x, 1.0, 0.5) title "alpha = 0.5, u = 1.0", \
+ cextreme(x, 0.0, 1.0) title "alpha = 1.0, u = 0.0"
+pause -1 "Hit return to continue"
+
+# F PDF and CDF
+df1 = 5.0; df2 = 9.0
+mu = df2 < 2.0 ? 1.0 : df2 / (df2 - 2.0)
+sigma = df2 < 4.0 ? 1.0 : mu * sqrt(2.0 * (df1 + df2 - 2.0) / (df1 * (df2 - 4.0)))
+xmin = mu - 3.0 * sigma
+xmin = xmin < 0 ? 0 : xmin
+xmax = int(mu + 3.0 * sigma)
+#Mode of F PDF used
+ymax = df1 < 3.0 ? 1.0 : 1.1 * f((df1 / 2.0 - 1.0) / (df1 / 2.0 + df1 / df2), df1, df2)
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100
+set title "F PDF"
+set key right box
+plot [xmin:xmax] [0:ymax] f(x, 5.0, 9.0) title "df1 = 5, df2 = 9", \
+ f(x, 7.0, 6.0) title "df1 = 7, df2 = 6"
+pause -1 "Hit return to continue"
+set title "F CDF"
+set key left box
+plot [xmin:xmax] [0:1.1] cf(x, 5.0, 9.0) title "df1 = 5, df2 = 9", \
+ cf(x, 7.0, 6.0) title "df1 = 7, df2 = 6"
+pause -1 "Hit return to continue"
+
+# Gamma PDF and incomplete gamma CDF
+rho = 1.0; lambda = 1.3
+mu = rho / lambda
+sigma = sqrt(rho) / lambda
+xmin = mu - 4.0 * sigma
+xmin = xmin < 0 ? 0 : xmin
+xmax = mu + 4.0 * sigma
+ymax = rho < 1.0 ? 2.0 : 1.1 * gmm((rho - 1.0) / lambda, rho, lambda) #Mode of gamma pdf used
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set title "gamma PDF"
+set key right
+r1 = 0.5; r2 = 1.0; r3 = 1.0; r4 = 1.3; r5 = 2.0; r6 = 4.0; r7 = 6.0
+l1 = 1.0; l2 = 1.0; l3 = 1.3; l4 = 1.3; l5 = 2.0; l6 = 2.0; l7 = 2.0
+set arrow 1 from 1,1.3 to 0.15,gmm(0.15,r1,l1)
+set label 1 "rho < 1, tends to infinity" at 1.1,1.3 left
+set arrow 2 from 1.15,1.1 to 0.35,gmm(0.35,r3,l3)
+set label 2 "rho = 1, finite, nonzero limit" at 1.25,1.1 left
+set arrow 3 from 1.5,0.9 to 1.0,gmm(1.0,r5,l5)
+set label 3 "rho > 1, tends to zero" at 1.6,0.9 left
+keystr(rho,lambda) = sprintf("rho = %0.1f, lambda = %0.1f", rho, lambda)
+plot [0:5] [0:1.5] rho = r1, lambda = l1, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r2, lambda = l2, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r3, lambda = l3, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r4, lambda = l4, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r5, lambda = l5, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r6, lambda = l6, gmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r7, lambda = l7, gmm(x, rho, lambda) title keystr(rho,lambda)
+pause -1 "Hit return to continue"
+unset label 1; unset label 2; unset label 3
+unset arrow 1; unset arrow 2; unset arrow 3
+set title "incomplete gamma CDF"
+set key right bottom
+plot [0:5] [0:1.1] rho = r1, lambda = l1, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r2, lambda = l2, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r3, lambda = l3, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r4, lambda = l4, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r5, lambda = l5, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r6, lambda = l6, cgmm(x, rho, lambda) title keystr(rho,lambda), \
+ rho = r7, lambda = l7, cgmm(x, rho, lambda) title keystr(rho,lambda)
+pause -1 "Hit return to continue"
+
+# Geometric PDF and CDF
+p = 0.4
+mu = (1.0 - p) / p
+sigma = sqrt(mu / p)
+xmin = int(mu - 4.0 * sigma)
+xmin = xmin < -1 ? -1 : xmin
+xmin = -1
+xmax = int(mu + 4.0 * sigma)
+ymax = 1.1 * geometric(0, p) #mode of geometric PDF used
+unset key
+unset zeroaxis
+set xrange [xmin : xmax]
+set yrange [0 : ymax]
+set xlabel "k ->"
+set ylabel "probability density ->"
+set ytics 0, ymax / 10, ymax
+set format x "%2.0f"
+set format y "%3.2f"
+set sample (xmax - xmin) + 1
+set title "geometric PDF with p = 0.4"
+plot geometric(x, p) with impulses
+pause -1 "Hit return to continue"
+set title "geometric CDF with p = 0.4"
+set yrange [0 : 1.1]
+set ytics 0, 0.1, 1.0
+plot cgeometric(x, p) with steps
+pause -1 "Hit return to continue"
+
+# Half normal PDF and CDF
+mu = sqrt2invpi
+sigma = 1.0
+s = sigma*sqrt(1.0 - 2.0/pi)
+xmin = -0.2
+xmax = mu + 4.0 * s
+ymax = 1.1 * halfnormal(0, sigma) #Mode of half normal PDF used
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [-0.1: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set parametric
+set trange [xmin:xmax]
+set title "half normal PDF, sigma = 1.0"
+set arrow 1 from 0.5,0.13 to 0.0,0.4
+set label 1 "Discontinuity achieved by plotting\ntwice with limited parametric ranges" at 0.2,0.1 left
+plot t<0?t:-eps, halfnormal(t<0?t:-eps, sigma) ls 1, t<0?0.0:t, halfnormal(t<0?0.0:t, sigma) ls 1
+pause -1 "Hit return to continue"
+set title "half normal CDF, sigma = 1.0"
+set yrange [-0.1:1.1]
+set arrow 1 from 0.45,0.1 to 0.05,0.01
+set label 1 "Cusp achieved by plotting twice\nwith limited parametric ranges" at 0.5,0.1 left
+plot t<0?t:-eps, chalfnormal(t<0?t:-eps, sigma) ls 1, t<0?0.0:t, chalfnormal(t<0?0.0:t, sigma) ls 1
+pause -1 "Hit return to continue"
+unset label 1
+unset arrow 1
+unset parametric
+
+# Hypergeometric PDF and CPF
+N = 75; C = 25; d = 10
+p = real(C) / N
+mu = d * p
+sigma = sqrt(real(N - d) / (N - 1.0) * d * p * (1.0 - p))
+xmin = int(mu - 4.0 * sigma)
+xmin = xmin < -1 ? -1 : xmin
+xmax = int(mu + 4.0 * sigma)
+xmax = xmax < d+1 ? d+1 : xmax
+ymax = 1.1 * hypgeo(int(mu),N,C,d) # approximate mode of hypergeometric PDF used
+unset key
+unset zeroaxis
+set xrange [xmin : xmax]
+set yrange [0 : ymax]
+set xlabel "k ->"
+set ylabel "probability density ->"
+set ytics 0, ymax / 10, ymax
+set format x "%2.0f"
+set format y "%3.2f"
+set sample (xmax - xmin) + 1
+set title "hypergeometric PDF with N = 75, C = 25, d = 10"
+plot hypgeo(x,N,C,d) with impulses
+pause -1 "Hit return to continue"
+set yrange [0 : 1.1]
+set ytics 0, 1.0 / 10.0, 1.1
+set title "hypergeometric CDF with N = 75, C = 25, d = 10"
+plot chypgeo(x,N,C,d) with steps
+pause -1 "Hit return to continue"
+
+# Laplace PDF
+mu = 0.0; b = 1.0
+sigma = sqrt(2.0) * b
+xmin = mu - 4.0 * sigma
+xmax = mu + 4.0 * sigma
+ymax = 1.1 * laplace(mu, mu, b) #Mode of laplace PDF used
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [0: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100+1
+set title "laplace (or double exponential) PDF with mu = 0, b = 1"
+set arrow 1 from -0.95,0.5 to -0.1,0.5
+set label 1 "Cusp achieved by selecting point\nas part of function samples" at -1.0,0.5 right
+plot laplace(x, mu, b)
+pause -1 "Hit return to continue"
+unset label 1
+unset arrow 1
+set title "laplace (or double exponential) CDF with mu = 0, b = 1"
+set yrange [0: 1.1]
+plot claplace(x, mu, b)
+pause -1 "Hit return to continue"
+
+# Logistic PDF and CDF
+a = 0.0; lambda = 2.0
+mu = a
+sigma = pi / (sqrt(3.0) * lambda)
+xmin = mu - 4.0 * sigma
+xmax = mu + 4.0 * sigma
+ymax = 1.1 * logistic(mu, a, lambda) #Mode of logistic PDF used
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [0: ymax]
+unset key
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set title "logistic PDF with a = 0, lambda = 2"
+plot logistic(x, a, lambda)
+pause -1 "Hit return to continue"
+set title "logistic CDF with a = 0, lambda = 2"
+set yrange [0: 1.1]
+plot clogistic(x, a, lambda)
+pause -1 "Hit return to continue"
+
+# Lognormal PDF and CDF
+mu = 1.0; sigma = 0.5
+m = exp(mu + 0.5 * sigma**2)
+s = sqrt(exp(2.0 * mu + sigma**2) * (2.0 * exp(sigma) - 1.0))
+xmin = m - 4.0 * s
+xmin = xmin < 0 ? 0 : xmin
+xmax = m + 4.0 * s
+ymax = 1.1 * lognormal(exp(mu - sigma**2), mu, sigma) #Mode of lognormal PDF used
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [0: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.2f"
+set format y "%.2f"
+set sample 100
+set title "lognormal PDF with mu = 1.0, sigma = 0.5"
+plot lognormal(x, mu, sigma)
+pause -1 "Hit return to continue"
+set title "lognormal CDF with mu = 1.0, sigma = 0.5"
+set yrange [0: 1.1]
+plot clognormal(x, mu, sigma)
+pause -1 "Hit return to continue"
+
+# Maxwell PDF
+a = 0.5
+mu = 2.0 / sqrt(pi) / a
+sigma = sqrt(3.0 - 8.0/pi) / a
+xmin = int(mu - 3.0 * sigma)
+xmin = xmin < 0 ? 0 : xmin
+xmax = int(mu + 3.0 * sigma)
+a = 1.5
+ymax = 1.1 * maxwell(1.0 / a, a) + 0.5 #Mode of maxwell PDF used
+ymax = int(ymax + 0.5)
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set title "maxwell PDF"
+set key right top box
+plot [xmin:xmax] [0:ymax] maxwell(x, 1.5) title "a = 1.5", \
+ maxwell(x, 1.0) title "a = 1.0", \
+ maxwell(x, 0.5) title "a = 0.5"
+pause -1 "Hit return to continue"
+set title "maxwell CDF"
+set key right bottom box
+plot [xmin:xmax] [0:1.1] cmaxwell(x, 1.5) title "a = 1.5", \
+ cmaxwell(x, 1.0) title "a = 1.0", \
+ cmaxwell(x, 0.5) title "a = 0.5"
+pause -1 "Hit return to continue"
+
+# Negative binomial PDF and CDF
+r = 8; p = 0.4
+mu = r * (1.0 - p) / p
+sigma = sqrt(mu / p)
+xmin = int(mu - 4.0 * sigma)
+xmin = xmin < 0 ? 0 : xmin
+xmax = int(mu + 4.0 * sigma)
+ymax = 1.1 * negbin(int(mu - (1.0-p)/p), r, p) #mode of gamma PDF used
+unset key
+unset zeroaxis
+set xrange [xmin-1 : xmax]
+set yrange [0 : ymax]
+set xlabel "k ->"
+set ylabel "probability density ->"
+set ytics 0, ymax / 10, ymax
+set format x "%2.0f"
+set format y "%3.2f"
+set sample (xmax - xmin+1) + 1
+set title "negative binomial (or pascal or polya) PDF with r = 8, p = 0.4"
+plot negbin(x, r, p) with impulses
+pause -1 "Hit return to continue"
+set yrange [0 : 1.1]
+set ytics 0, 0.1, 1.0
+set title "negative binomial (or pascal or polya) CDF with r = 8, p = 0.4"
+plot cnegbin(x, r, p) with steps
+pause -1 "Hit return to continue"
+
+# Negative exponential PDF and CDF
+lambda = 2.0
+mu = 1.0 / lambda
+sigma = 1.0 / lambda
+xmax = mu + 4.0 * sigma
+ymax = lambda #No mode
+unset key
+set zeroaxis
+set xrange [0: xmax]
+set yrange [0: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.2f"
+set format y "%.1f"
+set sample 100
+set title "negative exponential (or exponential) PDF with lambda = 2.0"
+plot nexp(x, lambda)
+pause -1 "Hit return to continue"
+set title "negative exponential (or exponential) CDF with lambda = 2.0"
+set yrange [0: 1.1]
+plot cnexp(x, lambda)
+pause -1 "Hit return to continue"
+
+# Normal PDF and CDF
+mu = 0.0; sigma = 1.0
+xmin = mu - 4.0 * sigma
+xmax = mu + 4.0 * sigma
+mu = 2.0; sigma = 0.5
+ymax = 1.1 * normal(mu, mu, sigma) #Mode of normal PDF used
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 100
+set title "normal (also called gauss or bell-curved) PDF"
+set key left top box
+plot [xmin:xmax] [0:ymax] normal(x, 0, 1.0) title "mu = 0, sigma = 1.0", \
+ normal(x, 2, 0.5) title "mu = 2, sigma = 0.5", \
+ normal(x, 1, 2.0) title "mu = 1, sigma = 2.0"
+pause -1 "Hit return to continue"
+set title "normal (also called gauss or bell-curved) CDF"
+set key left top box
+plot [xmin:xmax] [0:1.1] mu = 0, sigma = 1.0, cnormal(x, mu, sigma) title "mu = 0, sigma = 1.0", \
+ mu = 2, sigma = 0.5, cnormal(x, mu, sigma) title "mu = 2, sigma = 0.5", \
+ mu = 1, sigma = 2.0, cnormal(x, mu, sigma) title "mu = 1, sigma = 2.0"
+pause -1 "Hit return to continue"
+
+# Pareto PDF and CDF
+a = 1.0; b = 3.0
+mu = a * b / (b - 1.0)
+sigma = a * sqrt(b) / (sqrt(b - 2.0) * (b - 1.0))
+xmin = mu - 4.0 * sigma
+xmin = xmin < 0 ? 0 : xmin
+xmax = int(mu + 4.0 * sigma)
+ymax = 1.1 * pareto(a, a, b) #mode of pareto PDF used
+ymin = -0.1 * pareto(a, a, b)
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [ymin: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.1f"
+set sample 200+1
+set title "pareto PDF with a = 1, b = 3"
+# Discontinuity at a
+set parametric
+set trange [0:1-eps]
+x1(t) = -1 + 2*t
+x2(t) = 1 + 3*t
+set arrow 1 from 1.75,0.8 to 1.0,0.8
+set arrow 2 from 1.0,0.0 to 1.0,3.0 nohead lt 0
+set label 1 "Discontinuity achieved by plotting twice\nwith affine mapped parametric ranges" at 1.8,0.8 left
+plot x1(t), pareto(x1(t), a, b) ls 1, x2(t), pareto(x2(t), a, b) ls 1
+pause -1 "Hit return to continue"
+unset arrow 2
+set title "pareto CDF with a = 1, b = 3"
+unset parametric
+set yrange [-0.1: 1.1]
+set arrow 1 from 1.45,0.1 to 1.05,0.01
+set label 1 "Cusp achieved by selecting point\nas part of function samples" at 1.5,0.1 left
+plot cpareto(x, a, b)
+pause -1 "Hit return to continue"
+unset label 1
+unset arrow 1
+
+# Poisson PDF and CDF
+mu = 4.0
+sigma = sqrt(mu)
+xmin = int(mu - 4.0 * sigma)
+xmin = xmin < -1 ? -1 : xmin
+xmax = int(mu + 4.0 * sigma)
+ymax = 1.1 * poisson(mu, mu) #mode of poisson PDF used
+unset key
+set zeroaxis
+set xrange [xmin : xmax]
+set yrange [0 : ymax]
+set xlabel "k ->"
+set ylabel "probability density ->"
+set ytics 0, ymax / 10, ymax
+set format x "%2.0f"
+set format y "%3.2f"
+set sample (xmax - xmin) + 1
+set title "poisson PDF with mu = 4.0"
+plot poisson(x, mu) with impulses
+pause -1 "Hit return to continue"
+set yrange [-0.1 : 1.1]
+set ytics -0.1, 0.1, 1.1
+set title "poisson CDF with mu = 4.0"
+plot cpoisson(x, mu) with steps
+pause -1 "Hit return to continue"
+
+# Rayleigh PDF and CDF
+lambda = 2.0
+mu = 0.5 * sqrt(pi / lambda)
+sigma = sqrt((1.0 - pi / 4.0) / lambda)
+xmax = mu + 4.0 * sigma
+ymax = 1.1 * rayleigh(1.0 / sqrt(2.0 * lambda), lambda) #Mode of rayleigh PDF used
+unset key
+set zeroaxis
+set xrange [0: xmax]
+set yrange [0: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.2f"
+set format y "%.1f"
+set sample 100
+set title "rayleigh PDF with lambda = 2.0"
+plot rayleigh(x, lambda)
+pause -1 "Hit return to continue"
+set title "rayleigh CDF with lambda = 2.0"
+set yrange [0: 1.1]
+plot crayleigh(x, lambda)
+pause -1 "Hit return to continue"
+
+# Sine PDF and CDF
+a = 3.2; f = 2.6
+mu = a / 2.0
+sigma = sqrt(a * a / 3.0 * (1.0 - 3.0 / (2.0 * n * n * pi * pi)) - mu * mu)
+xmin = 0.0
+xmax = a - eps
+a = 2; f = 1.0
+ymax = 1.1 * 2.0 / a #Mode of sine PDF used
+set zeroaxis
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.2f"
+set format y "%.1f"
+set sample 250
+set title "sine PDF"
+set key bottom outside
+keystr(a, f) = sprintf("a = %0.1f, f = %0.1f", a, f)
+a1 = 2.0; a2 = 3.25; a3 = 2.75
+f1 = 1.0; f2 = 3.0; f3 = 2.6; f4 = 0.0
+plot [xmin:xmax] [0:ymax] a = a1, f = f1, sine(x, f, a) title keystr(a, f), \
+ a = a1, f = f2, sine(x, f, a) title keystr(a, f), \
+ a = a2, f = f3, sine(x, f, a) title keystr(a, f), \
+ a = a3, f = f4, sine(x, f, a) title keystr(a, f) with steps
+pause -1 "Hit return to continue"
+set title "sine CDF"
+set key top left
+plot [xmin:xmax] [0:1.1] a = a1, f = f1, csine(x, f, a) title keystr(a, f), \
+ a = a1, f = f2, csine(x, f, a) title keystr(a, f), \
+ a = a2, f = f3, csine(x, f, a) title keystr(a, f), \
+ a = a3, f = f4, csine(x, f, a) title keystr(a, f) with steps
+pause -1 "Hit return to continue"
+
+# t PDF and CDF
+nu = 20
+mu = 0
+sigma = nu > 2 ? sqrt(nu / (nu - 2.0)) : 1.0
+xmin = mu - 4.0 * sigma
+xmax = mu + 4.0 * sigma
+ymax = 1.1 * t(mu, nu) #Mode of t PDF used
+set key inside center left title "degrees of freedom"
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [0: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 100
+set title "t PDF (and Gaussian limit)"
+ks(nu) = sprintf("nu = %d", nu)
+plot t(x, 1) ti ks(1), t(x, 2) ti ks(2), t(x, 4) ti ks(4), t(x, 10) ti ks(10), \
+ t(x, 20) ti ks(20), normal(x, 0, 1) ti "normal"
+pause -1 "Hit return to continue"
+set title "t CDF (and Gaussian limit)"
+set yrange [0: 1.1]
+plot ct(x, 1) ti ks(1), ct(x, 2) ti ks(2), ct(x, 4) ti ks(4), ct(x, 10) ti ks(10), \
+ ct(x, 20) ti ks(20), cnormal(x, 0, 1) ti "normal"
+pause -1 "Hit return to continue"
+
+# Thanks to efrank@upenn5.hep.upenn.edu for telling us about this
+# triangular PDF and CDF
+m = 3.0
+g = 2.0
+mu = m
+sigma = g/sqrt(6.0)
+xmin = m - 1.1*g
+xmax = m + 1.1*g
+ymax = 1.1 * triangular(m, m, g) #Mode of triangular PDF used
+ymin = -ymax/11.0;
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [ymin: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.1f"
+set format y "%.2f"
+set sample 50*1.1*g+1
+set title "triangular PDF with m = 3.0, g = 2.0"
+plot triangular(x, m, g)
+pause -1 "Hit return to continue"
+set title "triangular CDF with m = 3.0, g = 2.0"
+set yrange [-0.1: 1.1]
+plot ctriangular(x, m, g)
+pause -1 "Hit return to continue"
+
+# Uniform PDF and CDF
+a = -2.0; b= 2.0
+mu = (a + b) / 2.0
+sigma = (b - a) / sqrt(12.0)
+xmin = a - 0.1*(b - a)
+xmax = b + 0.1*(b - a)
+ymax = 1.1 * uniform(mu, a, b) #No mode
+ymin = -0.1 * uniform(mu, a, b)
+unset key
+set zeroaxis
+set xrange [xmin: xmax]
+set yrange [ymin: ymax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%.2f"
+set format y "%.2f"
+set sample 120+1
+set title "uniform PDF with a = -2.0, b = 2.0"
+plot uniform(x, a, b) with steps
+pause -1 "Hit return to continue"
+set title "uniform CDF with a = -2.0, b = 2.0"
+set yrange [-0.1 : 1.1]
+plot cuniform(x, a, b)
+pause -1 "Hit return to continue"
+
+# Weibull PDF and CDF
+lambda = 1.0/5; a = 1.0
+mu = 1.0 / lambda * gamma(1.0 / a) / a
+sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2))
+xmin = mu - 4.0 * sigma
+xmin = xmin < 0 ? 0 : xmin
+#Mode of weibull PDF used
+ymax = 1.8 * (a >= 1.0 ? weibull(((a - 1.0) / a)**(1.0 / a) / lambda, a, lambda) : 2.0)
+lambda = 1.0/15; a = 10.0
+mu = 1.0 / lambda * gamma(1.0 / a) / a
+sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2))
+xmax = int(mu + 4.0 * sigma)
+set key on title "" inside top right
+set zeroaxis
+set grid
+set xrange [xmin : xmax]
+set xlabel "x ->"
+set ylabel "probability density ->"
+set xtics autofreq
+set ytics autofreq
+set format x "%g"
+set format y "%g"
+set sample 100
+set title "weibull PDF"
+ks(a,lambda) = sprintf("lambda = 1/%g, a = %0.1f", 1.0/lambda, a)
+a1 = 0.5; a2 = 1.0; a3 = 2.0; a4 = 10.0
+lambda1 = 1.0/5; lambda2 = 1.0/15
+set arrow 1 from 3.8,0.27 to 0.5,weibull(0.5,a1,lambda1)
+set label 1 "a < 1, rate descreasing over time" at 4,0.27 left
+set arrow 2 from 8,0.19 to 6.4,weibull(6.4,a3,lambda1)
+set arrow 3 from 10.5,0.19 to 13,weibull(13,a4,lambda2)
+set label 2 "a > 1, rate increasing over time" at 9,0.2 center
+plot [] [0:ymax] lambda = lambda1, a = a1, weibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda1, a = a2, weibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda1, a = a3, weibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda2, a = a4, weibull(x, a, lambda) ti ks(a, lambda)
+pause -1 "Hit return to continue"
+unset label 1; unset label 2
+unset arrow 1; unset arrow 2; unset arrow 3
+set key at 9,0.4 center
+set title "weibull CDF"
+plot [] [0:1.1] lambda = lambda1, a = a1, cweibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda1, a = a2, cweibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda1, a = a3, cweibull(x, a, lambda) ti ks(a, lambda), \
+ lambda = lambda2, a = a4, cweibull(x, a, lambda) ti ks(a, lambda)
+pause -1 "Hit return to continue"
+reset
projects / gnuplot / blobdiff