R = c(6.35,15.9,99,107,162,170)
mid.R = ( R + c(0,R[1:5])) / 2


numbers=c(13,4,18,1,20,5,12,9,14,11,8,16,7,19,3,17,2,15,10,6)

annulus= function(xy.matrix) {
	L = length(xy.matrix)/2
	result=c()
	for(i in 1:L){
		ring.no = 7;
		r = sqrt(xy.matrix[i,1]^2 + xy.matrix[i,2]^2)
		if(r <= R[3] && r > R[2]) ring.no=3   
	    if(r <= R[6] && r > R[5]) ring.no=6
	    if(r <= R[5] && r > R[4]) ring.no=5
	    if(r <= R[4] && r > R[3]) ring.no=4
	    if(r <= R[2] && r > R[1]) ring.no=2
	    if(r <  R[1])             ring.no=1
		result=c(result, ring.no)
		}
	return(result)
	}

#  xy=matrix(1:400,nrow=200) ; xy[,2] = 0; xy ; annulus(xy)

score = function(xy.matrix) {
	    ring.no = annulus(xy.matrix) ; ring.no
		r = sqrt(xy.matrix[,1]^2 + xy.matrix[,2]^2) ; r
		a = (360/(2*pi))*acos(abs(xy.matrix[,1]) / r ); a
		a =          a*(xy.matrix[,1] > 0 & xy.matrix[,2] > 0)   +
		       (180-a)*(xy.matrix[,1] < 0 & xy.matrix[,2] > 0)   + 
		       (180+a)*(xy.matrix[,1] < 0 & xy.matrix[,2] < 0)   + 
		       (360-a)*(xy.matrix[,1] > 0 & xy.matrix[,2] < 0) ; a
	    angle = (a > 7.5 & a <= 360) * (a-7.5) + (a <= 7.5) * (352.5+a) ; angle
	    slice.no = ceiling(angle/20) ; slice.no
	    score = 50*(ring.no==1) + 25*(ring.no==2) +
	            1*(ring.no==3)*numbers[slice.no]  +
	            3*(ring.no==4)*numbers[slice.no]  +
	            1*(ring.no==5)*numbers[slice.no]  +
	            2*(ring.no==6)*numbers[slice.no]  +
	            0*(ring.no==7)  
		return(score)
}
	   

pmf = function(sigma) {
	pmf=matrix(rep(10^(-20),122),nrow=2)
	CDF = pchisq( (c(R, 1000)/sigma)^2 , 2)
	pmf.annuli = ( CDF - c(0,CDF[1:6]) ) + 10^(-20)
	pmf[1,1:7]      = pmf.annuli # off the board
	
	pmf[2,c(50,25)] =                    pmf.annuli[1:2]
	pmf[2,1:20]     = pmf[2,1:20]     + (pmf.annuli[3]+pmf.annuli[5])/20
	pmf[2,3*(1:20)] = pmf[2,3*(1:20)] +  pmf.annuli[4]/20
	pmf[2,2*(1:20)] = pmf[2,2*(1:20)] +  pmf.annuli[6]/20
	pmf[2,61]       = pmf[1,7]; # off the board
	return( pmf )
	}

log.lik = function(sigma, freq.annuli, freq.scores ) {
	
		pmf.annuli = (pmf(sigma))[1,1:7];
		log.lik.annuli = sum(freq.annuli*log(pmf.annuli) )
		
		pmf.scores = (pmf(sigma))[2,];
		log.lik.scores = sum( freq.scores*log(pmf.scores) )
		
		return( c(log.lik.annuli, log.lik.scores) )
}


### E-M algorithm

