334 lines
10 KiB
Gnuplot
334 lines
10 KiB
Gnuplot
#
|
|
# $Id: prob2.dem,v 1.9 2006/06/14 03:24:09 sfeam Exp $
|
|
#
|
|
# Demo Statistical Approximations version 1.1
|
|
#
|
|
# Copyright (c) 1991, Jos van der Woude, jvdwoude@hut.nl
|
|
|
|
# History:
|
|
# -- --- 1991 Jos van der Woude: 1st version
|
|
# 06 Jun 2006 Dan Sebald: Added plot methods for better visual effect.
|
|
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print " Statistical Approximations, version 1.1"
|
|
print ""
|
|
print " Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl"
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print ""
|
|
print " NOTE: contains 10 plots and consequently takes some time to run"
|
|
print " Press Ctrl-C to exit right now"
|
|
print ""
|
|
pause -1 " Press Return to start demo ..."
|
|
|
|
load "stat.inc"
|
|
rnd(x) = floor(x+0.5)
|
|
r_xmin = -1
|
|
r_sigma = 4.0
|
|
|
|
# Binomial PDF using normal approximation
|
|
n = 25; p = 0.15
|
|
mu = n * p
|
|
sigma = sqrt(n * p * (1.0 - p))
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * binom(floor((n+1)*p), n, p) #mode of binomial PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "binomial PDF using normal approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot binom(rnd(x), n, p) with histeps, normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Binomial PDF using poisson approximation
|
|
n = 50; p = 0.1
|
|
mu = n * p
|
|
sigma = sqrt(mu)
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * binom(floor((n+1)*p), n, p) #mode of binomial PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample (xmax - xmin + 3)
|
|
set title "binomial PDF using poisson approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot binom(x, n, p) with histeps, poisson(x, mu) with histeps
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Geometric PDF using gamma approximation
|
|
p = 0.3
|
|
mu = (1.0 - p) / p
|
|
sigma = sqrt(mu / p)
|
|
lambda = p
|
|
rho = 1.0 - p
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * p
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "geometric PDF using gamma approximation"
|
|
set arrow from mu, 0 to mu, gmm(mu, rho, lambda) nohead
|
|
set arrow from mu, gmm(mu + sigma, rho, lambda) \
|
|
to mu + sigma, gmm(mu + sigma, rho, lambda) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, gmm(mu + sigma, rho, lambda)
|
|
plot geometric(rnd(x),p) with histeps, gmm(x, rho, lambda)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Geometric PDF using normal approximation
|
|
p = 0.3
|
|
mu = (1.0 - p) / p
|
|
sigma = sqrt(mu / p)
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * p
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "geometric PDF using normal approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot geometric(rnd(x),p) with histeps, normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Hypergeometric PDF using binomial approximation
|
|
nn = 75; mm = 25; n = 10
|
|
p = real(mm) / nn
|
|
mu = n * p
|
|
sigma = sqrt(real(nn - n) / (nn - 1.0) * n * p * (1.0 - p))
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * hypgeo(floor(mu), nn, mm, n) #mode of binom PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample (xmax - xmin + 3)
|
|
set title "hypergeometric PDF using binomial approximation"
|
|
set arrow from mu, 0 to mu, binom(floor(mu), n, p) nohead
|
|
set arrow from mu, binom(floor(mu + sigma), n, p) \
|
|
to mu + sigma, binom(floor(mu + sigma), n, p) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, binom(floor(mu + sigma), n, p)
|
|
plot hypgeo(x, nn, mm, n) with histeps, binom(x, n, p) with histeps
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Hypergeometric PDF using normal approximation
|
|
nn = 75; mm = 25; n = 10
|
|
p = real(mm) / nn
|
|
mu = n * p
|
|
sigma = sqrt(real(nn - n) / (nn - 1.0) * n * p * (1.0 - p))
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * hypgeo(floor(mu), nn, mm, n) #mode of binom PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "hypergeometric PDF using normal approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot hypgeo(rnd(x), nn, mm, n) with histeps, normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Negative binomial PDF using gamma approximation
|
|
r = 8; p = 0.6
|
|
mu = r * (1.0 - p) / p
|
|
sigma = sqrt(mu / p)
|
|
lambda = p
|
|
rho = r * (1.0 - p)
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * gmm((rho - 1) / lambda, rho, lambda) #mode of gamma PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "negative binomial PDF using gamma approximation"
|
|
set arrow from mu, 0 to mu, gmm(mu, rho, lambda) nohead
|
|
set arrow from mu, gmm(mu + sigma, rho, lambda) \
|
|
to mu + sigma, gmm(mu + sigma, rho, lambda) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, gmm(mu + sigma, rho, lambda)
|
|
plot negbin(rnd(x), r, p) with histeps, gmm(x, rho, lambda)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Negative binomial PDF using normal approximation
|
|
r = 8; p = 0.4
|
|
mu = r * (1.0 - p) / p
|
|
sigma = sqrt(mu / p)
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * negbin(floor((r-1)*(1-p)/p), r, p) #mode of gamma PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "negative binomial PDF using normal approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot negbin(rnd(x), r, p) with histeps, normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Normal PDF using logistic approximation
|
|
mu = 1.0; sigma = 1.5
|
|
a = mu
|
|
lambda = pi / (sqrt(3.0) * sigma)
|
|
xmin = mu - r_sigma * sigma
|
|
xmax = mu + r_sigma * sigma
|
|
ymax = 1.1 * logistic(mu, a, lambda) #mode of logistic PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin: xmax]
|
|
set yrange [0 : ymax]
|
|
set xlabel "x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%.1f"
|
|
set format y "%.2f"
|
|
set sample 200
|
|
set title "normal PDF using logistic approximation"
|
|
set arrow from mu,0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot logistic(x, a, lambda), normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
unset arrow
|
|
unset label
|
|
|
|
# Poisson PDF using normal approximation
|
|
mu = 5.0
|
|
sigma = sqrt(mu)
|
|
xmin = floor(mu - r_sigma * sigma)
|
|
xmin = xmin < r_xmin ? r_xmin : xmin
|
|
xmax = ceil(mu + r_sigma * sigma)
|
|
ymax = 1.1 * poisson(mu, mu) #mode of poisson PDF used
|
|
set key box
|
|
unset zeroaxis
|
|
set xrange [xmin - 1 : xmax + 1]
|
|
set yrange [0 : ymax]
|
|
set xlabel "k, x ->"
|
|
set ylabel "probability density ->"
|
|
set ytics 0, ymax / 10.0, ymax
|
|
set format x "%2.0f"
|
|
set format y "%3.2f"
|
|
set sample 200
|
|
set title "poisson PDF using normal approximation"
|
|
set arrow from mu, 0 to mu, normal(mu, mu, sigma) nohead
|
|
set arrow from mu, normal(mu + sigma, mu, sigma) \
|
|
to mu + sigma, normal(mu + sigma, mu, sigma) nohead
|
|
set label "mu" at mu + 0.5, ymax / 10
|
|
set label "sigma" at mu + 0.5 + sigma, normal(mu + sigma, mu, sigma)
|
|
plot poisson(rnd(x), mu) with histeps, normal(x, mu, sigma)
|
|
pause -1 "Hit return to continue"
|
|
reset
|