---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
 Examples of distributions

# References for the most common probability distributions are available here:
help("distribution")
#
# I toss up a coin 10 times. What is the law of distribution of the random variable:
# N = "number of heads" ?
# Experimental study:
n = 10000; da = NULL; for (i in 1:n) da[i]=sum(RUNIF(10, 0,1))
hist(da, seq(-0.5,10.5, 1),probability=TRUE,col="grey",main="",cex.axis=0.9)
BOX(); abovex("heads in 10 tosses")
statistics(da)
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.000   4.000   5.000   4.983   6.000  10.000

   
# The theoretical study (for choose see):
N=10; da = NULL; for (i in 0:N) da = c(da,rep(i,choose(10,i)))
hist(da, seq(-0.5,10.5, 1),probability=TRUE,col="grey",main="",cex.axis=0.9)
BOX(); abovex("heads in 10 tosses")
summary(da)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#      0       4       5       5       6      10
# As index of the data dispersion I can take 3^quartile-1^quartile, ie:
Percentile(da,75)-Percentile(da,25)
[1] 2
# or the standard deviation of "da" (see here)
Sd(da)
# 1.581139
# Using Var (= √Sd) I see that Sd(da) √2.5
Var(da)
# 2.5                   It is demonstrable that Var is √N/2
# I recall that these are weighted values. I can proceed also in this way (see here):
j=0:10; Wmean(j,choose(10,j)); WVar(j,choose(10,j))
#            5                    2.5

# This is an example of binomial distribution  Pr(N=k) = C(n,k)*p^k*(1p)^(n-k)
# when p = 1/2. In general: 
binom = function(m) choose(N,m)*p^m*(1-p)^(N-m)

# An example with p ≠ 1/2.
# An automatic oven produces on average 1 biscuit burnt every 8; biscuits are then
# automatically mixed and packaged in packs of 6. What is the distribution law of the
# random variable  N = "n. of burnt cookies in a pack" ?
# It is a discrete probability distribution.
N = 6; p = 1/8; x = 0:N; BC = binom(0:N)
Plane(0,N, 0,max(BC) ); POINT(x,BC, "blue")
# For a better view we add these commands and obtain the following graph, on the left:
polyl(x,BC, "blue"); segm(x,0, x,BC, "blue")
   
Wmean(x,BC); WVar(x,BC)
# 0.75      0.65625
# I can prove that I obtain the same valued with:
N*p; N*p*(1-p)
# 0.75  0.65625
 
# To understand this it is useful to interpret BC as X1+…+X6 where Xi is 1 with proba-
# bility 1/8, is 0 otherwise.  With this simulation we obtain the graph on the right.
n=10000; da=NULL; N=6
for(i in 1:n) {X=0; for(j in 1:N) if(runif(1, 0,1)<1/8) X=X+1; da[i]=X} 
hist(da, seq(-0.5,N+1/2, 1),probability=TRUE,col="grey",main="",cex.axis=0.9)
BOX(); abovex("burnt cookies in a pack")
 
# It's clear why the mean is N*p: if every Xi is 1 with probability p X1+… has mean N*p
 
# NOTE: y = dbinom(x, N, p) is equivalent to  x=0:N; ; binom = ...; y = binom(0:N)
# See here for an example of the multinomial distribution.
 
# When N is large (as in the following example) I can approximate the binomial
# distribution with the gaussian distribution, that we recall dawn.
# The following graph, on the left, is for N=20:
          
p=0.2; N=20
x=0:N; y=binom(0:N); Plane(0,N, 0,max(y))
POINT(x,y, "blue"); polyl(x,y, "blue"); segm(x,0, x,y, "blue")
m=N*p; v=N*p*(1-p); s=sqrt(v)
# 4      3.2        1.788854
 
# The distribution rule called normal or gaussian has the following density function f;
# it depends only on 2 parameters, the mean m and σ, called the standard deviation of f.
 
           
#      (σ is the Greek letter "sigma", which corresponds to our "s")
 
