#
# $Id: random.dem,v 1.12.2.2 2008/09/23 23:11:31 sfeam Exp $
#
# random.dem
#
# Lattice test for random numbers;
# If you can see any patterns in this plot, the random number generator
# is not very good.
#
# Copyright (c) 1991, Jos van der Woude, jvdwoude@hut.nl

# History:
#	-  6. 6. 2006 ds: added univariate and multivariate normal example
#	- 10. 5. 2006 ds: added univariate and multivariate normal example
#	-  ?. ?  1991 jvdw: 1st version

unset key
set xrange [0: 1]
set yrange [0: 1]
set zrange [0: 1]
set title "Lattice test for random numbers"
set xlabel "rand(n) ->"
set ylabel "rand(n + 1) ->"
set zlabel "rand(n + 2) ->"
set format x "%3.2f"
set format y "%3.2f"
set format z "%3.2f"
set tics
set sample 1000
set style function dots
set parametric
plot rand(0), rand(0)
pause -1 "Hit return to continue"
print "3D plot ahead, one moment please ..."
set sample 50
splot rand(0), rand(0), rand(0)
pause -1 "Hit return to continue"

print ""
print "Multivariate normal distribution"
print ""
print "The surface plot shows a two variable multivariate probability"
print "density function.  On the x-y plane are some samples of the random"
print "vector and a contour plot illustrating the correlation, which in"
print "this case is zero, i.e. a circle.  (Easier to view in map mode.)"
print ""
nsamp = 50
unset xlabel
unset ylabel
unset zlabel
set parametric
# A somewhat inelegant way of generating N random data points.  A future
# non-pressing plot-command feature addition may address this issue.
set samples nsamp
set format "%8.5g"
set table "random.tmp"
plot invnorm(rand(0)),invnorm(rand(0))
unset table
unset format
#
tstring(n) = sprintf("%d random samples from a 2D Gaussian PDF with\nunit variance, zero mean and no dependence", n)
set title tstring(nsamp)
unset key
set hidden3d
set contour
set view 68, 28, 1, 1
set cntrparam levels discrete 0.1
unset clabel
set xrange [-3:3]
set yrange [-3:3]
set zrange [-0.2:0.2]
set ztics 0,0.05
set urange [-3:3]
set vrange [-3:3]
set ticslevel 0
set isosamples 30
splot u,v,( 1/(2*pi) * exp(-0.5 * (u**2 + v**2)) ) with line lc rgb "black", \
   "random.tmp" using 1:2:(-0.2) with points pointtype 7 lc rgb "black"
pause -1 "Hit return to continue"
unset contour

load "stat.inc"

print ""
print "Simple Monte Carlo simulation"
print ""
print "The first curve is a histogram where the binned frequency of occurence"
print "of a pseudo random variable distributed according to the normal"
print "(Gaussian) law is scaled such that the histogram converges to the"
print "normal probability density function with increasing number of samples"
print "used in the Monte Carlo simulation.  The second curve is the normal"
print "probability density function with unit variance and zero mean."
print ""
nsamp = 5000
binwidth = 20
xlow = -3.0
xhigh = 3.0
scale = (binwidth/(xhigh-xlow))
# A somewhat inelegant way of generating N random data points.  A future
# non-pressing plot-command feature addition may address this issue.
set parametric
set samples nsamp
set format "%8.5g"
set table "random.tmp"
plot invnorm(rand(0)),(1.0*scale/nsamp)
unset table
unset format
#
unset parametric
set samples 200
tstring(n) = sprintf("Histogram of %d random samples from a univariate\nGaussian PDF with unit variance and zero mean", n)
set title tstring(nsamp)
set key
set grid
set xrange [-3:3]
set yrange [0:0.45]
bin(x) = (1.0/scale)*floor(x*scale)
plot "random.tmp" using (bin($1)):2 smooth frequency with steps \
                title "scaled bin frequency", \
                normal(x,0,1) with lines title "Gaussian p.d.f."
pause -1 "Hit return to continue"

print ""
print "Another Monte Carlo simulation"
print ""
print "This is similar to the previous simulation but uses multivariate"
print "zero mean, unit variance normal data by computing the distance "
print "each point is from the origin.  That distribution is known to fit"
print "the Maxwell probability law, as shown."
print ""
reset
# A somewhat inelegant way of generating N random data points.  A future
# non-pressing plot-command feature addition may address this issue.
nsamp = 3000
set parametric
set samples nsamp
set isosamples 2,2 # Smallest possible
set format "%8.5g"
set table "random.tmp"
splot invnorm(rand(0)),invnorm(rand(0)),invnorm(rand(0))
unset table
unset format
#
oneplot = 1
#
if (oneplot) set multiplot layout 1,2
#
unset key
rlow = -4.0
rhigh = 4.0
set xrange [rlow:rhigh]; set yrange [rlow:rhigh]; set zrange [rlow:rhigh]
set xtics axis nomirror; set ytics axis nomirror; set ztics axis nomirror;
set border 0
set xyplane at 0
set xzeroaxis lt -1
set yzeroaxis lt -1
set zzeroaxis lt -1
set view 68, 28, 1.4, 0.9
tstring(n) = sprintf("Gaussian 3D cloud of %d random samples\n", n)
set title tstring(nsamp) offset graph 0.15, graph -0.33
splot "random.tmp" every :::::0 with dots
if (!oneplot) pause -1 "Hit return to continue"

unset parametric
unset xzeroaxis; unset yzeroaxis;
set border
set grid
set samples 200
set size 0.47,0.72
set origin 0.44,0.18
tstring(n) = sprintf("Histogram of distance from origin of\n%d multivariate unit variance samples", n)
set title tstring(nsamp) offset graph 0, graph 0.15
set key bmargin right vertical
xlow = 0.0
xhigh = 4.5
binwidth = 20
scale = (binwidth/(xhigh-xlow))
set xrange [0:xhigh]
set yrange [0:0.65]
bin(x) = (1.0/scale)*floor(x*scale)
plot "random.tmp" using (bin(sqrt($1**2+$2**2+$3**2))):(1.0*scale/nsamp) every :::::0 smooth frequency with steps \
                title "scaled bin frequency", \
                maxwell(x, 1/sqrt(2)) with lines title "Maxwell p.d.f."
#
if (oneplot) unset multiplot
#
pause -1 "Hit return to continue"

unset key
unset grid
set xrange [-8:8]
set yrange [-8:8]
set size ratio 1.0
set zeroaxis
set border 3
set tics out scale 0.5
set xtics 1.0 border rangelimited nomirror
set ytics 1.0 border rangelimited nomirror
set tics scale 0.5
set format xy "%.0f"

set title 'Example of range-limited axes and tics'

A = pi/9.
plot 'random.tmp' using (1.5 + $1*cos(A)-2.*$2*sin(A)):(1.0 + $1*sin(A)+2.*$2*cos(A)) with dots

pause -1 "Hit return to continue"