git.maemo.org Git - gnuplot/blob - demo/prob2.dem

   1 #
   2 # $Id: prob2.dem,v 1.9 2006/06/14 03:24:09 sfeam Exp $
   3 #
   4 # Demo Statistical Approximations version 1.1
   5 #
   6 # Copyright (c) 1991, Jos van der Woude, jvdwoude@hut.nl
   7
   8 # History:
   9 #    -- --- 1991 Jos van der Woude:  1st version
  10 #    06 Jun 2006 Dan Sebald:  Added plot methods for better visual effect.
  11
  12 print ""
  13 print ""
  14 print ""
  15 print ""
  16 print ""
  17 print ""
  18 print "                        Statistical Approximations, version 1.1"
  19 print ""
  20 print "        Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl"
  21 print ""
  22 print ""
  23 print ""
  24 print ""
  25 print ""
  26 print ""
  27 print ""
  28 print ""
  29 print ""
  30 print ""
  31 print ""
  32 print "     NOTE: contains 10 plots and consequently takes some time to run"
  33 print "                      Press Ctrl-C to exit right now"
  34 print ""
  35 pause -1 "                      Press Return to start demo ..."
  36
  37 load "stat.inc"
  38 rnd(x) = floor(x+0.5)
  39 r_xmin = -1
  40 r_sigma = 4.0
  41
  42 # Binomial PDF using normal approximation
  43 n = 25; p = 0.15
  44 mu = n * p
  45 sigma = sqrt(n * p * (1.0 - p))
  46 xmin = floor(mu - r_sigma * sigma)
  47 xmin = xmin < r_xmin ? r_xmin : xmin
  48 xmax = ceil(mu + r_sigma * sigma)
  49 ymax = 1.1 * binom(floor((n+1)*p), n, p) #mode of binomial PDF used
  50 set key box
  51 unset zeroaxis
  52 set xrange [xmin - 1 : xmax + 1]
  53 set yrange [0 : ymax]
  54 set xlabel "k, x ->"
  55 set ylabel "probability density ->"
  56 set ytics 0, ymax / 10.0, ymax
  57 set format x "%2.0f"
  58 set format y "%3.2f"
  59 set sample 200
  60 set title "binomial PDF using normal approximation"
  61 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
  62 set arrow from mu, normal(mu + sigma, mu, sigma) \
  63           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
  64 set label "mu" at mu + 0.5, ymax / 10
  65 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
  66 plot binom(rnd(x), n, p) with histeps, normal(x, mu, sigma)
  67 pause -1 "Hit return to continue"
  68 unset arrow
  69 unset label
  70
  71 # Binomial PDF using poisson approximation
  72 n = 50; p = 0.1
  73 mu = n * p
  74 sigma = sqrt(mu)
  75 xmin = floor(mu - r_sigma * sigma)
  76 xmin = xmin < r_xmin ? r_xmin : xmin
  77 xmax = ceil(mu + r_sigma * sigma)
  78 ymax = 1.1 * binom(floor((n+1)*p), n, p) #mode of binomial PDF used
  79 set key box
  80 unset zeroaxis
  81 set xrange [xmin - 1 : xmax + 1]
  82 set yrange [0 : ymax]
  83 set xlabel "k ->"
  84 set ylabel "probability density ->"
  85 set ytics 0, ymax / 10.0, ymax
  86 set format x "%2.0f"
  87 set format y "%3.2f"
  88 set sample (xmax - xmin + 3)
  89 set title "binomial PDF using poisson approximation"
  90 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
  91 set arrow from mu, normal(mu + sigma, mu, sigma) \
  92           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
  93 set label "mu" at mu + 0.5, ymax / 10
  94 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
  95 plot binom(x, n, p) with histeps, poisson(x, mu) with histeps
  96 pause -1 "Hit return to continue"
  97 unset arrow
  98 unset label
  99
 100 # Geometric PDF using gamma approximation
 101 p = 0.3
 102 mu = (1.0 - p) / p
 103 sigma = sqrt(mu / p)
 104 lambda = p
 105 rho = 1.0 - p
 106 xmin = floor(mu - r_sigma * sigma)
 107 xmin = xmin < r_xmin ? r_xmin : xmin
 108 xmax = ceil(mu + r_sigma * sigma)
 109 ymax = 1.1 * p
 110 set key box
 111 unset zeroaxis
 112 set xrange [xmin - 1 : xmax + 1]
 113 set yrange [0 : ymax]
 114 set xlabel "k, x ->"
 115 set ylabel "probability density ->"
 116 set ytics 0, ymax / 10.0, ymax
 117 set format x "%2.0f"
 118 set format y "%3.2f"
 119 set sample 200
 120 set title "geometric PDF using gamma approximation"
 121 set arrow from mu, 0 to mu, gmm(mu, rho, lambda) nohead
 122 set arrow from mu, gmm(mu + sigma, rho, lambda) \
 123           to mu + sigma, gmm(mu + sigma, rho, lambda) nohead
 124 set label "mu" at mu + 0.5, ymax / 10
 125 set label "sigma" at mu + 0.5 + sigma, gmm(mu + sigma, rho, lambda)
 126 plot geometric(rnd(x),p) with histeps, gmm(x, rho, lambda)
 127 pause -1 "Hit return to continue"
 128 unset arrow
 129 unset label
 130
 131 # Geometric PDF using normal approximation
 132 p = 0.3
 133 mu = (1.0 - p) / p
 134 sigma = sqrt(mu / p)
 135 xmin = floor(mu - r_sigma * sigma)
 136 xmin = xmin < r_xmin ? r_xmin : xmin
 137 xmax = ceil(mu + r_sigma * sigma)
 138 ymax = 1.1 * p
 139 set key box
 140 unset zeroaxis
 141 set xrange [xmin - 1 : xmax + 1]
 142 set yrange [0 : ymax]
 143 set xlabel "k, x ->"
 144 set ylabel "probability density ->"
 145 set ytics 0, ymax / 10.0, ymax
 146 set format x "%2.0f"
 147 set format y "%3.2f"
 148 set sample 200
 149 set title "geometric PDF using normal approximation"
 150 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
 151 set arrow from mu, normal(mu + sigma, mu, sigma) \
 152           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
 153 set label "mu" at mu + 0.5, ymax / 10
 154 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
 155 plot geometric(rnd(x),p) with histeps, normal(x, mu, sigma)
 156 pause -1 "Hit return to continue"
 157 unset arrow
 158 unset label
 159
 160 # Hypergeometric PDF using binomial approximation
 161 nn = 75; mm = 25; n = 10
 162 p = real(mm) / nn
 163 mu = n * p
 164 sigma = sqrt(real(nn - n) / (nn - 1.0) * n * p * (1.0 - p))
 165 xmin = floor(mu - r_sigma * sigma)
 166 xmin = xmin < r_xmin ? r_xmin : xmin
 167 xmax = ceil(mu + r_sigma * sigma)
 168 ymax = 1.1 * hypgeo(floor(mu), nn, mm, n) #mode of binom PDF used
 169 set key box
 170 unset zeroaxis
 171 set xrange [xmin - 1 : xmax + 1]
 172 set yrange [0 : ymax]
 173 set xlabel "k ->"
 174 set ylabel "probability density ->"
 175 set ytics 0, ymax / 10.0, ymax
 176 set format x "%2.0f"
 177 set format y "%3.2f"
 178 set sample (xmax - xmin + 3)
 179 set title "hypergeometric PDF using binomial approximation"
 180 set arrow from mu, 0 to mu, binom(floor(mu), n, p) nohead
 181 set arrow from mu, binom(floor(mu + sigma), n, p) \
 182           to mu + sigma, binom(floor(mu + sigma), n, p) nohead
 183 set label "mu" at mu + 0.5, ymax / 10
 184 set label "sigma" at mu + 0.5 + sigma, binom(floor(mu + sigma), n, p)
 185 plot hypgeo(x, nn, mm, n) with histeps, binom(x, n, p) with histeps
 186 pause -1 "Hit return to continue"
 187 unset arrow
 188 unset label
 189
 190 # Hypergeometric PDF using normal approximation
 191 nn = 75; mm = 25; n = 10
 192 p = real(mm) / nn
 193 mu = n * p
 194 sigma = sqrt(real(nn - n) / (nn - 1.0) * n * p * (1.0 - p))
 195 xmin = floor(mu - r_sigma * sigma)
 196 xmin = xmin < r_xmin ? r_xmin : xmin
 197 xmax = ceil(mu + r_sigma * sigma)
 198 ymax = 1.1 * hypgeo(floor(mu), nn, mm, n) #mode of binom PDF used
 199 set key box
 200 unset zeroaxis
 201 set xrange [xmin - 1 : xmax + 1]
 202 set yrange [0 : ymax]
 203 set xlabel "k, x ->"
 204 set ylabel "probability density ->"
 205 set ytics 0, ymax / 10.0, ymax
 206 set format x "%2.0f"
 207 set format y "%3.2f"
 208 set sample 200
 209 set title "hypergeometric PDF using normal approximation"
 210 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
 211 set arrow from mu, normal(mu + sigma, mu, sigma) \
 212           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
 213 set label "mu" at mu + 0.5, ymax / 10
 214 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
 215 plot hypgeo(rnd(x), nn, mm, n) with histeps, normal(x, mu, sigma)
 216 pause -1 "Hit return to continue"
 217 unset arrow
 218 unset label
 219
 220 # Negative binomial PDF using gamma approximation
 221 r = 8; p = 0.6
 222 mu = r * (1.0 - p) / p
 223 sigma = sqrt(mu / p)
 224 lambda = p
 225 rho = r * (1.0 - p)
 226 xmin = floor(mu - r_sigma * sigma)
 227 xmin = xmin < r_xmin ? r_xmin : xmin
 228 xmax = ceil(mu + r_sigma * sigma)
 229 ymax = 1.1 * gmm((rho - 1) / lambda, rho, lambda) #mode of gamma PDF used
 230 set key box
 231 unset zeroaxis
 232 set xrange [xmin - 1 : xmax + 1]
 233 set yrange [0 : ymax]
 234 set xlabel "k, x ->"
 235 set ylabel "probability density ->"
 236 set ytics 0, ymax / 10.0, ymax
 237 set format x "%2.0f"
 238 set format y "%3.2f"
 239 set sample 200
 240 set title "negative binomial PDF using gamma approximation"
 241 set arrow from mu, 0 to mu, gmm(mu, rho, lambda) nohead
 242 set arrow from mu, gmm(mu + sigma, rho, lambda) \
 243           to mu + sigma, gmm(mu + sigma, rho, lambda) nohead
 244 set label "mu" at mu + 0.5, ymax / 10
 245 set label "sigma" at mu + 0.5 + sigma, gmm(mu + sigma, rho, lambda)
 246 plot negbin(rnd(x), r, p) with histeps, gmm(x, rho, lambda)
 247 pause -1 "Hit return to continue"
 248 unset arrow
 249 unset label
 250
 251 # Negative binomial PDF using normal approximation
 252 r = 8; p = 0.4
 253 mu = r * (1.0 - p) / p
 254 sigma = sqrt(mu / p)
 255 xmin = floor(mu - r_sigma * sigma)
 256 xmin = xmin < r_xmin ? r_xmin : xmin
 257 xmax = ceil(mu + r_sigma * sigma)
 258 ymax = 1.1 * negbin(floor((r-1)*(1-p)/p), r, p) #mode of gamma PDF used
 259 set key box
 260 unset zeroaxis
 261 set xrange [xmin - 1 : xmax + 1]
 262 set yrange [0 : ymax]
 263 set xlabel "k, x ->"
 264 set ylabel "probability density ->"
 265 set ytics 0, ymax / 10.0, ymax
 266 set format x "%2.0f"
 267 set format y "%3.2f"
 268 set sample 200
 269 set title "negative binomial PDF using normal approximation"
 270 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
 271 set arrow from mu, normal(mu + sigma, mu, sigma) \
 272           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
 273 set label "mu" at mu + 0.5, ymax / 10
 274 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
 275 plot negbin(rnd(x), r, p) with histeps, normal(x, mu, sigma)
 276 pause -1 "Hit return to continue"
 277 unset arrow
 278 unset label
 279
 280 # Normal PDF using logistic approximation
 281 mu = 1.0; sigma = 1.5
 282 a = mu
 283 lambda = pi / (sqrt(3.0) * sigma)
 284 xmin = mu - r_sigma * sigma
 285 xmax = mu + r_sigma * sigma
 286 ymax = 1.1 * logistic(mu, a, lambda) #mode of logistic PDF used
 287 set key box
 288 unset zeroaxis
 289 set xrange [xmin: xmax]
 290 set yrange [0 : ymax]
 291 set xlabel "x ->"
 292 set ylabel "probability density ->"
 293 set ytics 0, ymax / 10.0, ymax
 294 set format x "%.1f"
 295 set format y "%.2f"
 296 set sample 200
 297 set title "normal PDF using logistic approximation"
 298 set arrow from mu,0 to mu, normal(mu, mu, sigma) nohead
 299 set arrow from mu, normal(mu + sigma, mu, sigma) \
 300           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
 301 set label "mu" at mu + 0.5, ymax / 10
 302 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
 303 plot logistic(x, a, lambda), normal(x, mu, sigma)
 304 pause -1 "Hit return to continue"
 305 unset arrow
 306 unset label
 307
 308 # Poisson PDF using normal approximation
 309 mu = 5.0
 310 sigma = sqrt(mu)
 311 xmin = floor(mu - r_sigma * sigma)
 312 xmin = xmin < r_xmin ? r_xmin : xmin
 313 xmax = ceil(mu + r_sigma * sigma)
 314 ymax = 1.1 * poisson(mu, mu) #mode of poisson PDF used
 315 set key box
 316 unset zeroaxis
 317 set xrange [xmin - 1 : xmax + 1]
 318 set yrange [0 : ymax]
 319 set xlabel "k, x ->"
 320 set ylabel "probability density ->"
 321 set ytics 0, ymax / 10.0, ymax
 322 set format x "%2.0f"
 323 set format y "%3.2f"
 324 set sample 200
 325 set title "poisson PDF using normal approximation"
 326 set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
 327 set arrow from mu, normal(mu + sigma, mu, sigma) \
 328           to mu + sigma, normal(mu + sigma, mu, sigma) nohead
 329 set label "mu" at mu + 0.5, ymax / 10
 330 set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
 331 plot poisson(rnd(x), mu) with histeps, normal(x, mu, sigma)
 332 pause -1 "Hit return to continue"
 333 reset