#------------------------------------------
# Problem proposed 14 feb 2025
# 
# GdA, feb 2025
#-----------------------------------------

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

nevts = 100000 

#------------------------------------------------------------
# lambda:  Normal(10000, 100)
mu = 10000; sigma = 100
lambda <- rnorm(nevts, mu, sigma) 

hist(lambda, nc=100, prob=TRUE, col='cyan')
cat(sprintf('\n lambda:   Mean = %.0f,  std = %.0f\n',
            mean(lambda), sd(lambda) ))
print(summary(lambda))
pausa()

#-----------------------------------------------------------
# n: Poisson(lambda)
n <- rpois(nevts, lambda)
hist(n, nc=100, prob=TRUE, col='cyan')
cat(sprintf('\n n:   Mean = %.0f,  std = %.0f\n',
            mean(n), sd(n) ))
print(summary(n))    
pausa() 

#-----------------------------------------------------------
# x: Binomial(n, p) with certain p (-> 'p0')       
p0 <- 0.95          
x <- rbinom(nevts, n, p0)
hist(x, nc=100, prob=TRUE, col='cyan')
cat(sprintf('\n x:   Mean = %.0f,  std = %.0f\n',
            mean(x), sd(x) ))
print(summary(x))    

#-----------------------------------------------------------
# Variant with uncertain p, described by a Beta  
# having meean 0.95 e sd 0.05
mu.p    = 0.95
sigma.p = 0.05

# Calculate the Beta parameters r nd s 
# (in the literature also indicated by 'alpha' e 'beta'))
r <- mu.p * ( mu.p*(1-mu.p)/sigma.p^2 - 1 )
s <- r * (1-mu.p) / mu.p

cat(sprintf("\n =>  r = %.3f,  s = %.3f \n", r, s))

# check resulting expected value and standard deviation  
mu.r <- r /(r+s)
sigma.r <- sqrt( r*s / ((r+s+1)*(r+s)^2) ) 

cat(sprintf("\n Check:  =>  mu = %.3f,  sigma = %.3f \n\n",
            mu.r, sigma.r ) )

# generate p according to such Beta
p <- rbeta(nevts, r, s) 
hist(p, nc=100, prob=TRUE, col='cyan')
cat(sprintf('\n p:   Mean = %.3f,  std = %.3f\n',
            mean(p), sd(p) ))
print(summary(p))     

pausa()

# Generate x (-> 'x1') using the uncertain p 
x1 <- rbinom(nevts, n, p)
#  
hist(x, nc=40, prob=TRUE, col='cyan', xlim=c(7000,10500))
pausa()
hist(x1, nc=100, prob=TRUE, col='gray', add=TRUE)
cat(sprintf('\n x1:   Mean = %.0f,  std = %.0f\n',
            mean(x1), sd(x1) ))
print(summary(x1))    

pausa()

#-----------------------------------------------------------
# Variant with uncertain p, described by a Beta  
# having mean 0.95 e sd 0.025
# -> Note the substantial difference in the Beta!
mu.p    = 0.95
sigma.p = 0.025

# Calculate the Beta parameters r nd s 
# (in the literature also indicated by 'alpha' e 'beta'))
r <- mu.p * ( mu.p*(1-mu.p)/sigma.p^2 - 1 )
s <- r * (1-mu.p) / mu.p

cat(sprintf("\n =>  r = %.3f,  s = %.3f \n", r, s))

# check resulting expected value and standard deviation  
mu.r <- r /(r+s)
sigma.r <- sqrt( r*s / ((r+s+1)*(r+s)^2) ) 

cat(sprintf("\n Check:  =>  mu = %.3f,  sigma = %.3f \n\n",
            mu.r, sigma.r ) )

# generate p according to such Beta
p <- rbeta(nevts, r, s) 
hist(p, nc=100, prob=TRUE, col='cyan')
cat(sprintf('\n p:   Mean = %.3f,  std = %.3f\n',
            mean(p), sd(p) ))
print(summary(p))     

pausa()

# Generate x (-> 'x2') using the uncertain p 
x2 <- rbinom(nevts, n, p)
#  
hist(x, nc=40, prob=TRUE, col='cyan', xlim=c(7000,10500))
pausa()
hist(x2, nc=100, prob=TRUE, col='gray', add=TRUE)
cat(sprintf('\n x1:   Mean = %.0f,  std = %.0f\n',
            mean(x2), sd(x2) ))
print(summary(x2))