#---------------------------------------------
# inferring p, predicting n1 -> use hit/miss
#
# GdA, May 2020
#---------------------------------------------

#model parameters
n0 = 20
x0 = 10
n1 = 10

# MC parameters
N = 20000

# unnormalized distribution 
uf <- function(p, x1, n0, x0, n1) {
   if(p<0 || p>1)     return(0)
   if(x1 <0 || x1>n1) return(0)
   return( p^x0*(1-p)^(n0-x0) * p^x1*(1-p)^(n1-x1) /
           (factorial(x1)*factorial(n1-x1)) ) 
}

# find the maximum
p  <- seq(0,1,len=101)
x1 <- 0:n1
f.max = 0
for (i in 1:length(p)) {
    for (j in 1:length(x1)) {
      f =  uf(p[i], x1[j], n0, x0, n1)
    if( f > f.max ) f.max = f        
  } 
}
f.max <- f.max * 1.1   # exaggerated safety factor

# sample in 2D
#1. random choice in the plane (p,x1) 
p.r  <- runif(N)
x1.r <- sample(0:n1, N, rep=TRUE)
# 2. Calculate the function in correspondence of the extracted points
f <- numeric(N)
for (i in 1:N) {
  f[i] <- uf(p.r[i], x1.r[i], n0, x0, n1)   
}
# 3. accept events
hit <-  runif(N)*f.max <= f  # accepted events ('hit')
p   <- p.r[hit] 
x1  <- x1.r[hit]
# efficiency: sum(hit)/N

# analysis and plots of the chain
par(mfcol=c(2,2))
hist(p, nc=50, col='cyan', xlim=c(0,1))
barplot(table(x1)/N, col='cyan', xlab='x1')    #,xlim=c(0,n1))
plot(x1, p, col='blue', cex=1)
plot(NULL,  xaxt = 'n', yaxt = 'n', bty = 'n', pch = '',
     ylab = '', xlab = '', xlim=c(0,10), ylim=c(0,10))
text(4,9, sprintf("mean(p) = %.3f\n sigma(p) =  %.3f", mean(p), sd(p)), cex=1.5,  col='blue')
text(4,5, sprintf("mean(x1) = %.3f\n sigma(x1) =  %.3f", mean(x1), sd(x1)), cex=1.5, col='blue')
text(4,1, sprintf("rho(p,x1) = %.3f", cor(p,x1)), cex=1.5, col='blue')
par(mfcol=c(1,1))