# 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 |