##  2 : 136 split when denomiantors are 6:94


# large sample binomial

se = sqrt((2/138)*(136/138)/138) ; se
[1] 0.01017339


> P = 2/138 + 1.645*c(-1,0,1)*se


> 94*P/(6*(1-P))
[1] -0.03505355  0.23039216  0.50500884

> PU = 0.04 ; pbinom(2,138,PU)
[1] 0.0828316
> PU = 0.035 ; pbinom(2,138,PU)
[1] 0.1350635
> PU = 0.045 ; pbinom(2,138,PU)
[1] 0.04956551   # close enough


> PL = 0.005 ; 1-pbinom(1,138,PL)
[1] 0.1520662
> PL = 0.0075 ; 1-pbinom(1,138,PL)
[1] 0.2771612
> PL = 0.002 ; 1-pbinom(1,138,PL)
[1] 0.03160267
> PL = 0.003 ; 1-pbinom(1,138,PL)
[1] 0.06510314
> PL = 0.0025 ; 1-pbinom(1,138,PL)
[1] 0.04724291  # again close enough


> P=c(PL,PU);P
[1] 0.0025 0.0450

> 94*P/(6*(1-P))
[1] 0.03926483 0.73821990


# log.lik for theta

log.lik = function(theta) {
	P = 6*theta/(6*theta+94)
	return( 2 * log(P) + 136 * log(1-P) ) 
} 

theta.values = 2^seq(-7,10, 0.1) ; theta.values

max.log.lik = max(log.lik(theta.values))

theta.ML = theta.values[which(log.lik(theta.values)==max.log.lik)]
theta.ML

Lik.interval = theta.values[max.log.lik - log.lik(theta.values) <  1.353]

round( c(Lik.interval[1], Lik.interval[length(Lik.interval)]), 2)


par(mfrow=c(3,1),mar = c(4,4,3,1))

XLIM=c(0.001,5)

plot(theta.values,exp(log.lik(theta.values)), 
  xlim=XLIM,
  ylim=exp(c(max.log.lik- 3, max.log.lik)),
  ylab="lik[theta]")


# Posterior for theta if prior for theta is...

## theta ~ gamma(shape=100, scale=0.01)

SHAPE=2; SCALE = 1/2
prior = dgamma(theta.values,shape=SHAPE, scale=SCALE)
plot(theta.values, prior, xlim=XLIM)

prior.times.lik = prior * exp(log.lik(theta.values))

plot(theta.values, prior.times.lik,
  ylab="prior x Likelihood", xlim=XLIM,
  main="POSTERIOR, but needs rescaling")

# now integrate to get constant of integration, and rescale


########################################################

##### other better transforms of the theta parameter...

par(mfrow=c(2,2),mar = c(4,4,3,1))

XLIM = c(0.001, 1.5)

plot(theta.values,log.lik(theta.values), 
  xlim=XLIM,
  ylim=c(max.log.lik- 3, max.log.lik),
  ylab="log.lik[theta]")
abline(h=max.log.lik-1.353)

psi=log2(theta.values)
plot(psi,log.lik(theta.values),
  xlab="psi = log2[theta]",ylab="log.lik[psi]", 
  xlim=log2(XLIM) ,
  ylim=c(max.log.lik- 3, max.log.lik),)
abline(h=max.log.lik-1.353)

psi = theta.values^(1/3)   # power of 1/3 often works for (0,Inf) parameters
plot(psi,log.lik(theta.values),
  xlab="psi = theta^(1/3)",
  ylab="log.lik[psi]",
  xlim=XLIM^(1/3),
  ylim=c(max.log.lik- 3, max.log.lik), main="most 'normal-looking' of 3?" )
abline(h=max.log.lik-1.353)

plot(psi,exp(log.lik(theta.values)),
 xlab="psi = theta^(1/3)",
 ylab="lik[psi]",
 xlim=XLIM^(1/3) )