mle.EM = function( freq.annuli, starting.sigma ) {
 		
 	  n.throws = sum(freq.annuli) ;
 	  
      # work with (unobserved) r.v  = squared radius r2, expected value E.r2
  
 	  lower.r2 = (c(0,R))^2 ; upper.r2 = (c(R,999))^2 
 	  
 	  sigma.hat = starting.sigma 
 	  
 	  n.iterations = 0 ; change.in.sigma.hat = 9999
 	  
 	  while ( n.iterations < 100 & change.in.sigma.hat > 0.00001 ) {
 	
 	      E.r2 = 2 * sigma.hat^2
 	    
          mu.given.annulus = lower.r2 + E.r2 + (lower.r2-upper.r2) * exp(lower.r2/E.r2) / 
                                   (exp(upper.r2/E.r2) - exp(lower.r2/E.r2)) 	     
 	      new.sigma.hat = sqrt( 0.5 * sum(mu.given.annulus*freq.annuli)/n.throws )
 	     
 	      change.in.sigma.hat = abs(new.sigma.hat - sigma.hat)
 	     
 	      n.iterations = n.iterations + 1
 	      sigma.hat = new.sigma.hat
 	   }
 	   return ( c(sigma.hat, n.iterations) )
}
 	  

#####

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

n = 50

criterion = -2;

