2 # $Id: prob2.dem,v 1.9 2006/06/14 03:24:09 sfeam Exp $
4 # Demo Statistical Approximations version 1.1
6 # Copyright (c) 1991, Jos van der Woude, jvdwoude@hut.nl
9 # -- --- 1991 Jos van der Woude: 1st version
10 # 06 Jun 2006 Dan Sebald: Added plot methods for better visual effect.
18 print " Statistical Approximations, version 1.1"
20 print " Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl"
32 print " NOTE: contains 10 plots and consequently takes some time to run"
33 print " Press Ctrl-C to exit right now"
35 pause -1 " Press Return to start demo ..."
42 # Binomial PDF using normal approximation
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
52 set xrange [xmin - 1 : xmax + 1]
55 set ylabel "probability density ->"
56 set ytics 0, ymax / 10.0, ymax
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"
71 # Binomial PDF using poisson approximation
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
81 set xrange [xmin - 1 : xmax + 1]
84 set ylabel "probability density ->"
85 set ytics 0, ymax / 10.0, ymax
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"
100 # Geometric PDF using gamma approximation
106 xmin = floor(mu - r_sigma * sigma)
107 xmin = xmin < r_xmin ? r_xmin : xmin
108 xmax = ceil(mu + r_sigma * sigma)
112 set xrange [xmin - 1 : xmax + 1]
113 set yrange [0 : ymax]
115 set ylabel "probability density ->"
116 set ytics 0, ymax / 10.0, ymax
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"
131 # Geometric PDF using normal approximation
135 xmin = floor(mu - r_sigma * sigma)
136 xmin = xmin < r_xmin ? r_xmin : xmin
137 xmax = ceil(mu + r_sigma * sigma)
141 set xrange [xmin - 1 : xmax + 1]
142 set yrange [0 : ymax]
144 set ylabel "probability density ->"
145 set ytics 0, ymax / 10.0, ymax
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"
160 # Hypergeometric PDF using binomial approximation
161 nn = 75; mm = 25; n = 10
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
171 set xrange [xmin - 1 : xmax + 1]
172 set yrange [0 : ymax]
174 set ylabel "probability density ->"
175 set ytics 0, ymax / 10.0, ymax
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"
190 # Hypergeometric PDF using normal approximation
191 nn = 75; mm = 25; n = 10
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
201 set xrange [xmin - 1 : xmax + 1]
202 set yrange [0 : ymax]
204 set ylabel "probability density ->"
205 set ytics 0, ymax / 10.0, ymax
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"
220 # Negative binomial PDF using gamma approximation
222 mu = r * (1.0 - p) / 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
232 set xrange [xmin - 1 : xmax + 1]
233 set yrange [0 : ymax]
235 set ylabel "probability density ->"
236 set ytics 0, ymax / 10.0, ymax
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"
251 # Negative binomial PDF using normal approximation
253 mu = r * (1.0 - p) / 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
261 set xrange [xmin - 1 : xmax + 1]
262 set yrange [0 : ymax]
264 set ylabel "probability density ->"
265 set ytics 0, ymax / 10.0, ymax
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"
280 # Normal PDF using logistic approximation
281 mu = 1.0; sigma = 1.5
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
289 set xrange [xmin: xmax]
290 set yrange [0 : ymax]
292 set ylabel "probability density ->"
293 set ytics 0, ymax / 10.0, ymax
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"
308 # Poisson PDF using normal approximation
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
317 set xrange [xmin - 1 : xmax + 1]
318 set yrange [0 : ymax]
320 set ylabel "probability density ->"
321 set ytics 0, ymax / 10.0, ymax
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"