#   Homework 2 R Code    2/10/2011

pr1 = function(n=100,cc=seq(.8,1.4,,151)) {
    y = cc;  cc0 = cc <= 1;  cc1 = cc>1
    y[cc0] = n/((n-1)*cc[cc0]) - 1
    y[cc1] = 1 + (2*cc[cc1]^(-n) - n/cc[cc1])/(n-1)
    plot(cc,y,type="l",ylim=range(0,y),xlab="c",ylab="ISE(c)")
    abline(v=1,lty=2); abline(v=2^(1/(n-1)),col=2)
    text(1.15,.15,paste("n =",n))
}

pr5 = function() {
   logseq = function(a,b,len) { exp(seq(log(a),log(b),len=len)) }
   bestmin = function(z) { i = seq(nrow(z))[apply(z,1,"min")==min(z)]
      j = seq(ncol(z))[z[i,]==min(z)]; return(c(i,j)) }
   f.l2e = function(x,mu,sig) {
      nmu = length(mu);  nsig = length(sig)
      val = matrix(0,nmu,nsig)
      for(i in 1:nmu) { m = mu[i]
          for(j in 1:nsig) { s = sig[j]
	      val[i,j] = 1/(2*sqrt(pi)*s) - 2*mean(dnorm(x,m,s))
      }   }
      plot(range(mu),range(sig),type="n",log="y")
      contour(mu,sig,val,levels=seq(min(val),max(val),,45),add=T,drawlabels=F)
      xs=c(mean(x),sd(x)); points(xs[1],xs[2],pch=16,col=2,cex=2)
      ij=bestmin(val); points(mu[ij[1]],sig[ij[2]],pch=16,col=4,cex=2)
      browser()
      xx=seq(-3,6,,101); yy = dnorm(xx,mu[ij[1]],sig[ij[2]])
      hist(x,prob=T,col=grey(.8),ylim=range(yy)+c(0,.05))
         lines(xx,dnorm(xx,xs[1],xs[2]),col=2,lwd=3)
         lines(xx,yy,col=4,lwd=3)
      browser()
    }
    set.seed(234); x = rnorm(100)
        f.l2e(x,seq(-1,2,.05),logseq(.7,3,51))
    set.seed(245); x = c(rnorm(100), rnorm(25,5))
        f.l2e(x,seq(-1,2,.01),logseq(.7,3,51))
}