vcs.maemo.org Git - gnuplot/blob - demo/prob.dem

   1 #
   2 # $Id: prob.dem,v 1.9 2006/06/14 03:24:09 sfeam Exp $
   3 #
   4 # Demo Statistical Functions version 2.3
   5 #
   6 # Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl
   7
   8 # History:
   9 #    -- --- 1992 Jos van der Woude:  1st version
  10 #    06 Jun 2006 Dan Sebald:  Added some variety and plotting techniques for
  11 #                             better visual effect.  More tutorial in nature.
  12
  13 print "                   Statistical Library Demo, version 2.3"
  14 print "\n          Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl"
  15 print "\n\n\n\n\n\n\n"
  16 print "NOTE: contains 54 plots and consequently takes a lot of time to run"
  17 print "                      Press Ctrl-C to exit right now"
  18 pause -1 "                      Press Return to start demo ..."
  19
  20 load "stat.inc"
  21
  22 eps = 1.0e-10  # Supposed to be float resolution (nice if were defined internally)
  23
  24 ## Gamma function
  25 xmin = -5.5
  26 xmax = 5
  27 ymin = -10
  28 ymax = 10
  29 unset key
  30 set xzeroaxis
  31 gsampfunc(t,n) = t<0?0.5*1/(-t+1.0)**n:1.0-0.5*1/(t+1.0)**n
  32 set parametric
  33 set trange [-1:1]
  34 set sample 200
  35 set xrange [xmin : xmax]
  36 set yrange [ymin : ymax]
  37 set xlabel "x"
  38 set ylabel "gamma(x)"
  39 set arrow 1 from  0,ymin to  0,ymax nohead lt 0
  40 set arrow 2 from -1,ymin to -1,ymax nohead lt 0
  41 set arrow 3 from -2,ymin to -2,ymax nohead lt 0
  42 set arrow 4 from -3,ymin to -3,ymax nohead lt 0
  43 set arrow 5 from -4,ymin to -4,ymax nohead lt 0
  44 set arrow 6 from -5,ymin to -5,ymax nohead lt 0
  45 set title "gamma function, very useful function for probability"
  46 plot gsampfunc(5*t,5)-6, gamma(gsampfunc(5*t,5)-6) lt 1, \
  47      gsampfunc(5*t,5)-5, gamma(gsampfunc(5*t,5)-5) lt 1, \
  48      gsampfunc(5*t,4)-4, gamma(gsampfunc(5*t,4)-4) lt 1, \
  49      gsampfunc(5*t,3)-3, gamma(gsampfunc(5*t,3)-3) lt 1, \
  50      gsampfunc(5*t,2)-2, gamma(gsampfunc(5*t,2)-2) lt 1, \
  51      gsampfunc(5*t,1)-1, gamma(gsampfunc(5*t,1)-1) lt 1, \
  52      5*gsampfunc(5*t,2), gamma(5*gsampfunc(5*t,2)) lt 1
  53 pause -1 "Hit return to continue"
  54 ymin = ymin/2
  55 ymax = ymax/2
  56 set yrange [ymin:ymax]
  57 set ylabel "lgamma(x)"
  58 set arrow 1 from  0,ymin to  0,ymax nohead lt 0
  59 set arrow 2 from -1,ymin to -1,ymax nohead lt 0
  60 set arrow 3 from -2,ymin to -2,ymax nohead lt 0
  61 set arrow 4 from -3,ymin to -3,ymax nohead lt 0
  62 set arrow 5 from -4,ymin to -4,ymax nohead lt 0
  63 set arrow 6 from -5,ymin to -5,ymax nohead lt 0
  64 set title "log gamma function, similarly very useful function"
  65 plot gsampfunc(5*t,5)-6, lgamma(gsampfunc(5*t,5)-6) lt 1, \
  66      gsampfunc(5*t,5)-5, lgamma(gsampfunc(5*t,5)-5) lt 1, \
  67      gsampfunc(5*t,4)-4, lgamma(gsampfunc(5*t,4)-4) lt 1, \
  68      gsampfunc(5*t,3)-3, lgamma(gsampfunc(5*t,3)-3) lt 1, \
  69      gsampfunc(5*t,3)-2, lgamma(gsampfunc(5*t,3)-2) lt 1, \
  70      gsampfunc(5*t,3)-1, lgamma(gsampfunc(5*t,3)-1) lt 1, \
  71      5*gsampfunc(5*t,3), lgamma(5*gsampfunc(5*t,3)) lt 1
  72 pause -1 "Hit return to continue"
  73 reset
  74
  75 # Arcsinus PDF and CDF
  76 r = 2.0
  77 mu = 0.0
  78 sigma = r / sqrt2
  79 xmin = -(r+1)
  80 xmax = r+1
  81 unset key
  82 set zeroaxis
  83 set xrange [xmin : xmax]
  84 set yrange [-0.2 : 1.2]
  85 set xlabel "x ->"
  86 set ylabel "probability density ->"
  87 set xtics autofreq
  88 set ytics autofreq
  89 set format x "%.1f"
  90 set format y "%.1f"
  91 set sample 50*xmax+1
  92 set title "arcsin PDF with r = 2.0"
  93 plot arcsin(x, r)
  94 pause -1 "Hit return to continue"
  95 set title "arcsin CDF with r = 2.0"
  96 set yrange [-0.2 : 1.2]
  97 plot carcsin(x, r)
  98 pause -1 "Hit return to continue"
  99
 100 # Beta PDF and CDF
 101 p = 5.0; q = 3.0
 102 mu = p / (p + q)
 103 sigma = sqrt(p**q) / ((p + q ) * sqrt(p + q + 1.0))
 104 xmin = 0.0
 105 xmax = 1.0
 106 #Mode of beta PDF used
 107 ymax = (p < 1.0 || q < 1.0) ? 2.0 : 1.4 * beta((p - 1.0)/(p + q - 2.0), p, q)
 108 set key right box
 109 set zeroaxis
 110 set xrange [xmin : xmax]
 111 set yrange [0 : ymax]
 112 set xlabel "x ->"
 113 set ylabel "probability density ->"
 114 set xtics autofreq
 115 set ytics autofreq
 116 set format x "%.1f"
 117 set format y "%.1f"
 118 set sample 200
 119 set title "beta PDF"
 120 plot beta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \
 121      beta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \
 122      beta(x, 0.5, 2.5) title "p = 0.5, q = 2.5"
 123 pause -1 "Hit return to continue"
 124 set yrange [0:1.1]
 125 set title "incomplete beta CDF"
 126 set key left box
 127 plot cbeta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \
 128      cbeta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \
 129      cbeta(x, 0.5, 2.5) title "p = 0.5, q = 2.5"
 130 pause -1 "Hit return to continue"
 131
 132 # Binomial PDF and CDF
 133 n = 25; p = 0.15
 134 mu = n * p
 135 sigma = sqrt(n * p * (1.0 - p))
 136 xmin = int(mu - 4.0 * sigma)
 137 xmin = xmin < -2 ? -2 : xmin
 138 xmax = int(mu + 4.0 * sigma)
 139 xmax = xmax < n+2 ? n+2 : xmax
 140 ymax = 1.1 * binom(int((n+1)*p), n, p) #Mode of binomial PDF used
 141 unset key
 142 unset zeroaxis
 143 set xrange [xmin : xmax]
 144 set yrange [0 : ymax]
 145 set xlabel "k ->"
 146 set ylabel "probability density ->"
 147 set ytics 0, ymax / 10, ymax
 148 set format x "%2.0f"
 149 set format y "%3.2f"
 150 set sample (xmax - xmin) + 1
 151 set title "binomial PDF with n = 25, p = 0.15"
 152 plot binom(x, n, p) with impulses
 153 pause -1 "Hit return to continue"
 154 set ytics autofreq
 155 set xzeroaxis
 156 set title "binomial CDF with n = 25, p = 0.15"
 157 set yrange [-0.1 : 1.1]
 158 set ytics 0, 0.1, 1.0
 159 plot cbinom(x, n, p) with steps
 160 pause -1 "Hit return to continue"
 161
 162 # Cauchy PDF and CDF
 163 a = 0.0; b = 2.0
 164 #cauchy PDF has no moments
 165 xmin = a - 5.0 * b
 166 xmax = a + 5.0 * b
 167 ymax = 1.1 * cauchy(a, a, b) #Mode of cauchy PDF used
 168 set key left box
 169 set zeroaxis
 170 set xlabel "x ->"
 171 set ylabel "probability density ->"
 172 set xtics autofreq
 173 set ytics autofreq
 174 set format x "%.1f"
 175 set format y "%.2f"
 176 set sample 100
 177 set title "cauchy PDF"
 178 a=0
 179 b=2
 180 plot [xmin:xmax] [0:ymax] cauchy(x, 0, 2) title "a = 0, b = 2", \
 181                           cauchy(x, 0, 4) title "a = 0, b = 4"
 182 pause -1 "Hit return to continue"
 183 set title "cauchy CDF"
 184 plot [xmin:xmax] [0:1.0] ccauchy(x, 0, 2) title "a = 0, b = 2", \
 185                          ccauchy(x, 0, 4) title "a = 0, b = 4"
 186 pause -1 "Hit return to continue"
 187
 188 # Chi-square PDF and CDF
 189 k = 4.0
 190 mu = k
 191 sigma = sqrt(2.0 * k)
 192 xmin = mu - 4.0 * sigma
 193 xmin = xmin < 0 ? 0 : xmin
 194 xmax = int(mu + 4.0 * sigma)
 195 k = 2.0
 196 ymax = (k > 2.0 ? 1.1*chisq(k - 2.0, k) : 0.5) #Mode of chi PDF used
 197 set key right box
 198 set zeroaxis
 199 set xrange [xmin+eps : xmax] #Discontinuity at zero for k < 2
 200 set yrange [0:ymax]
 201 set xlabel "x ->"
 202 set ylabel "probability density ->"
 203 set xtics autofreq
 204 set ytics autofreq
 205 set format x "%.1f"
 206 set format y "%.2f"
 207 set sample 100
 208 set title "chi-square PDF"
 209 set key right box
 210 set samples 15*20+1
 211 keystr(k) = sprintf("k = %d", k)
 212 plot k = 1, x==0?1/0:chisq(x, k) title keystr(k), \
 213      k = 2, x==0?1/0:chisq(x, k) title keystr(k), \
 214      k = 3, chisq(x, k) title keystr(k), \
 215      k = 4, chisq(x, k) title keystr(k), \
 216      k = 5, chisq(x, k) title keystr(k), \
 217      k = 6, chisq(x, k) title keystr(k), \
 218      k = 7, chisq(x, k) title keystr(k), \
 219      k = 8, chisq(x, k) title keystr(k)
 220 pause -1 "Hit return to continue"
 221 set yrange [0:1.1]
 222 set key bottom right box
 223 set title "chi-square CDF"
 224 plot k = 1, cchisq(x, k) title keystr(k), \
 225      k = 2, cchisq(x, k) title keystr(k), \
 226      k = 3, cchisq(x, k) title keystr(k), \
 227      k = 4, cchisq(x, k) title keystr(k), \
 228      k = 5, cchisq(x, k) title keystr(k), \
 229      k = 6, cchisq(x, k) title keystr(k), \
 230      k = 7, cchisq(x, k) title keystr(k), \
 231      k = 8, cchisq(x, k) title keystr(k)
 232 pause -1 "Hit return to continue"
 233
 234 # Erlang PDF and CDF
 235 lambda = 1.0; n = 2.0
 236 mu = n / lambda
 237 sigma = sqrt(n) / lambda
 238 xmax = int(mu + 5.0 * sigma)
 239 n = 1.0
 240 ymax = n < 2.0 ? 1.0 : 1.1 * erlang((n - 1.0) / lambda, n, lambda) #Mode of erlang PDF used
 241 set zeroaxis
 242 set xlabel "x ->"
 243 set ylabel "probability density ->"
 244 set xtics autofreq
 245 set ytics autofreq
 246 set format x "%.1f"
 247 set format y "%.1f"
 248 set sample 100
 249 set title "erlang PDF"
 250 set key top right box
 251 l1 = 1.0; l2 = 0.5
 252 set arrow 1 from 2,0.8 to 0.33,erlang(0.33,1,l1)
 253 set arrow 2 from 2,0.8 to 0.33,erlang(0.33,1,l2)
 254 set label 1 "n = 1, exponential r.v." at 2.1,0.8 left
 255 keystr(n,lambda) = sprintf("lambda = %0.1f, n = %d", lambda, n)
 256 plot [0:xmax] [0:ymax] n = 1, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \
 257                        n = 1, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda), \
 258                        n = 2, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \
 259                        n = 2, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda)
 260 pause -1 "Hit return to continue"
 261 unset label 1
 262 unset arrow 1; unset arrow 2
 263 set title "erlang CDF"
 264 set key bottom right box
 265 plot [0:xmax] [0:1.1] n = 1, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \
 266                       n = 1, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda), \
 267                       n = 2, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \
 268                       n = 2, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda)
 269 pause -1 "Hit return to continue"
 270
 271 # Thanks to mrb2j@kelvin.seas.Virginia.EDU for telling us about this.
 272 # Extreme (Gumbel extreme value) PDF and CDF
 273 alpha = 1.0; u = 0.0
 274 mu = u + (0.577215665/alpha)   # Euler's constant
 275 sigma = pi/(sqrt(6.0)*alpha)
 276 xmin = mu - 6.0 * sigma
 277 xmax = mu + 6.0 * sigma
 278 ymax = 1.1 * extreme(u, u, alpha) #Mode of extreme PDF used
 279 ymax = int(10*ymax)/10.0
 280 set zeroaxis
 281 set xlabel "x ->"
 282 set ylabel "probability density ->"
 283 set xtics autofreq
 284 set ytics autofreq
 285 set format x "%.1f"
 286 set format y "%.2f"
 287 set sample 100
 288 set title "extreme PDF"
 289 set key top left box
 290 plot [xmin:xmax] [0:ymax] extreme(x, 1.0, 0.5) title "alpha = 0.5, u = 1.0", \
 291                           extreme(x, 0.0, 1.0) title "alpha = 1.0, u = 0.0"
 292 pause -1 "Hit return to continue"
 293 set title "extreme CDF"
 294 plot [xmin:xmax] [0:1.1] cextreme(x, 1.0, 0.5) title "alpha = 0.5, u = 1.0", \
 295                          cextreme(x, 0.0, 1.0) title "alpha = 1.0, u = 0.0"
 296 pause -1 "Hit return to continue"
 297
 298 # F PDF and CDF
 299 df1 = 5.0; df2 = 9.0
 300 mu = df2 < 2.0 ? 1.0 : df2 / (df2 - 2.0)
 301 sigma = df2 < 4.0 ? 1.0 : mu * sqrt(2.0 * (df1 + df2 - 2.0) / (df1 * (df2 - 4.0)))
 302 xmin = mu - 3.0 * sigma
 303 xmin = xmin < 0 ? 0 : xmin
 304 xmax = int(mu + 3.0 * sigma)
 305 #Mode of F PDF used
 306 ymax = df1 < 3.0 ? 1.0 : 1.1 * f((df1 / 2.0 - 1.0) / (df1 / 2.0 + df1 / df2), df1, df2)
 307 set zeroaxis
 308 set xlabel "x ->"
 309 set ylabel "probability density ->"
 310 set xtics autofreq
 311 set ytics autofreq
 312 set format x "%.1f"
 313 set format y "%.2f"
 314 set sample 100
 315 set title "F PDF"
 316 set key right box
 317 plot [xmin:xmax] [0:ymax] f(x, 5.0, 9.0) title "df1 = 5, df2 = 9", \
 318                           f(x, 7.0, 6.0) title "df1 = 7, df2 = 6"
 319 pause -1 "Hit return to continue"
 320 set title "F CDF"
 321 set key left box
 322 plot [xmin:xmax] [0:1.1] cf(x, 5.0, 9.0) title "df1 = 5, df2 = 9", \
 323                          cf(x, 7.0, 6.0) title "df1 = 7, df2 = 6"
 324 pause -1 "Hit return to continue"
 325
 326 # Gamma PDF and incomplete gamma CDF
 327 rho = 1.0; lambda = 1.3
 328 mu = rho / lambda
 329 sigma = sqrt(rho) / lambda
 330 xmin = mu - 4.0 * sigma
 331 xmin = xmin < 0 ? 0 : xmin
 332 xmax = mu + 4.0 * sigma
 333 ymax = rho < 1.0 ? 2.0 : 1.1 * gmm((rho - 1.0) / lambda, rho, lambda) #Mode of gamma pdf used
 334 set zeroaxis
 335 set xlabel "x ->"
 336 set ylabel "probability density ->"
 337 set xtics autofreq
 338 set ytics autofreq
 339 set format x "%.1f"
 340 set format y "%.1f"
 341 set sample 100
 342 set title "gamma PDF"
 343 set key right
 344 r1 = 0.5; r2 = 1.0; r3 = 1.0; r4 = 1.3; r5 = 2.0; r6 = 4.0; r7 = 6.0
 345 l1 = 1.0; l2 = 1.0; l3 = 1.3; l4 = 1.3; l5 = 2.0; l6 = 2.0; l7 = 2.0
 346 set arrow 1 from 1,1.3 to 0.15,gmm(0.15,r1,l1)
 347 set label 1 "rho < 1, tends to infinity" at 1.1,1.3 left
 348 set arrow 2 from 1.15,1.1 to 0.35,gmm(0.35,r3,l3)
 349 set label 2 "rho = 1, finite, nonzero limit" at 1.25,1.1 left
 350 set arrow 3 from 1.5,0.9 to 1.0,gmm(1.0,r5,l5)
 351 set label 3 "rho > 1, tends to zero" at 1.6,0.9 left
 352 keystr(rho,lambda) = sprintf("rho = %0.1f, lambda = %0.1f", rho, lambda)
 353 plot [0:5] [0:1.5] rho = r1, lambda = l1, gmm(x, rho, lambda) title keystr(rho,lambda), \
 354                    rho = r2, lambda = l2, gmm(x, rho, lambda) title keystr(rho,lambda), \
 355                    rho = r3, lambda = l3, gmm(x, rho, lambda) title keystr(rho,lambda), \
 356                    rho = r4, lambda = l4, gmm(x, rho, lambda) title keystr(rho,lambda), \
 357                    rho = r5, lambda = l5, gmm(x, rho, lambda) title keystr(rho,lambda), \
 358                    rho = r6, lambda = l6, gmm(x, rho, lambda) title keystr(rho,lambda), \
 359                    rho = r7, lambda = l7, gmm(x, rho, lambda) title keystr(rho,lambda)
 360 pause -1 "Hit return to continue"
 361 unset label 1; unset label 2; unset label 3
 362 unset arrow 1; unset arrow 2; unset arrow 3
 363 set title "incomplete gamma CDF"
 364 set key right bottom
 365 plot [0:5] [0:1.1] rho = r1, lambda = l1, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 366                    rho = r2, lambda = l2, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 367                    rho = r3, lambda = l3, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 368                    rho = r4, lambda = l4, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 369                    rho = r5, lambda = l5, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 370                    rho = r6, lambda = l6, cgmm(x, rho, lambda) title keystr(rho,lambda), \
 371                    rho = r7, lambda = l7, cgmm(x, rho, lambda) title keystr(rho,lambda)
 372 pause -1 "Hit return to continue"
 373
 374 # Geometric PDF and CDF
 375 p = 0.4
 376 mu = (1.0 - p) / p
 377 sigma = sqrt(mu / p)
 378 xmin = int(mu - 4.0 * sigma)
 379 xmin = xmin < -1 ? -1 : xmin
 380 xmin = -1
 381 xmax = int(mu + 4.0 * sigma)
 382 ymax = 1.1 * geometric(0, p) #mode of geometric PDF used
 383 unset key
 384 unset zeroaxis
 385 set xrange [xmin : xmax]
 386 set yrange [0 : ymax]
 387 set xlabel "k ->"
 388 set ylabel "probability density ->"
 389 set ytics 0, ymax / 10, ymax
 390 set format x "%2.0f"
 391 set format y "%3.2f"
 392 set sample (xmax - xmin) + 1
 393 set title "geometric PDF with p = 0.4"
 394 plot geometric(x, p) with impulses
 395 pause -1 "Hit return to continue"
 396 set title "geometric CDF with p = 0.4"
 397 set yrange [0 : 1.1]
 398 set ytics 0, 0.1, 1.0
 399 plot cgeometric(x, p) with steps
 400 pause -1 "Hit return to continue"
 401
 402 # Half normal PDF and CDF
 403 mu = sqrt2invpi
 404 sigma = 1.0
 405 s = sigma*sqrt(1.0 - 2.0/pi)
 406 xmin = -0.2
 407 xmax = mu + 4.0 * s
 408 ymax = 1.1 * halfnormal(0, sigma) #Mode of half normal PDF used
 409 unset key
 410 set zeroaxis
 411 set xrange [xmin: xmax]
 412 set yrange [-0.1: ymax]
 413 set xlabel "x ->"
 414 set ylabel "probability density ->"
 415 set xtics autofreq
 416 set ytics autofreq
 417 set format x "%.1f"
 418 set format y "%.1f"
 419 set sample 100
 420 set parametric
 421 set trange [xmin:xmax]
 422 set title "half normal PDF, sigma = 1.0"
 423 set arrow 1 from 0.5,0.13 to 0.0,0.4
 424 set label 1 "Discontinuity achieved by plotting\ntwice with limited parametric ranges" at 0.2,0.1 left
 425 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
 426 pause -1 "Hit return to continue"
 427 set title "half normal CDF, sigma = 1.0"
 428 set yrange [-0.1:1.1]
 429 set arrow 1 from 0.45,0.1 to 0.05,0.01
 430 set label 1 "Cusp achieved by plotting twice\nwith limited parametric ranges" at 0.5,0.1 left
 431 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
 432 pause -1 "Hit return to continue"
 433 unset label 1
 434 unset arrow 1
 435 unset parametric
 436
 437 # Hypergeometric PDF and CPF
 438 N = 75; C = 25; d = 10
 439 p = real(C) / N
 440 mu = d * p
 441 sigma = sqrt(real(N - d) / (N - 1.0) * d * p * (1.0 - p))
 442 xmin = int(mu - 4.0 * sigma)
 443 xmin = xmin < -1 ? -1 : xmin
 444 xmax = int(mu + 4.0 * sigma)
 445 xmax = xmax < d+1 ? d+1 : xmax
 446 ymax = 1.1 * hypgeo(int(mu),N,C,d) # approximate mode of hypergeometric PDF used
 447 unset key
 448 unset zeroaxis
 449 set xrange [xmin : xmax]
 450 set yrange [0 : ymax]
 451 set xlabel "k ->"
 452 set ylabel "probability density ->"
 453 set ytics 0, ymax / 10, ymax
 454 set format x "%2.0f"
 455 set format y "%3.2f"
 456 set sample (xmax - xmin) + 1
 457 set title "hypergeometric PDF with N = 75, C = 25, d = 10"
 458 plot hypgeo(x,N,C,d) with impulses
 459 pause -1 "Hit return to continue"
 460 set yrange [0 : 1.1]
 461 set ytics 0, 1.0 / 10.0, 1.1
 462 set title "hypergeometric CDF with N = 75, C = 25, d = 10"
 463 plot chypgeo(x,N,C,d) with steps
 464 pause -1 "Hit return to continue"
 465
 466 # Laplace PDF
 467 mu = 0.0; b = 1.0
 468 sigma = sqrt(2.0) * b
 469 xmin = mu - 4.0 * sigma
 470 xmax = mu + 4.0 * sigma
 471 ymax = 1.1 * laplace(mu, mu, b) #Mode of laplace PDF used
 472 unset key
 473 set zeroaxis
 474 set xrange [xmin: xmax]
 475 set yrange [0: ymax]
 476 set xlabel "x ->"
 477 set ylabel "probability density ->"
 478 set xtics autofreq
 479 set ytics autofreq
 480 set format x "%.1f"
 481 set format y "%.2f"
 482 set sample 100+1
 483 set title "laplace (or double exponential) PDF with mu = 0, b = 1"
 484 set arrow 1 from -0.95,0.5 to -0.1,0.5
 485 set label 1 "Cusp achieved by selecting point\nas part of function samples" at -1.0,0.5 right
 486 plot laplace(x, mu, b)
 487 pause -1 "Hit return to continue"
 488 unset label 1
 489 unset arrow 1
 490 set title "laplace (or double exponential) CDF with mu = 0, b = 1"
 491 set yrange [0: 1.1]
 492 plot claplace(x, mu, b)
 493 pause -1 "Hit return to continue"
 494
 495 # Logistic PDF and CDF
 496 a = 0.0; lambda = 2.0
 497 mu = a
 498 sigma = pi / (sqrt(3.0) * lambda)
 499 xmin = mu - 4.0 * sigma
 500 xmax = mu + 4.0 * sigma
 501 ymax = 1.1 * logistic(mu, a, lambda) #Mode of logistic PDF used
 502 unset key
 503 set zeroaxis
 504 set xrange [xmin: xmax]
 505 set yrange [0: ymax]
 506 unset key
 507 set zeroaxis
 508 set xlabel "x ->"
 509 set ylabel "probability density ->"
 510 set xtics autofreq
 511 set ytics autofreq
 512 set format x "%.1f"
 513 set format y "%.1f"
 514 set sample 100
 515 set title "logistic PDF with a = 0, lambda = 2"
 516 plot logistic(x, a, lambda)
 517 pause -1 "Hit return to continue"
 518 set title "logistic CDF with a = 0, lambda = 2"
 519 set yrange [0: 1.1]
 520 plot clogistic(x, a, lambda)
 521 pause -1 "Hit return to continue"
 522
 523 # Lognormal PDF and CDF
 524 mu = 1.0; sigma = 0.5
 525 m = exp(mu + 0.5 * sigma**2)
 526 s = sqrt(exp(2.0 * mu + sigma**2) * (2.0 * exp(sigma) - 1.0))
 527 xmin = m - 4.0 * s
 528 xmin = xmin < 0 ? 0 : xmin
 529 xmax = m + 4.0 * s
 530 ymax = 1.1 * lognormal(exp(mu - sigma**2), mu, sigma) #Mode of lognormal PDF used
 531 unset key
 532 set zeroaxis
 533 set xrange [xmin: xmax]
 534 set yrange [0: ymax]
 535 set xlabel "x ->"
 536 set ylabel "probability density ->"
 537 set xtics autofreq
 538 set ytics autofreq
 539 set format x "%.2f"
 540 set format y "%.2f"
 541 set sample 100
 542 set title "lognormal PDF with mu = 1.0, sigma = 0.5"
 543 plot lognormal(x, mu, sigma)
 544 pause -1 "Hit return to continue"
 545 set title "lognormal CDF with mu = 1.0, sigma = 0.5"
 546 set yrange [0: 1.1]
 547 plot clognormal(x, mu, sigma)
 548 pause -1 "Hit return to continue"
 549
 550 # Maxwell PDF
 551 a = 0.5
 552 mu = 2.0 / sqrt(pi) / a
 553 sigma = sqrt(3.0 - 8.0/pi) / a
 554 xmin = int(mu - 3.0 * sigma)
 555 xmin = xmin < 0 ? 0 : xmin
 556 xmax = int(mu + 3.0 * sigma)
 557 a = 1.5
 558 ymax = 1.1 * maxwell(1.0 / a, a) + 0.5 #Mode of maxwell PDF used
 559 ymax = int(ymax + 0.5)
 560 set zeroaxis
 561 set xlabel "x ->"
 562 set ylabel "probability density ->"
 563 set xtics autofreq
 564 set ytics autofreq
 565 set format x "%.1f"
 566 set format y "%.1f"
 567 set sample 100
 568 set title "maxwell PDF"
 569 set key right top box
 570 plot [xmin:xmax] [0:ymax] maxwell(x, 1.5) title "a = 1.5", \
 571                    maxwell(x, 1.0) title "a = 1.0", \
 572                    maxwell(x, 0.5) title "a = 0.5"
 573 pause -1 "Hit return to continue"
 574 set title "maxwell CDF"
 575 set key right bottom box
 576 plot [xmin:xmax] [0:1.1] cmaxwell(x, 1.5) title "a = 1.5", \
 577                          cmaxwell(x, 1.0) title "a = 1.0", \
 578                          cmaxwell(x, 0.5) title "a = 0.5"
 579 pause -1 "Hit return to continue"
 580
 581 # Negative binomial PDF and CDF
 582 r = 8; p = 0.4
 583 mu = r * (1.0 - p) / p
 584 sigma = sqrt(mu / p)
 585 xmin = int(mu - 4.0 * sigma)
 586 xmin = xmin < 0 ? 0 : xmin
 587 xmax = int(mu + 4.0 * sigma)
 588 ymax = 1.1 * negbin(int(mu - (1.0-p)/p), r, p) #mode of gamma PDF used
 589 unset key
 590 unset zeroaxis
 591 set xrange [xmin-1 : xmax]
 592 set yrange [0 : ymax]
 593 set xlabel "k ->"
 594 set ylabel "probability density ->"
 595 set ytics 0, ymax / 10, ymax
 596 set format x "%2.0f"
 597 set format y "%3.2f"
 598 set sample (xmax - xmin+1) + 1
 599 set title "negative binomial (or pascal or polya) PDF with r = 8, p = 0.4"
 600 plot negbin(x, r, p) with impulses
 601 pause -1 "Hit return to continue"
 602 set yrange [0 : 1.1]
 603 set ytics 0, 0.1, 1.0
 604 set title "negative binomial (or pascal or polya) CDF with r = 8, p = 0.4"
 605 plot cnegbin(x, r, p) with steps
 606 pause -1 "Hit return to continue"
 607
 608 # Negative exponential PDF and CDF
 609 lambda = 2.0
 610 mu = 1.0 / lambda
 611 sigma = 1.0 / lambda
 612 xmax =  mu + 4.0 * sigma
 613 ymax = lambda #No mode
 614 unset key
 615 set zeroaxis
 616 set xrange [0: xmax]
 617 set yrange [0: ymax]
 618 set xlabel "x ->"
 619 set ylabel "probability density ->"
 620 set xtics autofreq
 621 set ytics autofreq
 622 set format x "%.2f"
 623 set format y "%.1f"
 624 set sample 100
 625 set title "negative exponential (or exponential) PDF with lambda = 2.0"
 626 plot nexp(x, lambda)
 627 pause -1 "Hit return to continue"
 628 set title "negative exponential (or exponential) CDF with lambda = 2.0"
 629 set yrange [0: 1.1]
 630 plot cnexp(x, lambda)
 631 pause -1 "Hit return to continue"
 632
 633 # Normal PDF and CDF
 634 mu = 0.0; sigma = 1.0
 635 xmin = mu - 4.0 * sigma
 636 xmax = mu + 4.0 * sigma
 637 mu = 2.0; sigma = 0.5
 638 ymax = 1.1 * normal(mu, mu, sigma) #Mode of normal PDF used
 639 set zeroaxis
 640 set xlabel "x ->"
 641 set ylabel "probability density ->"
 642 set xtics autofreq
 643 set ytics autofreq
 644 set format x "%.1f"
 645 set format y "%.1f"
 646 set sample 100
 647 set title "normal (also called gauss or bell-curved) PDF"
 648 set key left top box
 649 plot [xmin:xmax] [0:ymax] normal(x, 0, 1.0) title "mu = 0, sigma = 1.0", \
 650                           normal(x, 2, 0.5) title "mu = 2, sigma = 0.5", \
 651                           normal(x, 1, 2.0) title "mu = 1, sigma = 2.0"
 652 pause -1 "Hit return to continue"
 653 set title "normal (also called gauss or bell-curved) CDF"
 654 set key left top box
 655 plot [xmin:xmax] [0:1.1] mu = 0, sigma = 1.0, cnormal(x, mu, sigma) title "mu = 0, sigma = 1.0", \
 656                          mu = 2, sigma = 0.5, cnormal(x, mu, sigma) title "mu = 2, sigma = 0.5", \
 657                          mu = 1, sigma = 2.0, cnormal(x, mu, sigma) title "mu = 1, sigma = 2.0"
 658 pause -1 "Hit return to continue"
 659
 660 # Pareto PDF and CDF
 661 a = 1.0; b = 3.0
 662 mu = a * b / (b - 1.0)
 663 sigma = a * sqrt(b) / (sqrt(b - 2.0) * (b - 1.0))
 664 xmin = mu - 4.0 * sigma
 665 xmin = xmin < 0 ? 0 : xmin
 666 xmax = int(mu + 4.0 * sigma)
 667 ymax = 1.1 * pareto(a, a, b) #mode of pareto PDF used
 668 ymin = -0.1 * pareto(a, a, b)
 669 unset key
 670 set zeroaxis
 671 set xrange [xmin: xmax]
 672 set yrange [ymin: ymax]
 673 set xlabel "x ->"
 674 set ylabel "probability density ->"
 675 set xtics autofreq
 676 set ytics autofreq
 677 set format x "%.1f"
 678 set format y "%.1f"
 679 set sample 200+1
 680 set title "pareto PDF with a = 1, b = 3"
 681 # Discontinuity at a
 682 set parametric
 683 set trange [0:1-eps]
 684 x1(t) = -1 + 2*t
 685 x2(t) =  1 + 3*t
 686 set arrow 1 from 1.75,0.8 to 1.0,0.8
 687 set arrow 2 from 1.0,0.0 to 1.0,3.0 nohead lt 0
 688 set label 1 "Discontinuity achieved by plotting twice\nwith affine mapped parametric ranges" at 1.8,0.8 left
 689 plot x1(t), pareto(x1(t), a, b) ls 1, x2(t), pareto(x2(t), a, b) ls 1
 690 pause -1 "Hit return to continue"
 691 unset arrow 2
 692 set title "pareto CDF with a = 1, b = 3"
 693 unset parametric
 694 set yrange [-0.1: 1.1]
 695 set arrow 1 from 1.45,0.1 to 1.05,0.01
 696 set label 1 "Cusp achieved by selecting point\nas part of function samples" at 1.5,0.1 left
 697 plot cpareto(x, a, b)
 698 pause -1 "Hit return to continue"
 699 unset label 1
 700 unset arrow 1
 701
 702 # Poisson PDF and CDF
 703 mu = 4.0
 704 sigma = sqrt(mu)
 705 xmin = int(mu - 4.0 * sigma)
 706 xmin = xmin < -1 ? -1 : xmin
 707 xmax = int(mu + 4.0 * sigma)
 708 ymax = 1.1 * poisson(mu, mu) #mode of poisson PDF used
 709 unset key
 710 set zeroaxis
 711 set xrange [xmin : xmax]
 712 set yrange [0 : ymax]
 713 set xlabel "k ->"
 714 set ylabel "probability density ->"
 715 set ytics 0, ymax / 10, ymax
 716 set format x "%2.0f"
 717 set format y "%3.2f"
 718 set sample (xmax - xmin) + 1
 719 set title "poisson PDF with mu = 4.0"
 720 plot poisson(x, mu) with impulses
 721 pause -1 "Hit return to continue"
 722 set yrange [-0.1 : 1.1]
 723 set ytics -0.1, 0.1, 1.1
 724 set title "poisson CDF with mu = 4.0"
 725 plot cpoisson(x, mu) with steps
 726 pause -1 "Hit return to continue"
 727
 728 # Rayleigh PDF and CDF
 729 lambda = 2.0
 730 mu = 0.5 * sqrt(pi / lambda)
 731 sigma = sqrt((1.0 - pi / 4.0) / lambda)
 732 xmax = mu + 4.0 * sigma
 733 ymax = 1.1 * rayleigh(1.0 / sqrt(2.0 * lambda), lambda) #Mode of rayleigh PDF used
 734 unset key
 735 set zeroaxis
 736 set xrange [0: xmax]
 737 set yrange [0: ymax]
 738 set xlabel "x ->"
 739 set ylabel "probability density ->"
 740 set xtics autofreq
 741 set ytics autofreq
 742 set format x "%.2f"
 743 set format y "%.1f"
 744 set sample 100
 745 set title "rayleigh PDF with lambda = 2.0"
 746 plot rayleigh(x, lambda)
 747 pause -1 "Hit return to continue"
 748 set title "rayleigh CDF with lambda = 2.0"
 749 set yrange [0: 1.1]
 750 plot crayleigh(x, lambda)
 751 pause -1 "Hit return to continue"
 752
 753 # Sine PDF and CDF
 754 a = 3.2; f = 2.6
 755 mu = a / 2.0
 756 sigma = sqrt(a * a / 3.0 * (1.0 - 3.0 / (2.0 * n * n * pi * pi)) - mu * mu)
 757 xmin = 0.0
 758 xmax = a - eps
 759 a = 2; f = 1.0
 760 ymax = 1.1 * 2.0 / a #Mode of sine PDF used
 761 set zeroaxis
 762 set xlabel "x ->"
 763 set ylabel "probability density ->"
 764 set xtics autofreq
 765 set ytics autofreq
 766 set format x "%.2f"
 767 set format y "%.1f"
 768 set sample 250
 769 set title "sine PDF"
 770 set key bottom outside
 771 keystr(a, f) = sprintf("a = %0.1f, f = %0.1f", a, f)
 772 a1 = 2.0; a2 = 3.25; a3 = 2.75
 773 f1 = 1.0; f2 = 3.0; f3 = 2.6; f4 = 0.0
 774 plot [xmin:xmax] [0:ymax] a = a1, f = f1, sine(x, f, a) title keystr(a, f), \
 775                           a = a1, f = f2, sine(x, f, a) title keystr(a, f), \
 776                           a = a2, f = f3, sine(x, f, a) title keystr(a, f), \
 777                           a = a3, f = f4, sine(x, f, a) title keystr(a, f) with steps
 778 pause -1 "Hit return to continue"
 779 set title "sine CDF"
 780 set key top left
 781 plot [xmin:xmax] [0:1.1] a = a1, f = f1, csine(x, f, a) title keystr(a, f), \
 782                          a = a1, f = f2, csine(x, f, a) title keystr(a, f), \
 783                          a = a2, f = f3, csine(x, f, a) title keystr(a, f), \
 784                          a = a3, f = f4, csine(x, f, a) title keystr(a, f) with steps
 785 pause -1 "Hit return to continue"
 786
 787 # t PDF and CDF
 788 nu = 20
 789 mu = 0
 790 sigma = nu > 2 ? sqrt(nu / (nu - 2.0)) : 1.0
 791 xmin = mu - 4.0 * sigma
 792 xmax = mu + 4.0 * sigma
 793 ymax = 1.1 * t(mu, nu) #Mode of t PDF used
 794 set key inside center left title "degrees of freedom"
 795 set zeroaxis
 796 set xrange [xmin: xmax]
 797 set yrange [0: ymax]
 798 set xlabel "x ->"
 799 set ylabel "probability density ->"
 800 set xtics autofreq
 801 set ytics autofreq
 802 set format x "%.1f"
 803 set format y "%.2f"
 804 set sample 100
 805 set title "t PDF (and Gaussian limit)"
 806 ks(nu) = sprintf("nu = %d", nu)
 807 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), \
 808      t(x, 20) ti ks(20), normal(x, 0, 1) ti "normal"
 809 pause -1 "Hit return to continue"
 810 set title "t CDF (and Gaussian limit)"
 811 set yrange [0: 1.1]
 812 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), \
 813      ct(x, 20) ti ks(20), cnormal(x, 0, 1) ti "normal"
 814 pause -1 "Hit return to continue"
 815
 816 # Thanks to efrank@upenn5.hep.upenn.edu for telling us about this
 817 # triangular PDF and CDF
 818 m = 3.0
 819 g = 2.0
 820 mu = m
 821 sigma = g/sqrt(6.0)
 822 xmin = m - 1.1*g
 823 xmax = m + 1.1*g
 824 ymax = 1.1 * triangular(m, m, g) #Mode of triangular PDF used
 825 ymin = -ymax/11.0;
 826 unset key
 827 set zeroaxis
 828 set xrange [xmin: xmax]
 829 set yrange [ymin: ymax]
 830 set xlabel "x ->"
 831 set ylabel "probability density ->"
 832 set xtics autofreq
 833 set ytics autofreq
 834 set format x "%.1f"
 835 set format y "%.2f"
 836 set sample 50*1.1*g+1
 837 set title "triangular PDF with m = 3.0, g = 2.0"
 838 plot triangular(x, m, g)
 839 pause -1 "Hit return to continue"
 840 set title "triangular CDF with m = 3.0, g = 2.0"
 841 set yrange [-0.1: 1.1]
 842 plot ctriangular(x, m, g)
 843 pause -1 "Hit return to continue"
 844
 845 # Uniform PDF and CDF
 846 a = -2.0; b= 2.0
 847 mu = (a + b) / 2.0
 848 sigma = (b - a) / sqrt(12.0)
 849 xmin = a - 0.1*(b - a)
 850 xmax = b + 0.1*(b - a)
 851 ymax = 1.1 * uniform(mu, a, b) #No mode
 852 ymin = -0.1 * uniform(mu, a, b)
 853 unset key
 854 set zeroaxis
 855 set xrange [xmin: xmax]
 856 set yrange [ymin: ymax]
 857 set xlabel "x ->"
 858 set ylabel "probability density ->"
 859 set xtics autofreq
 860 set ytics autofreq
 861 set format x "%.2f"
 862 set format y "%.2f"
 863 set sample 120+1
 864 set title "uniform PDF with a = -2.0, b = 2.0"
 865 plot uniform(x, a, b) with steps
 866 pause -1 "Hit return to continue"
 867 set title "uniform CDF with a = -2.0, b = 2.0"
 868 set yrange [-0.1 : 1.1]
 869 plot cuniform(x, a, b)
 870 pause -1 "Hit return to continue"
 871
 872 # Weibull PDF and CDF
 873 lambda = 1.0/5; a = 1.0
 874 mu = 1.0 / lambda * gamma(1.0 / a) / a
 875 sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2))
 876 xmin = mu - 4.0 * sigma
 877 xmin = xmin < 0 ? 0 : xmin
 878 #Mode of weibull PDF used
 879 ymax = 1.8 * (a >= 1.0 ? weibull(((a - 1.0) / a)**(1.0 / a) / lambda, a, lambda) : 2.0)
 880 lambda = 1.0/15; a = 10.0
 881 mu = 1.0 / lambda * gamma(1.0 / a) / a
 882 sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2))
 883 xmax = int(mu + 4.0 * sigma)
 884 set key on title "" inside top right
 885 set zeroaxis
 886 set grid
 887 set xrange [xmin : xmax]
 888 set xlabel "x ->"
 889 set ylabel "probability density ->"
 890 set xtics autofreq
 891 set ytics autofreq
 892 set format x "%g"
 893 set format y "%g"
 894 set sample 100
 895 set title "weibull PDF"
 896 ks(a,lambda) = sprintf("lambda = 1/%g, a = %0.1f", 1.0/lambda, a)
 897 a1 = 0.5; a2 = 1.0; a3 = 2.0; a4 = 10.0
 898 lambda1 = 1.0/5; lambda2 = 1.0/15
 899 set arrow 1 from 3.8,0.27 to 0.5,weibull(0.5,a1,lambda1)
 900 set label 1 "a < 1, rate descreasing over time" at 4,0.27 left
 901 set arrow 2 from 8,0.19 to 6.4,weibull(6.4,a3,lambda1)
 902 set arrow 3 from 10.5,0.19 to 13,weibull(13,a4,lambda2)
 903 set label 2 "a > 1, rate increasing over time" at 9,0.2 center
 904 plot [] [0:ymax] lambda = lambda1, a = a1, weibull(x, a, lambda) ti ks(a, lambda), \
 905                  lambda = lambda1, a = a2, weibull(x, a, lambda) ti ks(a, lambda), \
 906                  lambda = lambda1, a = a3, weibull(x, a, lambda) ti ks(a, lambda), \
 907                  lambda = lambda2, a = a4, weibull(x, a, lambda) ti ks(a, lambda)
 908 pause -1 "Hit return to continue"
 909 unset label 1; unset label 2
 910 unset arrow 1; unset arrow 2; unset arrow 3
 911 set key at 9,0.4 center
 912 set title "weibull CDF"
 913 plot [] [0:1.1] lambda = lambda1, a = a1, cweibull(x, a, lambda) ti ks(a, lambda), \
 914                 lambda = lambda1, a = a2, cweibull(x, a, lambda) ti ks(a, lambda), \
 915                 lambda = lambda1, a = a3, cweibull(x, a, lambda) ti ks(a, lambda), \
 916                 lambda = lambda2, a = a4, cweibull(x, a, lambda) ti ks(a, lambda)
 917 pause -1 "Hit return to continue"
 918 reset