COLS=c("black","black","black") 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, 10, 20, 40, 65, 80) ) { 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 = exp(seq(-0.5,0.5,0.001)) * SIGMA ; 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.05,0), col="white",log="x") ; #points(sigma.range,LL.annuli, type="l", col=COL[2], lwd=9) ; #points(sigma.range,LL.scores, type="l", col=COL[3], lwd=3) max.ll.raw.data = - n * log(mle.annuli^2) - n LL.raw.data = (- n * log(sigma.range^2) - n * (mle.annuli/sigma.range)^2) - max.ll.raw.data points(sigma.range[LL.raw.data>criterion], LL.raw.data[LL.raw.data>criterion], type="l", col=COLS[1], lwd=2) #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) ; 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=COL[3], lwd=3) f.annuli=(pmf(SIGMA))[1,1:7]; #points(mle.annuli+(1:7),f.annuli+criterion, type="h", col=COL[2], lwd=1) Z = 0.2 s.min = mle.annuli / exp( Z * sqrt(1/I.scores) ) ; s.min s.max = mle.annuli * exp( Z * sqrt(1/I.scores) ) ; s.max edge = 1.75 F = 3.25 a = ( (86:89)+0.99 ) /(360/(2*pi)) y.y = sin( a ) ; y.y x.x = cos( a ) ; x.x a = c(1,1, 1,1, 46,46) LTY=c(rep(3,5),1) COL="grey40"; LWD=c(rep(0.75,5), 1) for (i in 1:6) { A = ( (a[i]:89)+0.99 ) /(360/(2*pi)) y.y = sin( A ) ; y.y x.x = cos( A ) ; x.x lines( mle.annuli*exp( x.x*3.25*sqrt(1/I.scores)*(R[i]/R[6] ) ), criterion + c(edge* (R[i]/R[6] ) * y.y) , lwd=LWD[i], col=COL, lty = LTY[i] ) lines( mle.annuli*exp(-x.x*3.25*sqrt(1/I.scores)*(R[i]/R[6] ) ), criterion + c(edge* (R[i]/R[6] ) * y.y) , lwd=LWD[i], col=COL , lty = LTY[i] ) lines( mle.annuli*exp(-x.x*3.25*sqrt(1/I.scores)*(R[i]/R[6] ) ), criterion - c(edge* (R[i]/R[6] ) * y.y) , lwd=LWD[i], col=COL , lty = LTY[i] ) lines( mle.annuli*exp( x.x*3.25*sqrt(1/I.scores)*(R[i]/R[6] ) ), criterion - c(edge* (R[i]/R[6] ) * y.y) , lwd=LWD[i], col=COL, lty = LTY[i] ) if(i>1) text( mle.annuli, criterion+ edge*mid.R[i]/R[6], toString( round(1000*f.annuli[i],0)) , cex=0.75, adj=c(0.5,0.5), col=COL, font=2) if(i==1) text( mle.annuli, criterion, toString( round(1000*f.annuli[i],0)) , cex=0.75, adj=c(0.5,0.5), col=COL, font=2) } text( mle.annuli, criterion+ edge, toString( round(1000*f.annuli[7],0)) , cex=0.75, adj= c(0.5,-0.5), col=COL , font=2) points(sigma.range,LL.annuli, type="l", col=COLS[2], lwd=2, lty=2) ; points(sigma.range,LL.scores, type="l", col=COLS[3], lwd=2, lty=3) text(exp(1.7/sqrt(I.scores))*mle.annuli,criterion-0.1,expression(sigma), cex=1.25 , adj=c(0.5,-0.1)) text(sigma.range[1],-0.01,expression(paste(plain(log),"Lik(", sigma, ")")) , cex=1.25, font=4 , adj=c(0.1,0.5)) text(sigma.range[length(sigma.range)], 0, "Information", adj=c(1,0.5), cex=1.25, font=2) yy=-0.11 if(SIGMA==5 ) segments( sigma.range[1], yy, sigma.range[50], -0.11, col=COLS[1], lwd=2) text(sigma.range[length(sigma.range)], yy, "locations: 200", adj=c(1,0.5), col=COLS[1] , cex=1.0) yy=-0.22 if(SIGMA==5 ) segments( sigma.range[1], yy, sigma.range[50], yy, col=COLS[2], lwd=2, lty=2) text(sigma.range[length(sigma.range)], yy, paste("rings:",toString( round(I.annuli,0)) ), adj=c(1,0.5), col=COLS[2] , cex=1.0) yy=-0.33 if(SIGMA==5 ) segments( sigma.range[1], yy, sigma.range[50], yy, col=COLS[3], lwd=2, lty=3) text(sigma.range[length(sigma.range)], yy, paste("scores:",toString( round(I.scores,0)) ), adj=c(1,0.5), col=COLS[3] , cex=1.0) efficiency = paste( toString( round(100*I.annuli / I.scores,0) ), "%" , sep="") ; } # 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* ######## unequal variances x=seq(0.01, 5.99, 0.01) f1 = dchisq(x,1) f2 = 0.5*dchisq(x,1) y = convolve(f1, f2, type="open") ; n=length(y) plot(1:n,y)