# slide 4 (.3-.1)-.2 abs((.3-.1)-.2)<1e-10 identical(.3-.1,.2) isTRUE(all.equal(.3-.1,.2)) .Machine$double.eps # slide 5 ## In math, x*(sqrt(x+1)-sqrt(x)) = x/(sqrt(x+1)+sqrt(x)) because ## ( sqrt(x+1) )^2 - ( sqrt(x) )^2 = 1 x <- 500009999998 x*(sqrt(x+1)-sqrt(x)) x/(sqrt(x+1)+sqrt(x)) # slide 13 f <- function(x) (x-2)*(1+x^2-0.3*(x-1)^3) A <- 1.4 B <- 2.9 for (i in 1:20) { m <- (A+B)/2 if (f(A)*f(m)>0) A <- m else B <- m print(A,digits=12) } f <- function(x) (x-2)*(1+x^2-0.3*(x-1)^3) f1 <- function(x) 1+x^2-0.3*(x-1)^3+(x-2)*(2*x-0.9*(x-1)^2) g <- function(x) x-f(x)/f1(x) xold <- 2.9 for (i in 1:20){ xnew=g(xold) xold=xnew print(xnew,digits=12) } # slide 17 uniroot(f,c(1,3)) # slide 20 f <- function(x) (x-2)*(1+x^2-0.3*(x-1)^3) F <- function(x) x^2/2+x^4/4-0.3*x*(x-1)^4/4-0.3*(x-1)^5/20-2*x-2*x^3/3+0.3*(x-1)^4/2 nlm(F,2.9) optim(2.9,F,method="BFGS") optim(2.9,F,f,method="BFGS") # slide 21 F <- function(x) x^2/2+x^4/4-0.3*x*(x-1)^4/4-0.3*(x-1)^5/20-2*x -2*x^3/3+0.3*(x-1)^4/2 optim(2.9,F,method="Nelder-Mead") # slide 25 x <- rcauchy(100) thetap <- seq(-10,10,length=100) loglik <- function(x,theta) dcauchy(x,theta,log=T) plot(thetap,apply(outer(x,thetap,loglik),2,sum),type="l") thetap <- seq(-50,50,length=100) plot(thetap,apply(outer(x,thetap,loglik),2,sum),type="l") score <- function(x,theta) 2*(x-theta)/(1+(x-theta)^2) plot(thetap,apply(outer(x,thetap,score),2,sum),type="l") abline(h=0) par(mfrow=c(1,2)) x <- rcauchy(100) thetap<- seq(from=min(x)-10,to=max(x)+10,length=100) plot(thetap,apply(outer(x,thetap,loglik),2,sum),type="l") plot(thetap,apply(outer(x,thetap,score),2,sum),type="l") abline(h=0) # slide 26 set.seed(1244) x <- rcauchy(25) q <- seq(min(x),max(x),by=0.01) minloglik <- function(q,x) apply(log(1+outer(x,q,"-")^2),2,sum) minloglik1 <- function(q,x) { v <- outer(x,q,"-") 2*apply(v/(1+v^2),2,sum) } uniroot(minloglik1,median(x)+c(-3,3),x=x) optim(median(x),method="BFGS",minloglik,x=x,hessian=TRUE) optim(median(x),lower=median(x)-3,upper=median(x)+3,method="Brent",minloglik,x=x,hessian=TRUE) # slide 28 optim(-70,method="BFGS",minloglik,x=x) # slide 29 set.seed(1244) x <- rcauchy(25) minloglik <- function(par,x){ q=par[1] s=par[2] length(x)*log(s)+apply(log(1+outer(x,q,"-")^2/s^2),2,sum) } optim(c(median(x),mad(x)),minloglik,x=x) # Default method: Nelder-Mead OP <- optim(c(median(x),mad(x)),method="BFGS",minloglik,x=x) OP # Quasi Newton sqrt(diag(solve(OP$hessian)))