# In R the normal density of mean m and standard deviation s is called
# dnorm(x, mean=m, sd=s). How to obtain the graph above on the right:
 
z = function(x) dnorm(x, mean=m, sd=s)
x=0:N; y=binom(0:N); Plane(0,N, 0,max(y))
POINT(x,y, "blue")
graph2(z, -20,20, "brown")
 
# One of the (not many) examples of gaussian phenomenon is the height of the people of
# a particular sex, age, place of origin.  Let's use the data already considered here
# (medical examination for a call-up - Italy, 1987):
 
da = read.table("http://macosa.dima.unige.it/R/callup.txt",skip=8, sep=";",header=TRUE)
str(da)
#  'data.frame':   4170 obs. of  4 variables:
#   $ chest : int  91 86 89 88 92 95 91 83 87 88 ...
#   $ weight: int  62 64 64 68 70 77 65 62 55 60 ...
#   $ region: int  15 15 15 15 15 15 15 15 15 15 ...
#   $ height: int  73 70 74 79 66 69 66 79 75 65 ...
# We add 1/2 to truncated data, so we do not have to add 1/2 to the average values
da$chest = da$chest+1/2; da$weight = da$weight+1/2; da$height = da$height+1/2
H = da$height; min(H); max(H)
#              47.5    96.5
m=mean(H); s=Sd(H)
z = function(x) dnorm(x, mean=m, sd=s)
BF=3;HF=2.1
graph2F(z, 47,97, "brown")
hist(H,n=20,right=FALSE,probability=TRUE,col="grey",add=TRUE)
graph2(z,47,97, "brown"); abovex("data")
POINT(m+s,z(m+s),"blue"); POINT(m-s,z(m-s),"blue")
# The data fit in the graph: see below on the left.
# I can simulate the phenomenon with rnorm(n=N, mean=m, sd=s) that generates N random
# data with normal distribution; obviously also the simulation fit (graph on the right):
# 
H1 = rnorm(n=4000, mean=m, sd=s)
graph2F(z, 47,97, "brown")
hist(H1,n=20,right=FALSE,probability=TRUE,col="grey",add=TRUE)
graph2(z,47,97, "brown"); abovex("random")
POINT(m+s,z(m+s),"blue"); POINT(m-s,z(m-s),"blue")

   

# We pointed out (in blue) in the graphs also m ± σ. One can prove that these points
# are flex points. In general, whatever the gaussian, we have:
integral(z,   m-s,m+s)        #  0.6826895
integral(z, m-s*2,m+s*2)      #  0.9544997
integral(z, m-s*3,m+s*3)      #  0.9973002
# and obviously:
integral(z,  -Inf,Inf)        #  1
 
# The variable "height" is virtually continuous: unlike the number of coins or biscuits
# which only assumed integer values, I can think that a length is a real number.
# Increasing the number of intervals the upper edge of the histogram tends to a curve;
# instead of the sum of the areas of the columns of a histogram between x=a and x=b, I
# have the integral of a function between x=a and x=b.
# In the discrete case I represent the frequency density (ie the unitary percentage
# frequency), in the continuous one I represent the density function.
 
                      ---------- ---------- ---------- ----------
 
# We note that x1k, x2k, x3k, … (if k ≠ 1) have a different kind of distribution of x1,
# x2, x3, …   An example:  x with gaussian distribution and x^3.
 
x <- rnorm(n=10000,mean=5,sd=1); x3=x^3
Plane(min(x),max(x),0,max(hist(x,plot=FALSE)$density))
hist(x,col="grey",probability=TRUE,add=TRUE)
statistics(x)
# median = 5.005     mean = 5.002
Plane(min(x3),max(x3),0,max(hist(x3,plot=FALSE)$density))
hist(x3,col="grey",probability=TRUE,add=TRUE)
statistics(x3)
# median = 125.400   mean = 140.200
        x    x^3
 
