#----------------------------------------------------------------------------
# pdf of sum of two variables following asymmetrics triangular distributions,
# based on http://www.roma1.infn.it/~dagos/R/triang.R,
# in order to reproduce the example
# given in http://arxiv.org/abs/physics/0403086, Fig. 1
#
# GdA 30/1/2016
#-----------------------------------------------------------------------------

source("triang.R")

x1 <- rtriang(1000000, 0.5, 1.5, 0.5)
x2 <- rtriang(1000000, 0.5, 1.5, 0.5)

s <- x1 + x2

salva = TRUE

if(salva) png(filename="somma_triangolari_asimmetriche_sampling.png",
              height=600, width=1000, bg="white")

par(mfrow = c(3,1))

h1 <- hist(x1, nc=200, xlim=c(-2,2), col='cyan')
h2 <- hist(x2, nc=200, xlim=c(-2,2), col='cyan')
hs <- hist(s, nc=200, xlim=c(-2,2), col='magenta')
# (details of the histograms are used below)

par(mfrow = c(1,1))

if(salva) dev.off()

# means and standard deviations
m1 = mean(x1)
s1 = sd(x1)
m2 = mean(x2)
s2 = sd(x2)
ms = mean(s)
ss = sd(s)

cat( sprintf("\nMeans: x1: %.3f;  x2:  %.3f;  sum: %.3f\n", m1, m2, ms))
cat( sprintf("StDev: x1: %.3f;  x2:  %.3f;  sum: %.3f\n\n", s1, s2, ss))

# medians 
med1 = median(x1)
med2 = median(x2)
meds = median(s)

cat( sprintf("Medians: x1: %.3f;  x2:  %.3f;  sum: %.3f\n\n",
             med1, med2, meds))

# very 'brutal' extimate of the modes from histograms
mod1 <- h1$mids[ which( h1$counts==max(h1$counts) ) ]
mod2 <- h2$mids[ which( h2$counts==max(h2$counts) ) ]
mods <- hs$mids[ which( hs$counts==max(hs$counts) ) ]

cat( sprintf("Hist. modes (att!): x1: %.3f;  x2:  %.3f;  sum: %.3f\n\n",
             mod1, mod2, mods))