for (SIGMA in c(5, 20, 40, 60, 80, 100) ) {
	 


pmf.annuli          = (pmf(SIGMA))[1,1:7];

E.freq.annuli.SIGMA = n * pmf.annuli
E.freq.annuli.SIGMA[E.freq.annuli.SIGMA < 10^(-15)] = 0
E.freq.annuli.SIGMA
sum(E.freq.annuli.SIGMA)

pmf.scores          = (pmf(SIGMA))[2,];
E.freq.scores.SIGMA = n * pmf.scores
E.freq.scores.SIGMA[E.freq.scores.SIGMA < 10^(-10)] = 0
E.freq.scores.SIGMA
sum(E.freq.scores.SIGMA)


sigma = seq(1,150,0.1) ; # sigma
ll.annuli= c(); ll.scores=c()
for ( i in 1:length(sigma) ) {
	ll = log.lik(sigma[i], E.freq.annuli.SIGMA, E.freq.scores.SIGMA  )
	ll.annuli=c(ll.annuli,ll[1])
	ll.scores=c(ll.scores, ll[2])
}

max.ll.annuli = max(ll.annuli) ; ll.0.annuli = ll.annuli - max.ll.annuli
max.ll.scores = max(ll.scores) ; ll.0.scores = ll.scores - max.ll.scores
c(max.ll.annuli,max.ll.scores)

mle.annuli = sigma[ll.0.annuli==0] ; mle.annuli
mle.scores = sigma[ll.0.scores==0] ; mle.scores

criterion = -2;
sigma.range = sigma[ll.0.scores > criterion | ll.0.annuli > criterion ];
LL.annuli = ll.0.annuli[ll.0.scores > criterion | ll.0.annuli > criterion]
LL.annuli[ LL.annuli <=  criterion ] = NA
LL.scores = ll.0.scores[ll.0.scores > criterion | ll.0.annuli > criterion]

plot(sigma.range,LL.annuli, type="l", ylim=c(criterion,0), col="blue",log="x", lwd=3) ; 
points(sigma.range,LL.scores, type="l", col="red", lwd=3)

#text(log(SIGMA), log.lik(SIGMA, E.freq.annuli.SIGMA, E.freq.scores.SIGMA ), "x" )

delta=0.001

d.0.annuli=
(
(log.lik( exp(log(mle.annuli)+delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1] - 
(log.lik( exp(log(mle.annuli)-delta/2),  E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1] 
)/delta
d.0.annuli

d.L.annuli=
(
(log.lik( exp(log(mle.annuli)-1*delta/2), E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1] - 
(log.lik( exp(log(mle.annuli)-3*delta/2 ) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1]
)/delta
d.L.annuli

d.R.annuli=
(
(log.lik( exp(log(mle.annuli)+3*delta/2), E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1] - 
(log.lik( exp(log(mle.annuli)+1*delta/2 ) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[1]
)/delta
d.R.annuli


I.annuli = -(d.R.annuli - d.L.annuli) / (2*delta) ; I.annuli ; 

d.0.scores=
(
(log.lik( exp(log(mle.scores)+delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2] - 
(log.lik( exp(log(mle.scores)-delta/2) ,          E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2]
)/delta
d.0.scores

d.L.scores=
(
(log.lik( exp(log(mle.scores)-1*delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2] - 
(log.lik( exp(log(mle.scores)-3*delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2]
)/delta
d.L.scores

d.R.scores=
(
(log.lik( exp(log(mle.scores)+3*delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2] - 
(log.lik( exp(log(mle.scores)+1*delta/2) , E.freq.annuli.SIGMA, E.freq.scores.SIGMA) )[2]
)/delta
d.R.scores

I.scores = -(d.R.scores - d.L.scores) / (2*delta) ; I.scores ; 

c(I.annuli, I.scores) ; 

text(mle.annuli, 0, toString( round(I.annuli,1) ), adj=c(-1,0), col="blue" , cex=1.25) 

text(mle.scores, 0, toString( round(I.scores,1) ), adj=c(2,0), col="red" , cex=1.25) 

efficiency =  paste( toString( round(100*I.annuli / I.scores,0) ), "%" , sep="") ;

text(sigma.range[length(sigma.range)], -0.1, efficiency, cex=2, adj=c(1,1)) 


ci.annuli=round(exp( log(mle.annuli) + c(-1.96,0,1.96)*sqrt(1/I.annuli) ),1)
ci.annuli
ci.scores=round(exp( log(mle.scores) + c(-1.96,0,1.96)*sqrt(1/I.scores) ),1)
ci.scores

log(ci.annuli[2])-log(ci.annuli[1])

log(ci.scores[2])-log(ci.scores[1])

print( mle.EM(E.freq.annuli.SIGMA, 50) )

f.scores=pmf(SIGMA)[2,]; 
# plot((1:61)-0.18,f.scores, type="h", ylim=c(0,1), col="red", lwd=3)
f.annuli=(pmf(SIGMA))[1,1:7]; 
#points(mle.annuli+(1:7),f.annuli+criterion, type="h", col="blue", lwd=1)

fx=0.975 ; edge = 1.75

segments( fx*mle.annuli, criterion + edge, (1/fx)*mle.annuli, criterion + edge )
segments( fx*mle.annuli, criterion , (1/fx)*mle.annuli, criterion  )
segments( fx*mle.annuli, criterion , (fx)*mle.annuli, criterion + edge )
segments( (1/fx)*mle.annuli, criterion , (1/fx)*mle.annuli, criterion + edge )

for (i in 1:6) {
 if(i<6) segments( fx*mle.annuli, criterion+ edge*R[i]/R[6], (1/fx)*mle.annuli, criterion+ edge*R[i]/R[6] )
 text( mle.annuli, criterion+ edge*mid.R[i]/R[6],  
 toString( round(1000*f.annuli[i],0)) , cex=1)
}
text( mle.annuli, criterion+ edge,  
 toString( round(1000*f.annuli[7],0)) , cex=1, adj= c(0.5,-0.5) )
 
} 

# end loop

# simulated throws

SIGMA = 80

simulated.throws  = function(n,sigma)  matrix(rnorm(2*n,0,sigma),nrow=n)
	
n=200 ; data.n =simulated.throws(n,SIGMA) ; data.n

annuli = annulus( data.n ) ; annuli; summary(annuli )

S = score(data.n) ;  noquote(paste(S,collapse=","))

SIGMA = seq(5,100,0.1) ; 
ll = log.lik(SIGMA, annuli) ;
plot(SIGMA,ll, type="l")

ML = max(ll)
MLE = SIGMA[ML == max(ll)








mu = 0.3 ; L=0.4 ; U=1.6

L+mu+(L-U)*exp(L/mu)/ (exp(U/mu) - exp(L/mu))

Y=rexp(100000,1/mu) ; #mean(Y)

Y.1.2 = Y[ Y >=L & Y < U]; mean(Y.1.2)

p.total = exp(-L/mu) -  exp(-U/mu) ; p.total

p.rectangle = (U-L)* (1/mu)*exp(-U/mu) ; p.rectangle

p.triangle = p.total - p.rectangle ; p.triangle

w.triangle  = p.triangle / (p.triangle + p.rectangle)
w.rectangle = p.rectangle / (p.triangle + p.rectangle)

c(w.triangle,w.rectangle)

mean.L.U = w.rectangle*(L+U)/2 + w.triangle*