n1 = c(5,6,4)
n2 = c(5,4,6)
y = c(7,5,6)
L =  apply(cbind(0,y-n2),1,max)
U = apply(cbind(y,n1),1,min)

likelihood = function(theta1,theta2){
sum(dbinom(L[1]:U[1],n1[1],theta1)*dbinom(y[1]-L[1]:U[1],n2[1],theta2))*
sum(dbinom(L[2]:U[2],n1[2],theta1)*dbinom(y[2]-L[2]:U[2],n2[2],theta2))*
sum(dbinom(L[3]:U[3],n1[3],theta1)*dbinom(y[3]-L[3]:U[3],n2[3],theta2))
}

N = 50
theta1 = seq(0.01,0.99,length=N)
theta2 = seq(0.01,0.99999,length=N)
gamma1 = log(theta1/(1-theta1))
gamma2 = log(theta2/(1-theta2))
like = matrix(0,N,N)
for (i in 1:N)
  for (j in 1:N)
       like[i,j] = likelihood(theta1[i],theta2[j])
  
par(mfrow=c(1,2))
contour(theta1,theta2,like,nlevels=20,xlab=expression(theta[1]),ylab=expression(theta[2]))
title("Likelihood function")
contour(gamma1,gamma2,like,nlevels=20)

persp(theta1,theta2,like, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed", cex.axis = 0.8)
persp(gamma1,gamma2,like, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed", cex.axis = 0.8)

M         = 50000
thetas1.d = runif(M)
thetas2.d = runif(M)
gammas1.d = log(thetas1.d/(1-thetas1.d))
gammas2.d = log(thetas2.d/(1-thetas2.d))

w = rep(0,M)
for (i in 1:M)
  w[i]  = likelihood(thetas1.d[i],thetas2.d[i])
ind = sample(1:M,size=M,replace=TRUE,prob=w)
thetas1 = thetas1.d[ind]
thetas2 = thetas2.d[ind]
gammas1 = log(thetas1/(1-thetas1))
gammas2 = log(thetas2/(1-thetas2))

par(mfrow=c(1,2))
plot(thetas1.d,thetas2.d)
points(thetas1,thetas2,col=2)
contour(theta1,theta2,like,col=3,add=TRUE,nlevels=30)
plot(gammas1,gammas2)
contour(gamma1,gamma2,like,nlevels=20,add=TRUE,col=3)


par(mfrow=c(2,2))
hist(thetas1,prob=TRUE)
abline(h=1,lwd=2)
hist(thetas2,prob=TRUE)
abline(h=1,lwd=2)
hist(gammas1,prob=TRUE)
lines(density(gammas1.d),lwd=2)
hist(gammas2,prob=TRUE)
lines(density(gammas2.d),lwd=2)