# x -> x^3 is a growing function; so the median of x^3 is the cube of the median of x;
# this does not happen for the mean!
 
                      ---------- ---------- ---------- ----------
 
# Generalizing what we just saw (here), we have:
# if I sum up N random variables that have the same distribution, I obtain a random
# variable that (as N -> ∞) tends toward a gaussian distribution. An example:
tot = 2000; x = NULL               #                N = 1
for (i in 1:tot) x[i] = sqrt(runif(1))+3*runif(1)^2
dev.new(); hist(x, probability=TRUE, col="grey",main="",cex.axis=0.9)
BOX(); abovex("n=1")
#                                                   N = 20
for (i in 1:tot) {x[i]=0; for (j in 1:20) x[i]=x[i]+sqrt(runif(1))+3*runif(1)^2}
dev.new(); hist(x, probability=TRUE, col="grey",main="",cex.axis=0.9)
BOX(); abovex("n=20")
    
 
# That's the reason why the gaussian distribution is important:  there are a few
# gaussian phenomena, but if I calculate the mean of N values with the same
# distribution (ie if I calculate their sum divided for N), this mean has a gaussian
# distribution.  This fact goes under the name of central limit theorem:
 
# If U1,…,Un are random indipendent variables with the same law of distribution, and if
# m and σ0 are its mean and standard deviation, when n grows the variable (U1+  +Un)/n
# tends to have a normal distribution with mean m and standard deviation σ0/√n.
 
# Let's consider as an example of non-Gaussian distribution the weights of the
# population of twenty years considered above.
da = read.table("http://macosa.dima.unige.it/R/callup.txt",skip=8, sep=";",header=TRUE)
da$chest = da$chest+1/2; da$weight <- da$weight+1/2; da$height <- da$height+1/2
W = da$weight; min(W); max(W)
#              48.5    126.5
BF=3.5; HF=2.5
max(hist(W, probability=TRUE, col="grey",main="",cex.axis=0.9)$density)
#   0.03856115
Plane(48,127,0,0.04)
hist(W, probability=TRUE, col="grey",cex.axis=0.9,add=TRUE)
abovex("weight - call-up - Italy, 1987"); POINT(mean(W),0,"blue")
 
      
length(W); Sd(W)/sqrt(length(W))
#  4170     0.1630721
E=3*Sd(W)/sqrt(length(W)); mean(W)-E; mean(W)+E
#           70.48752   71.46595
# The mean weight of the twenty years italian men was 71 ± 0.5 kg
 
# But the mean and the standard deviation that we considered are not the mean and the
# standard deviation of the population of all the twenty years italian men (in 1987),
# they are those of the sample. It can be proved that, while the statistical mean is a
# good approximation of the theoretical one, I have to consider (insted of the
# theoretic standard deviation, Sd) the experimental standard deviation sd, where n is
# replaced by n-1:
#                 Sd  is  √( ( (X1-M)^2 + … + (Xn-M)^2 ) / n )
#                 sd  is  √( ( (X1-M)^2 + … + (Xn-M)^2 ) / (n-1) )  [ > Sd]
# So, (U1+  +Un)/n tends to have a normal distribution which have the experimental
# mean as a mean and sd/√n as a standard deviation. I can calculare this value with
# SdM (standard deviation of the mean).  See here
# When n is large the difference between the two values is negligible:
 
Sd(W)/sqrt(4170); sd(W)/sqrt(4170); SdM(W)
#  0.1630721       0.1630916        0.1630916
# It is due to the replacement of √n with √(n-1):
n=4170; sqrt(n/(n-1)); n=20; sqrt(n/(n-1)); n=10; sqrt(n/(n-1))
#       1.00012              1.025978             1.054093
 
# Therefore, more correct calculations in our case would be:
E=3*SdM(W); mean(W)-E; mean(W)+E
#           70.48746   71.46601
 
                      ---------- ---------- ---------- ----------
 
