# 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 " Statistical Approximations, version 1.1" print "" print " Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl" print "" save_encoding = GPVAL_ENCODING set encoding utf8 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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot binom(rnd(x), n, p) with histeps, normal(x, mu, sigma) t "normal(x, μ, σ)" |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot binom(x, n, p) with histeps, poisson(x, mu) with histeps t "poisson(x, μ)" |
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 dt 3 set arrow from mu, gmm(mu + sigma, rho, lambda) \ to mu + sigma, gmm(mu + sigma, rho, lambda) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, gmm(mu + sigma, rho, lambda) offset 1, 0 plot geometric(rnd(x),p) with histeps, gmm(x, rho, lambda) t "gmm(x, ρ, λ)" |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot geometric(rnd(x), p) with histeps, normal(x, mu, sigma) t "normal(x, μ, σ)" |
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 dt 3 set arrow from mu, binom(floor(mu + sigma), n, p) \ to mu + sigma, binom(floor(mu + sigma), n, p) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, binom(floor(mu + sigma), n, p) offset 1, 0 plot hypgeo(x, nn, mm, n) with histeps, binom(x, n, p) with histeps |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot hypgeo(rnd(x), nn, mm, n) with histeps, normal(x, mu, sigma) t "normal(x, μ, σ)" |
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 dt 3 set arrow from mu, gmm(mu + sigma, rho, lambda) \ to mu + sigma, gmm(mu + sigma, rho, lambda) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, gmm(mu + sigma, rho, lambda) offset 1, 0 plot negbin(rnd(x), r, p) with histeps, gmm(x, rho, lambda) t "gmm(x, ρ, λ)" |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot negbin(rnd(x), r, p) with histeps, normal(x, mu, sigma) t "normal(x, μ, σ)" |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot logistic(x, a, lambda) t "logistic(x, a, λ)", normal(x, mu, sigma) t "normal(x, μ, σ)" |
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 dt 3 set arrow from mu, normal(mu + sigma, mu, sigma) \ to mu + sigma, normal(mu + sigma, mu, sigma) nohead dt 3 set label "μ" at mu, ymax / 10 offset 1, 0 set label "σ" at mu + sigma, normal(mu + sigma, mu, sigma) offset 1, 0 plot \ poisson(rnd(x), mu) with histeps t "poisson(rnd(x), μ)", \ normal(x, mu, sigma) t "normal(x, μ, σ)" |
reset set encoding save_encoding |