# Demo Statistical Functions version 2.3 # # Copyright (c) 1991, 1992 Jos van der Woude, jvdwoude@hut.nl # History: # -- --- 1992 Jos van der Woude: 1st version # 06 Jun 2006 Dan Sebald: Added some variety and plotting techniques for # better visual effect. More tutorial in nature. # 2019 BM: Use Unicode, loops and arrays print " Statistical Library Demo, version 2.3" print "\n Copyright (c) 1991, 1992, Jos van de Woude, jvdwoude@hut.nl" print "\n" print " Press Ctrl-C to exit right now" load "stat.inc" save_encoding = GPVAL_ENCODING set encoding utf8 eps = 1.0e-10 # Supposed to be float resolution (nice if were defined internally) ## Gamma function xmin = -5.5 xmax = 5 ymin = -10 ymax = 10 unset key set xzeroaxis gsampfunc(t,n) = t<0?0.5*1/(-t+1.0)**n:1.0-0.5*1/(t+1.0)**n set parametric set trange [-1:1] set sample 200 set xrange [xmin : xmax] set yrange [ymin : ymax] set xlabel "x" set ylabel "Γ(x)" set for [i=1:6] arrow i from 1-i,ymin to 1-i,ymax nohead lt 0 set title "Gamma function Γ, very useful function for probability" plot \ gsampfunc(5*t,5)-6, gamma(gsampfunc(5*t,5)-6) lt 1, \ for [i=5:1:-1] \ gsampfunc(5*t,i)-i, gamma(gsampfunc(5*t,i)-i) lt 1, \ 5*gsampfunc(5*t,2), gamma(5*gsampfunc(5*t,2)) lt 1 |
ymin = ymin/2 ymax = ymax/2 set yrange [ymin:ymax] set ylabel "lgamma(x)" set for [i=1:6] arrow i from 1-i,ymin to 1-i,ymax nohead lt 0 set title "log gamma function, similarly very useful function" plot gsampfunc(5*t,5)-6, lgamma(gsampfunc(5*t,5)-6) lt 1, \ gsampfunc(5*t,5)-5, lgamma(gsampfunc(5*t,5)-5) lt 1, \ gsampfunc(5*t,4)-4, lgamma(gsampfunc(5*t,4)-4) lt 1, \ gsampfunc(5*t,3)-3, lgamma(gsampfunc(5*t,3)-3) lt 1, \ gsampfunc(5*t,3)-2, lgamma(gsampfunc(5*t,3)-2) lt 1, \ gsampfunc(5*t,3)-1, lgamma(gsampfunc(5*t,3)-1) lt 1, \ 5*gsampfunc(5*t,3), lgamma(5*gsampfunc(5*t,3)) lt 1 |
unset for [i=1:6] arrow i unset sample unset parametric unset xzeroaxis set key # Arcsinus PDF and CDF r = 2.0 mu = 0.0 sigma = r / sqrt2 xmin = -(r+1) xmax = r+1 unset key set zeroaxis set xrange [xmin : xmax] set yrange [-0.2 : 1.2] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 50*xmax+1 set title "arcsin PDF with r = 2.0" plot arcsin(x, r) |
set title "arcsin CDF with r = 2.0" set yrange [-0.2 : 1.2] plot carcsin(x, r) |
# Beta PDF and CDF p = 5.0; q = 3.0 mu = p / (p + q) sigma = sqrt(p**q) / ((p + q ) * sqrt(p + q + 1.0)) xmin = 0.0 xmax = 1.0 #Mode of beta PDF used ymax = (p < 1.0 || q < 1.0) ? 2.0 : 1.4 * beta((p - 1.0)/(p + q - 2.0), p, q) set key right box set zeroaxis set xrange [xmin : xmax] set yrange [0 : ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 200 set title "beta PDF" plot beta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \ beta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \ beta(x, 0.5, 2.5) title "p = 0.5, q = 2.5" |
set yrange [0:1.1] set title "incomplete beta CDF" set key left box plot cbeta(x, 0.5, 0.7) title "p = 0.5, q = 0.7", \ cbeta(x, 5.0, 3.0) title "p = 5.0, q = 3.0", \ cbeta(x, 0.5, 2.5) title "p = 0.5, q = 2.5" |
# Binomial PDF and CDF n = 25; p = 0.15 mu = n * p sigma = sqrt(n * p * (1.0 - p)) xmin = int(mu - 4.0 * sigma) xmin = xmin < -2 ? -2 : xmin xmax = int(mu + 4.0 * sigma) xmax = xmax < n+2 ? n+2 : xmax ymax = 1.1 * binom(int((n+1)*p), n, p) #Mode of binomial PDF used unset key unset zeroaxis set xrange [xmin : xmax] set yrange [0 : ymax] set xlabel "k" set ylabel "probability density" set ytics 0, ymax / 10, ymax set format x "%2.0f" set format y "%3.2f" set sample (xmax - xmin) + 1 set title "binomial PDF with n = 25, p = 0.15" plot binom(x, n, p) with impulses |
set ytics autofreq set xzeroaxis set title "binomial CDF with n = 25, p = 0.15" set yrange [-0.1 : 1.1] set ytics 0, 0.1, 1.0 plot cbinom(x, n, p) with steps |
# Cauchy PDF and CDF a = 0.0; b = 2.0 #cauchy PDF has no moments xmin = a - 5.0 * b xmax = a + 5.0 * b ymax = 1.1 * cauchy(a, a, b) #Mode of cauchy PDF used set key left box set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100 set title "Cauchy PDF" a=0 b=2 plot [xmin:xmax] [0:ymax] cauchy(x, 0, 2) title "a = 0, b = 2", \ cauchy(x, 0, 4) title "a = 0, b = 4" |
set title "Cauchy CDF" plot [xmin:xmax] [0:1.0] ccauchy(x, 0, 2) title "a = 0, b = 2", \ ccauchy(x, 0, 4) title "a = 0, b = 4" |
# Chi-square PDF and CDF k = 4.0 mu = k sigma = sqrt(2.0 * k) xmin = mu - 4.0 * sigma xmin = xmin < 0 ? 0 : xmin xmax = int(mu + 4.0 * sigma) k = 2.0 ymax = (k > 2.0 ? 1.1*chisq(k - 2.0, k) : 0.5) #Mode of chi PDF used set key right box set zeroaxis set xrange [xmin+eps : xmax] #Discontinuity at zero for k < 2 set yrange [0:ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100 set title "Chi-square χ^2 PDF" set key right box set samples 15*20+1 keystr(k) = sprintf("k = %d", k) plot k = 1, x==0?1/0:chisq(x, k) title keystr(k), \ k = 2, x==0?1/0:chisq(x, k) title keystr(k), \ for [k = 3:8] chisq(x, k) title keystr(k) |
set yrange [0:1.1] set key bottom right box set title "Chi-square χ^2 CDF" plot for [k = 1:8] cchisq(x, k) title keystr(k) |
# Erlang PDF and CDF lambda = 1.0; n = 2.0 mu = n / lambda sigma = sqrt(n) / lambda xmax = int(mu + 5.0 * sigma) n = 1.0 ymax = n < 2.0 ? 1.0 : 1.1 * erlang((n - 1.0) / lambda, n, lambda) #Mode of erlang PDF used set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set title "Erlang PDF" set key top right box l1 = 1.0; l2 = 0.5 set arrow 1 from 2,0.8 to 0.33,erlang(0.33,1,l1) set arrow 2 from 2,0.8 to 0.33,erlang(0.33,1,l2) set label 1 "n = 1, exponential r.v." at 2.1,0.8 left keystr(n,lambda) = sprintf("λ = %0.1f, n = %d", lambda, n) plot [0:xmax] [0:ymax] n = 1, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \ n = 1, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda), \ n = 2, lambda = l1, erlang(x, n, lambda) title keystr(n,lambda), \ n = 2, lambda = l2, erlang(x, n, lambda) title keystr(n,lambda) |
unset label 1 unset arrow 1; unset arrow 2 set title "Erlang CDF" set key bottom right box plot [0:xmax] [0:1.1] n = 1, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \ n = 1, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda), \ n = 2, lambda = l1, cerlang(x, n, lambda) title keystr(n,lambda), \ n = 2, lambda = l2, cerlang(x, n, lambda) title keystr(n,lambda) |
# Thanks to mrb2j@kelvin.seas.Virginia.EDU for telling us about this. # Extreme (Gumbel extreme value) PDF and CDF alpha = 1.0; u = 0.0 mu = u + (0.577215665/alpha) # Euler's constant sigma = pi/(sqrt(6.0)*alpha) xmin = mu - 6.0 * sigma xmax = mu + 6.0 * sigma ymax = 1.1 * extreme(u, u, alpha) #Mode of extreme PDF used ymax = int(10*ymax)/10.0 set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100 set title "extreme PDF" set key top left box plot [xmin:xmax] [0:ymax] extreme(x, 1.0, 0.5) title "α = 0.5, u = 1.0", \ extreme(x, 0.0, 1.0) title "α = 1.0, u = 0.0" |
set title "extreme CDF" plot [xmin:xmax] [0:1.1] cextreme(x, 1.0, 0.5) title "α = 0.5, u = 1.0", \ cextreme(x, 0.0, 1.0) title "α = 1.0, u = 0.0" |
# F PDF and CDF df1 = 5.0; df2 = 9.0 mu = df2 < 2.0 ? 1.0 : df2 / (df2 - 2.0) sigma = df2 < 4.0 ? 1.0 : mu * sqrt(2.0 * (df1 + df2 - 2.0) / (df1 * (df2 - 4.0))) xmin = mu - 3.0 * sigma xmin = xmin < 0 ? 0 : xmin xmax = int(mu + 3.0 * sigma) #Mode of F PDF used ymax = df1 < 3.0 ? 1.0 : 1.1 * f((df1 / 2.0 - 1.0) / (df1 / 2.0 + df1 / df2), df1, df2) set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100 set title "F PDF" set key right box plot [xmin:xmax] [0:ymax] f(x, 5.0, 9.0) title "d_1 = 5, d_2 = 9", \ f(x, 7.0, 6.0) title "d_1 = 7, d_2 = 6" |
set title "F CDF" set key left box plot [xmin:xmax] [0:1.1] cf(x, 5.0, 9.0) title "d_1 = 5, d_2 = 9", \ cf(x, 7.0, 6.0) title "d_1 = 7, d_2 = 6" |
# Gamma PDF and incomplete gamma CDF rho = 1.0; lambda = 1.3 mu = rho / lambda sigma = sqrt(rho) / lambda xmin = mu - 4.0 * sigma xmin = xmin < 0 ? 0 : xmin xmax = mu + 4.0 * sigma ymax = rho < 1.0 ? 2.0 : 1.1 * gmm((rho - 1.0) / lambda, rho, lambda) #Mode of gamma pdf used set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set title "Gamma Γ PDF" set key right array r[7] = [0.5, 1.0, 1.0, 1.3, 2.0, 4.0, 6.0] array l[7] = [1.0, 1.0, 1.3, 1.3, 2.0, 2.0, 2.0] set arrow 1 from 1,1.3 to 0.15,gmm(0.15,r[1],l[1]) set label 1 "ρ < 1, tends to infinity" at 1.1,1.3 left set arrow 2 from 1.15,1.1 to 0.35,gmm(0.35,r[3],l[3]) set label 2 "ρ = 1, finite, nonzero limit" at 1.25,1.1 left set arrow 3 from 1.5,0.9 to 1.0,gmm(1.0,r[5],l[5]) set label 3 "ρ > 1, tends to zero" at 1.6,0.9 left keystr(rho,lambda) = sprintf("ρ = %0.1f, λ = %0.1f", rho, lambda) plot [0:5] [0:1.5] \ for [i = 1:7] rho = r[i], lambda = l[i], gmm(x, rho, lambda) title keystr(rho,lambda) |
unset label 1; unset label 2; unset label 3 unset arrow 1; unset arrow 2; unset arrow 3 set title "incomplete gamma CDF" set key right bottom plot [0:5] [0:1.1] \ for [i = 1:7] rho = r[i], lambda = l[i], cgmm(x, rho, lambda) title keystr(rho,lambda) |
undefine r l # Geometric PDF and CDF p = 0.4 mu = (1.0 - p) / p sigma = sqrt(mu / p) xmin = int(mu - 4.0 * sigma) xmin = xmin < -1 ? -1 : xmin xmin = -1 xmax = int(mu + 4.0 * sigma) ymax = 1.1 * geometric(0, p) #mode of geometric PDF used unset key unset zeroaxis set xrange [xmin : xmax] set yrange [0 : ymax] set xlabel "k" set ylabel "probability density" set ytics 0, ymax / 10, ymax set format x "%2.0f" set format y "%3.2f" set sample (xmax - xmin) + 1 set title "geometric PDF with p = 0.4" plot geometric(x, p) with impulses |
set title "geometric CDF with p = 0.4" set yrange [0 : 1.1] set ytics 0, 0.1, 1.0 plot cgeometric(x, p) with steps |
# Half normal PDF and CDF mu = sqrt2invpi sigma = 1.0 s = sigma*sqrt(1.0 - 2.0/pi) xmin = -0.2 xmax = mu + 4.0 * s ymax = 1.1 * halfnormal(0, sigma) #Mode of half normal PDF used unset key set zeroaxis set xrange [xmin: xmax] set yrange [-0.1: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set parametric set trange [xmin:xmax] set title "half normal PDF, σ = 1.0" set arrow 1 from 0.5,0.13 to 0.0,0.4 set label 1 "Discontinuity achieved by plotting\ntwice with limited parametric ranges" at 0.2,0.1 left 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 |
set title "half normal CDF, σ = 1.0" set yrange [-0.1:1.1] set arrow 1 from 0.45,0.1 to 0.05,0.01 set label 1 "Cusp achieved by plotting twice\nwith limited parametric ranges" at 0.5,0.1 left 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 |
unset label 1 unset arrow 1 unset parametric # Hypergeometric PDF and CPF N = 75; C = 25; d = 10 p = real(C) / N mu = d * p sigma = sqrt(real(N - d) / (N - 1.0) * d * p * (1.0 - p)) xmin = int(mu - 4.0 * sigma) xmin = xmin < -1 ? -1 : xmin xmax = int(mu + 4.0 * sigma) xmax = xmax < d+1 ? d+1 : xmax ymax = 1.1 * hypgeo(int(mu),N,C,d) # approximate mode of hypergeometric PDF used unset key unset zeroaxis set xrange [xmin : xmax] set yrange [0 : ymax] set xlabel "k" set ylabel "probability density" set ytics 0, ymax / 10, ymax set format x "%2.0f" set format y "%3.2f" set sample (xmax - xmin) + 1 set title "hypergeometric PDF with N = 75, C = 25, d = 10" plot hypgeo(x,N,C,d) with impulses |
set yrange [0 : 1.1] set ytics 0, 1.0 / 10.0, 1.1 set title "hypergeometric CDF with N = 75, C = 25, d = 10" plot chypgeo(x,N,C,d) with steps |
# Laplace PDF mu = 0.0; b = 1.0 sigma = sqrt(2.0) * b xmin = mu - 4.0 * sigma xmax = mu + 4.0 * sigma ymax = 1.1 * laplace(mu, mu, b) #Mode of laplace PDF used unset key set zeroaxis set xrange [xmin: xmax] set yrange [0: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100+1 set title "Laplace (or double exponential) PDF with µ = 0, b = 1" set arrow 1 from -0.95,0.5 to -0.1,0.5 set label 1 "Cusp achieved by selecting point\nas part of function samples" at -1.0,0.5 right plot laplace(x, mu, b) |
unset label 1 unset arrow 1 set title "Laplace (or double exponential) CDF with µ = 0, b = 1" set yrange [0: 1.1] plot claplace(x, mu, b) |
# Logistic PDF and CDF a = 0.0; lambda = 2.0 mu = a sigma = pi / (sqrt(3.0) * lambda) xmin = mu - 4.0 * sigma xmax = mu + 4.0 * sigma ymax = 1.1 * logistic(mu, a, lambda) #Mode of logistic PDF used unset key set zeroaxis set xrange [xmin: xmax] set yrange [0: ymax] unset key set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set title "logistic PDF with a = 0, λ = 2" plot logistic(x, a, lambda) |
set title "logistic CDF with a = 0, λ = 2" set yrange [0: 1.1] plot clogistic(x, a, lambda) |
# Lognormal PDF and CDF mu = 1.0; sigma = 0.5 m = exp(mu + 0.5 * sigma**2) s = sqrt(exp(2.0 * mu + sigma**2) * (2.0 * exp(sigma) - 1.0)) xmin = m - 4.0 * s xmin = xmin < 0 ? 0 : xmin xmax = m + 4.0 * s ymax = 1.1 * lognormal(exp(mu - sigma**2), mu, sigma) #Mode of lognormal PDF used unset key set zeroaxis set xrange [xmin: xmax] set yrange [0: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.2f" set format y "%.2f" set sample 100 set title "lognormal PDF with µ = 1.0, σ = 0.5" plot lognormal(x, mu, sigma) |
set title "lognormal CDF with µ = 1.0, σ = 0.5" set yrange [0: 1.1] plot clognormal(x, mu, sigma) |
# Maxwell PDF a = 0.5 mu = 2.0 / sqrt(pi) / a sigma = sqrt(3.0 - 8.0/pi) / a xmin = int(mu - 3.0 * sigma) xmin = xmin < 0 ? 0 : xmin xmax = int(mu + 3.0 * sigma) a = 1.5 ymax = 1.1 * maxwell(1.0 / a, a) + 0.5 #Mode of maxwell PDF used ymax = int(ymax + 0.5) set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set title "Maxwell PDF" set key right top box plot [xmin:xmax] [0:ymax] maxwell(x, 1.5) title "a = 1.5", \ maxwell(x, 1.0) title "a = 1.0", \ maxwell(x, 0.5) title "a = 0.5" |
set title "Maxwell CDF" set key right bottom box plot [xmin:xmax] [0:1.1] cmaxwell(x, 1.5) title "a = 1.5", \ cmaxwell(x, 1.0) title "a = 1.0", \ cmaxwell(x, 0.5) title "a = 0.5" |
# Negative binomial PDF and CDF r = 8; p = 0.4 mu = r * (1.0 - p) / p sigma = sqrt(mu / p) xmin = int(mu - 4.0 * sigma) xmin = xmin < 0 ? 0 : xmin xmax = int(mu + 4.0 * sigma) ymax = 1.1 * negbin(int(mu - (1.0-p)/p), r, p) #mode of gamma PDF used unset key unset zeroaxis set xrange [xmin-1 : xmax] set yrange [0 : ymax] set xlabel "k" set ylabel "probability density" set ytics 0, ymax / 10, ymax set format x "%2.0f" set format y "%3.2f" set sample (xmax - xmin+1) + 1 set title "negative binomial (or Pascal or Polya) PDF with r = 8, p = 0.4" plot negbin(x, r, p) with impulses |
set yrange [0 : 1.1] set ytics 0, 0.1, 1.0 set title "negative binomial (or Pascal or Polya) CDF with r = 8, p = 0.4" plot cnegbin(x, r, p) with steps |
# Negative exponential PDF and CDF lambda = 2.0 mu = 1.0 / lambda sigma = 1.0 / lambda xmax = mu + 4.0 * sigma ymax = lambda #No mode unset key set zeroaxis set xrange [0: xmax] set yrange [0: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.2f" set format y "%.1f" set sample 100 set title "negative exponential (or exponential) PDF with λ = 2.0" plot nexp(x, lambda) |
set title "negative exponential (or exponential) CDF with λ = 2.0" set yrange [0: 1.1] plot cnexp(x, lambda) |
# Normal PDF and CDF mu = 0.0; sigma = 1.0 xmin = mu - 4.0 * sigma xmax = mu + 4.0 * sigma mu = 2.0; sigma = 0.5 ymax = 1.1 * normal(mu, mu, sigma) #Mode of normal PDF used set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 100 set title "normal (also called Gauss or bell-curved) PDF" set key left top box plot [xmin:xmax] [0:ymax] normal(x, 0, 1.0) title "µ = 0, σ = 1.0", \ normal(x, 2, 0.5) title "µ = 2, σ = 0.5", \ normal(x, 1, 2.0) title "µ = 1, σ = 2.0" |
set title "normal (also called Gauss or bell-curved) CDF" set key left top box plot [xmin:xmax] [0:1.1] mu = 0, sigma = 1.0, cnormal(x, mu, sigma) title "µ = 0, σ = 1.0", \ mu = 2, sigma = 0.5, cnormal(x, mu, sigma) title "µ = 2, σ = 0.5", \ mu = 1, sigma = 2.0, cnormal(x, mu, sigma) title "µ = 1, σ = 2.0" |
# Pareto PDF and CDF a = 1.0; b = 3.0 mu = a * b / (b - 1.0) sigma = a * sqrt(b) / (sqrt(b - 2.0) * (b - 1.0)) xmin = mu - 4.0 * sigma xmin = xmin < 0 ? 0 : xmin xmax = int(mu + 4.0 * sigma) ymax = 1.1 * pareto(a, a, b) #mode of pareto PDF used ymin = -0.1 * pareto(a, a, b) unset key set zeroaxis set xrange [xmin: xmax] set yrange [ymin: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.1f" set sample 200+1 set title "Pareto PDF with a = 1, b = 3" # Discontinuity at a set parametric set trange [0:1-eps] x1(t) = -1 + 2*t x2(t) = 1 + 3*t set arrow 1 from 1.75,0.8 to 1.0,0.8 set arrow 2 from 1.0,0.0 to 1.0,3.0 nohead lt 0 set label 1 "Discontinuity achieved by plotting twice\nwith affine mapped parametric ranges" at 1.8,0.8 left plot x1(t), pareto(x1(t), a, b) ls 1, x2(t), pareto(x2(t), a, b) ls 1 |
unset arrow 2 set title "Pareto CDF with a = 1, b = 3" unset parametric set yrange [-0.1: 1.1] set arrow 1 from 1.45,0.1 to 1.05,0.01 set label 1 "Cusp achieved by selecting point\nas part of function samples" at 1.5,0.1 left plot cpareto(x, a, b) |
unset label 1 unset arrow 1 # Poisson PDF and CDF mu = 4.0 sigma = sqrt(mu) xmin = int(mu - 4.0 * sigma) xmin = xmin < -1 ? -1 : xmin xmax = int(mu + 4.0 * sigma) ymax = 1.1 * poisson(mu, mu) #mode of poisson PDF used unset key set zeroaxis set xrange [xmin : xmax] set yrange [0 : ymax] set xlabel "k" set ylabel "probability density" set ytics 0, ymax / 10, ymax set format x "%2.0f" set format y "%3.2f" set sample (xmax - xmin) + 1 set title "Poisson PDF with µ = 4.0" plot poisson(x, mu) with impulses |
set yrange [-0.1 : 1.1] set ytics -0.1, 0.1, 1.1 set title "Poisson CDF with µ = 4.0" plot cpoisson(x, mu) with steps |
# Rayleigh PDF and CDF lambda = 2.0 mu = 0.5 * sqrt(pi / lambda) sigma = sqrt((1.0 - pi / 4.0) / lambda) xmax = mu + 4.0 * sigma ymax = 1.1 * rayleigh(1.0 / sqrt(2.0 * lambda), lambda) #Mode of rayleigh PDF used unset key set zeroaxis set xrange [0: xmax] set yrange [0: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.2f" set format y "%.1f" set sample 100 set title "Rayleigh PDF with λ = 2.0" plot rayleigh(x, lambda) |
set title "Rayleigh CDF with λ = 2.0" set yrange [0: 1.1] plot crayleigh(x, lambda) |
# Sine PDF and CDF a = 3.2; f = 2.6 mu = a / 2.0 sigma = sqrt(a * a / 3.0 * (1.0 - 3.0 / (2.0 * n * n * pi * pi)) - mu * mu) xmin = 0.0 xmax = a - eps a = 2; f = 1.0 ymax = 1.1 * 2.0 / a #Mode of sine PDF used set zeroaxis set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.2f" set format y "%.1f" set sample 250 set title "sine PDF" set key bottom outside keystr(a, f) = sprintf("a = %0.1f, f = %0.1f", a, f) array aa[4] = [2.0, 2.0, 3.25, 2.75] array fa[4] = [1.0, 3.0, 2.6, 0.0] a1 = 2.0; a2 = 3.25; a3 = 2.75 f1 = 1.0; f2 = 3.0; f3 = 2.6; f4 = 0.0 plot [xmin:xmax] [0:ymax] a = a1, f = f1, sine(x, f, a) title keystr(a, f), \ a = a1, f = f2, sine(x, f, a) title keystr(a, f), \ a = a2, f = f3, sine(x, f, a) title keystr(a, f), \ a = a3, f = f4, sine(x, f, a) title keystr(a, f) with steps |
set title "sine CDF" set key top left plot [xmin:xmax] [0:1.1] a = a1, f = f1, csine(x, f, a) title keystr(a, f), \ a = a1, f = f2, csine(x, f, a) title keystr(a, f), \ a = a2, f = f3, csine(x, f, a) title keystr(a, f), \ a = a3, f = f4, csine(x, f, a) title keystr(a, f) with steps |
# t PDF and CDF nu = 20 mu = 0 sigma = nu > 2 ? sqrt(nu / (nu - 2.0)) : 1.0 xmin = mu - 4.0 * sigma xmax = mu + 4.0 * sigma ymax = 1.1 * t(mu, nu) #Mode of t PDF used set key inside center left title "degrees of freedom" set zeroaxis set xrange [xmin: xmax] set yrange [0: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 100 set title "t PDF (and Gaussian limit)" ks(nu) = sprintf("ν = %d", nu) array nu[6] = [1, 2, 3, 4, 10, 20] plot \ for [i=1:6] t(x, nu[i]) ti ks(nu[i]), \ normal(x, 0, 1) ti "normal" |
set title "t CDF (and Gaussian limit)" set yrange [0: 1.1] plot \ for [i=1:6] ct(x, nu[i]) ti ks(nu[i]), \ cnormal(x, 0, 1) ti "normal" |
undefine nu # Thanks to efrank@upenn5.hep.upenn.edu for telling us about this # triangular PDF and CDF m = 3.0 g = 2.0 mu = m sigma = g/sqrt(6.0) xmin = m - 1.1*g xmax = m + 1.1*g ymax = 1.1 * triangular(m, m, g) #Mode of triangular PDF used ymin = -ymax/11.0; unset key set zeroaxis set xrange [xmin: xmax] set yrange [ymin: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.1f" set format y "%.2f" set sample 50*1.1*g+1 set title "triangular PDF with m = 3.0, g = 2.0" plot triangular(x, m, g) |
set title "triangular CDF with m = 3.0, g = 2.0" set yrange [-0.1: 1.1] plot ctriangular(x, m, g) |
# Uniform PDF and CDF a = -2.0; b= 2.0 mu = (a + b) / 2.0 sigma = (b - a) / sqrt(12.0) xmin = a - 0.1*(b - a) xmax = b + 0.1*(b - a) ymax = 1.1 * uniform(mu, a, b) #No mode ymin = -0.1 * uniform(mu, a, b) unset key set zeroaxis set xrange [xmin: xmax] set yrange [ymin: ymax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%.2f" set format y "%.2f" set sample 120+1 set title "uniform PDF with a = -2.0, b = 2.0" plot uniform(x, a, b) with steps |
set title "uniform CDF with a = -2.0, b = 2.0" set yrange [-0.1 : 1.1] plot cuniform(x, a, b) |
# Weibull PDF and CDF lambda = 1.0/5; a = 1.0 mu = 1.0 / lambda * gamma(1.0 / a) / a sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2)) xmin = mu - 4.0 * sigma xmin = xmin < 0 ? 0 : xmin #Mode of weibull PDF used ymax = 1.8 * (a >= 1.0 ? weibull(((a - 1.0) / a)**(1.0 / a) / lambda, a, lambda) : 2.0) lambda = 1.0/15; a = 10.0 mu = 1.0 / lambda * gamma(1.0 / a) / a sigma = sqrt(lambda**(-2.0) * (2.0 * gamma(2.0 / a) / a - (gamma(1.0 / a) / a)**2)) xmax = int(mu + 4.0 * sigma) set key on title "" inside top right set zeroaxis set grid set xrange [xmin : xmax] set xlabel "x" set ylabel "probability density" set xtics autofreq set ytics autofreq set format x "%g" set format y "%g" set sample 100 set title "Weibull PDF" ks(a,lambda) = sprintf("λ = 1/%g, a = %0.1f", 1.0/lambda, a) a1 = 0.5; a2 = 1.0; a3 = 2.0; a4 = 10.0 lambda1 = 1.0/5; lambda2 = 1.0/15 set arrow 1 from 3.8,0.27 to 0.5,weibull(0.5,a1,lambda1) set label 1 "a < 1, rate decreasing over time" at 4,0.27 left set arrow 2 from 8,0.19 to 6.4,weibull(6.4,a3,lambda1) set arrow 3 from 10.5,0.19 to 13,weibull(13,a4,lambda2) set label 2 "a > 1, rate increasing over time" at 9,0.2 center plot [] [0:ymax] lambda = lambda1, a = a1, weibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda1, a = a2, weibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda1, a = a3, weibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda2, a = a4, weibull(x, a, lambda) ti ks(a, lambda) |
unset label 1; unset label 2 unset arrow 1; unset arrow 2; unset arrow 3 set key at 9,0.4 center set title "Weibull CDF" plot [] [0:1.1] lambda = lambda1, a = a1, cweibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda1, a = a2, cweibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda1, a = a3, cweibull(x, a, lambda) ti ks(a, lambda), \ lambda = lambda2, a = a4, cweibull(x, a, lambda) ti ks(a, lambda) |
reset set encoding save_encoding |