# In the discrete case the probability that U is between a and b is Pr(U=a)+…+Pr(U=b);
# in the continuous case the probability that the density function U is between a and b
# (where a and b may be infinite) is [a,b] U(x) dx.  In the discrete case if U has values
# in {v1,v2,…} the mean is v1·Pr(U=v1) + v2·Pr(U=v2) + …;  in the continuous case if U
# varies in I it is I x·U(x) dx;  so, if m is the mean, the variance is I (x-m)²·U(x) dx

# If I sum up only two random variables uniformly distributed on [0,1] I have (just
# like we saw here for the launch of two balanced dice):
HF=2; BF=4
PLANE(0,2, 0,1)
n = 1e6; U1 = runif(n,0,1); U2 = runif(n,0,1)
hist(U1+U2, probability=TRUE, col="grey90", add=TRUE, nclass=40)   # or:  n=40
# the density function (from x=0 to x=2) is:
f = function(x) -abs(x-1)+1
graph1(f, 0,2, "blue"); abovex("runif+runif")
# and (as it should be) I have:
integral(f, 0,2)
# 1
# and the mean:
g = function(x) x*f(x); integral(g, 0,2)
# 1
# The graph of the probability that U1+U2 is less than x (see here):
Gintegra(f,0,2,"black")
    
 
# This function (the function with black graph) is the cumulative distribution function
# of U1+U2.
 
# In the case of the gaussian distribution I have:
 
z = function(x) dnorm(x, mean=0, sd=1)
HF=2; BF=3
graphF(z, -4,4, "blue")
# On the right the graph of the cumulative distribution function:
Plane(-4,4, 0,1)
Gintegra(z,-10,10, "brown")
# I should have drawn the graph of ∫(-∞,x] z, but
# in Gintegra(f,A,B,col) I must use A and B finite. However ∫(-10,x] z is practically the
# same; indeed:
integral(z, -Inf,-10)
# 7.619853e-24
   
 
                      ---------- ---------- ---------- ----------
 
# The examples of simulation considered here employ a procedure that is based on
# central limit theorem and the methods for studying the binomial distribution
# considered at the beginning of this document, but merely compute the mean and the
# standard deviation  [ in this case N=1 and Var=p*(1-p) ]
 
# Let's remember the first example:
# the probability that throwing a pair of balanced dice their difference is 1
# (we - really - know it is 5/18 = 0.2777…)
 
Event <- function() abs( RUNIF(1, 1,6)-RUNIF(1, 1,6) ) == 1
 
PR(1e4)
# 27.9 +/- 1.34559 
PR(1e5)
# 27.737 +/- 0.4247283 
PR(1e6)
# 27.8421 +/- 0.1344666
 
# As n -> ∞ the standard deviation σ/√n -> 0 as 1/√n: the "error" of PR(N) is divided
# by 10 if N multiplies by 100.
 
# It is basically for this reason that (while analyzing statistical data) if the number
# of data is multiplied by 100 I can round the average with a further figure.
# If I have at least 10 data which are rounded to the thents I can round the mean to
# the cents.  If I have at least 1000 data which are rounded to the thents I can round
# the mean to the thousandths.  If I have at least 10 data which are rounded to the
# units I can round the mean to the tens.   Verification:
 
# Some data rounded to the units:
data=runif(11, 50,100);     c( mean(data),  mean(round(data)) )
# 77.32430  77.27273
data=runif(101, 50,100);    c( mean(data),  mean(round(data)) )
# 75.86425  75.83168
data=runif(1001, 50,100);   c( mean(data),  mean(round(data)) )
# 74.99541  74.99800
data=runif(10001, 50,100);  c( mean(data),  mean(round(data)) )
# 75.18941 75.18628
data=runif(100001, 50,100); c( mean(data),  mean(round(data)) )
# 74.94769  74.94847
 
# Data must be rounded!!!  See here.  See here too.

Other distributions

Other examples of use