1 / 73

Assignment 1

Assignment 1. Team members: Yubin Hou , Xuan Li, Nanwei Wang. Key steps. 1.Constructing functions 2.Graphing 3.Using different methods to find the root 4.Discussion. Functions. Graph. cauchy.log(theta,1) x <- seq (-10,10,0.01) y <- numeric(2001) for (i in 1:2001){

thuy
Download Presentation

Assignment 1

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Assignment 1 Team members: YubinHou, Xuan Li, Nanwei Wang

  2. Key steps 1.Constructing functions 2.Graphing 3.Using different methods to find the root 4.Discussion

  3. Functions

  4. Graph cauchy.log(theta,1) x <- seq(-10,10,0.01) y <- numeric(2001) for (i in 1:2001){ y[i]<- -n*log(pi)-sum(log(1+(z-x[i])^2)) } plot(x,y,type="l")

  5. Graphs

  6. Graphs #cauchy.log.derivative x <- seq(-10,10,0.01) y <- numeric(2001) for (i in 1:2001){ y[i]<- clogp(x[i]) } plot(x,y,type="l") par(new="T") y0 <- numeric(2001) plot(x,y0,type="l",col="red")

  7. Graphs

  8. Graphs #cauchy.log.derivativeII x <- seq(-10,10,0.01) y <- numeric(2001) for (i in 1:2001){ y[i]<- clogpp(x[i]) } par(new="F") plot(x,y,type="l")

  9. Graphs

  10. Newton Method diff <- 4 iter <- 0 x <- mean(z) while((diff>0.001) && (iter<30)){ oldx <- x x <- x-clogp(x)/clogpp(x) diff <- abs(x-oldx) iter <- iter+1 print(c(iter,diff)) } x;clogp(x)

  11. Newton Method x1 <- c(-11,-1,0,1.5,4,4.7,7,8,38) oldx1 <- numeric(length(x1)) diff <- 4 iter <- 0 for(i in 1:9){ print(x1[i]) diff <- 4 iter <- 0 while((diff>0.001) && (iter<30)){ oldx1[i] <- x1[i] x1[i] x1[i] <- x1[i]-clogp(x1[i])/clogpp(x1[i]) k <- -clogp(x1[i])/clogpp(x1[i]) diff <- abs(x1[i]-oldx1[i]) iter <- iter+1 print(c(iter,diff,k)) } print(x1[i]) print(clog(x1[i])) print('end') }

  12. Newton Method Start with mean= 5.106 [1] 1.00000 12.19122 [1] 2.00000 14.56064 [1] 3.00000 16.11695 [1] 4.000000 7.899357 [1] 5.000000 1.482321 [1] 6.000000 2.533977 [1] 7.00000000 0.05435966 [1] 8.0000000000 0.0002499021 Root= 54.87662

  13. Newton Method Start with -11 [1] 1.000000 5.493309 [1] 2.000000 8.797295 …. [1] 29 1681784824 [1] 30 3363569649 cannot find the root Start with -1 [1] 1.0000000 0.7585504 [1] 2.0000000 0.0499397 [1] 3.0000000000 0.0007765955 Root = -0.1922865

  14. Newton Method Start with 0 [1] 1.0000000 0.1963366 [1] 2.000000000 0.004054071 [1] 3.000000e+00 4.088535e-06 Root = -0.1922865 Start with 1.5 [1] 1.0000000 0.1756107 [1] 2.00000000 0.03670085 [1] 3.000000000 0.001273843 [1] 4.000000e+00 1.394325e-06 Root = 1.713587

  15. Newton Method Start with 4 [1] 1.000000 1.061099 [1] 2.000000 0.104381 [1] 3.00000000 0.01662494 [1] 4.0000000000 0.0004225948 Root= 2.817472 Start with 4.7 [1] 1.000000 5.446616 [1] 2.0000000 0.6614478 [1] 3.0000000 0.1064711 [1] 4.0000000000 0.0006475279 Root= -0.1922865

  16. Newton Method Start with 7 [1] 1.000000 4.835393 [1] 2.000000 9.700744 [1] 3.00000 17.41333 [1] 4.000000 3.048335 [1] 5.000000 1.358983 [1] 6.0000000 0.4326834 [1] 7.00000000 0.03051532 [1] 8.0000000000 0.0001414995 Root= 41.04085

  17. Newton Method Start with 8 [1] 1.00000 5.80202 [1] 2.00000 11.56717 … [1] 29 1172669668 [1] 30 2345339335 Cannot find the root Start with 38 [1] 1.000000 4.957268 [1] 2.0000000 0.1401467 [1] 3.00000000 0.02118844 [1] 4.0000000000 0.0005549434 Root = 42.79538

  18. Method comparisons

  19. Bisection method ∙code: diff <- 4 > iter <- 0 > a <- -1 > b <- 1 > x <- (a + b)/2 > while ((diff > 1e-04) && (iter < 30)) { + oldx <- x + olda <- a + oldb <- b + if (fq(a) * fq(x) <= 0) { + a <- a + b <- x + x <- (a + b)/2 + } else { + a <- x + b <- b + x <- (a + b)/2 + } + diff <- abs(x - oldx) + iter <- iter + 1 + print(c(iter, diff))

  20. Bisection method [a,b]=[-1,1] [1] 1.0 0.5 [1] 2.00 0.25 [1] 3.000 0.125 [1] 4.0000 0.0625 [1] 5.00000 0.03125 [1] 6.000000 0.015625 [1] 7.0000000 0.0078125 [1] 8.00000000 0.00390625 [1] 9.000000000 0.001953125 [1] 1.000000e+01 9.765625e-04 [1] 1.100000e+01 4.882812e-04 [1] 1.200000e+01 2.441406e-04 [1] 1.300000e+01 1.220703e-04 [1] 1.400000e+01 6.103516e-05

  21. Bisection >x [1] -0.1923218 > f(x) [1] -72.91582 > fq(x) [1] 0.0001083593

  22. Bisection method [a,b]=[2,3] [1] 1.00 0.25 [1] 2.000 0.125 [1] 3.0000 0.0625 [1] 4.00000 0.03125 [1] 5.000000 0.015625 [1] 6.0000000 0.0078125 [1] 7.00000000 0.00390625 [1] 8.000000000 0.001953125 [1] 9.0000000000 0.0009765625 [1] 1.000000e+01 4.882812e-04 [1] 1.100000e+01 2.441406e-04 [1] 1.200000e+01 1.220703e-04 [1] 1.300000e+01 6.103516e-05

  23. Bisection method [a,b]=[2,3] > x [1] 2.817444 > f(x) [1] -74.36046 > fq(x) [1] 6.217143e-05

  24. Bisection method [a,b]=[1,2] [1] 1.00 0.25 [1] 2.000 0.125 [1] 3.0000 0.0625 [1] 4.00000 0.03125 [1] 5.000000 0.015625 [1] 6.0000000 0.0078125 [1] 7.00000000 0.00390625 [1] 8.000000000 0.001953125 [1] 9.0000000000 0.0009765625 [1] 1.000000e+01 4.882812e-04 [1] 1.100000e+01 2.441406e-04 [1] 1.200000e+01 1.220703e-04 [1] 1.300000e+01 6.103516e-05

  25. Bisection method [a,b]=[1,2] > x [1] 1.713562 > f(x) [1] -74.64202 > fq(x) [1] -2.204954e-05

  26. Fixed-point method Code: diff <- 4 x <- -1 iter <- 0 r <- 0.25 while ((diff > 0.001) && (iter < 30)) { oldx <- x x <- r * fq(x) + x diff <- abs(x - oldx) iter <- iter + 1 print(c(iter, diff)) } x fq(x) f(x)

  27. Fixed-point method starting at -1 R=0.25 [1] 1.0000000 0.4803539 [1] 2.0000000 0.2196615 [1] 3.00000000 0.08015055 [1] 4.00000000 0.02106847 [1] 5.000000000 0.004983352 [1] 6.000000000 0.001151817 [1] 7.0000000000 0.0002648679

  28. Fixed point method starting at -1 R=0.25 > x [1] -0.1923656 > fq(x) [1] 0.0002433502 > f(x) [1] -72.91582

  29. R=0.64 [1] 1.000000 1.229706 [1] 2.0000000 0.7803109 [1] 3.0000000 0.6061322 … … [1] 28.0000000 0.1247265 [1] 29.0000000 0.1205457 [1] 30.0000000 0.1163133

  30. Fixed point method starting at -1 R=0.64 > x [1] -0.2502889 > fq(x) [1] 0.1757899 > f(x) [1] -72.92095

  31. R=1 [1] 1.000000 1.921416 [1] 2.000000 1.247424 [1] 3.0000000 0.3938901 … … [1] 29.00000 1.35571 [1] 30.000000 1.585071

  32. Fixed point method starting at 1 R=1 > x [1] -0.9775387 > fq(x) [1] 1.865197 > f(x)

  33. Secant method diff<-4 iter<-0 x1<- -10 x2<- 10 while ((diff>0.0001) && (iter<30)) { oldx1<-x1 oldx2<-x2 x2<-x2-fq(x2)*(x2-x1)/(fq(x2)-fq(x1)) x1<-oldx2 diff<-abs(x2-oldx2) iter<-iter+1 print(c(iter,diff)) } x2 fq(x2) f(x2)

  34. Secant method (-2, -1) [1] 1.0000000 0.5587942 [1] 2.0000000 0.3166416 [1] 3.00000000 0.07383809 [1] 4.000000000 0.006036734 [1] 5.00000e+00 7.90164e-05 > x2 [1] -0.1922865 > fq(x2) [1] -3.682938e-07 > f(x2) [1] -72.91582

  35. Secant method (-3, 3) [1] 1.0000000 0.3142115 [1] 2.000000 0.097743 [1] 3.0000000 0.0425518 [1] 4.000000000 0.009069993 [1] 5.0000000000 0.0004528824 [1] 6.000000e+00 5.943044e-06 > x2 [1] 2.817472 > fq(x2) [1] -9.092855e-09 > f(x2) [1] -74.36046

  36. Secant method (-10, 10) [1] 1.00000 11.69632 [1] 2.000000 5.975252 [1] 3.000000 3.941006 [1] 4.0000000 0.8404016 [1] 5.0000000 0.3122807 [1] 6.000000000 0.002378905 [1] 7.0000000000 0.0002878937 [1] 8.000000e+00 1.452314e-07 > x2 [1] -0.1922866 > fq(x2) [1] -3.130409e-11 > f(x2) [1] -72.91582

  37. Method comparisons

  38. Implementation in normal distribution

  39. ●get the random data with standard normal distribution: the code: >data1 <- rnorm(20, mean = 0, sd = 1) >data1 >mean(data1)

  40. Plot the loglikelihood function: Code: f <- function(x) { q <- (-10) * log(2 * pi) - sum((data1 - x)^2/2) return(q) } fq<- function(x) { q <- sum(data1 - x) return(q) } fqq<- function(x) { q <- -20 return(q) } x <- (-500:500)/100 y <- numeric(1001) for (i in 1:1001) { y[i] <- f(x[i]) } plot(x, y, type = "l")

  41. The loglikelihood function

  42. Graphs Plot the first derivative loglikelihood: Code: yq <- numeric(1001) for (i in 1:1001) { yq[i] <- fq(x[i]) } plot(x, yq, type = "l")

  43. The first derivative loglikelihood

  44. Graphs The second derivative loglikelihood: code: yqq <- numeric(1001) for (i in 1:1001) { yqq[i] <- fqq(x[i]) } plot(x, yqq, type = "l")

  45. The second derivative loglikelihood

  46. Newton method Code: diff <- 4 iter <- 0 x1 <- mean(data1) while ((diff > 0.001) && (iter < 30)) { oldx <- x1 x1 <- x1 - fq(x1)/fqq(x1) diff <- abs(x1 - oldx) iter <- iter + 1 print(c(iter, diff)) }

  47. Newton method ●Start with mean(data1): [1] 1.000000e+00 2.775558e-17 > x [1] -0.06505894 > fq(x) [1] 0.07091004

  48. Newton method ●start with -11: [1] 1.00000 10.77859 [1] 2.000000e+00 2.775558e-17 > x [1] -0.06505894 > fq(x) [1] 0.07091004

  49. Newton method ●start with -1: [1] 1.0000000 0.7785914 [1] 2.000000e+00 2.775558e-17 > x [1] -0.06505894 > fq(x) [1] 0.07091004

  50. Newton method start with 0: [1] 1.0000000 0.2214086 [1] 2 0 >x [1] -0.06505894 > fq(x) [1] 0.07091004

More Related