计算机代考程序代写 Practical 10 Solution

Practical 10 Solution

Practical 10 Solution

Libaries and least square estimate needed for Question 1

rm(list=ls())
library(boot)
z<-matrix(0,13,2) z[,1]<-c(0.01,0.48,0.71,0.95,1.19,0.01,0.48,1.44,0.71,1.96,0.01,1.44,1.96) z[,2]<-c(127.6,124.0,110.8,103.9,101.5,130.1,122.0,92.3,113.1,83.7,128.0,91.4,86.2) temp<-lm(z[,2]~z[,1]) temp$coef[2]/temp$coef[1] ## z[, 1] ## -0.1850722 Question 1.1 lm1.bt=function(x,i) { temp<-lm(x[i,2]~x[i,1])$coef ratio<-temp[2]/temp[1] return(ratio) } set.seed(5); boot1<-boot(data=z, stat=lm1.bt, R=20) print(boot1) ## ## ORDINARY NONPARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = z, statistic = lm1.bt, R = 20) ## ## ## Bootstrap Statistics : ## original bias std. error ## t1* -0.1850722 0.000468238 0.006386187 mean(boot1$t)-boot1$t0 ## x[i, 1] ## 0.000468238 sd(boot1$t) ## [1] 0.006386187 boot.ci(boot1, type = "norm") 1 ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 20 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = boot1, type = "norm") ## ## Intervals : ## Level Normal ## 95% (-0.1981, -0.1730 ) ## Calculations and Intervals on Original Scale ba<-boot.array(boot1, indices=T) res<-rep(0,5) for (i in 1:5) { z2<-z[ba[i,],1:2] temp<-lm(z2[,2]~z2[,1]) res[i]<-temp$coef[2]/temp$coef[1] } res-boot1$t ## [,1] ## [1,] 0.000000000 ## [2,] 0.000000000 ## [3,] 0.000000000 ## [4,] 0.000000000 ## [5,] 0.000000000 ## [6,] -0.015447394 ## [7,] -0.001291957 ## [8,] -0.007105341 ## [9,] 0.002143070 ## [10,] -0.013688289 ## [11,] -0.015129832 ## [12,] 0.006220850 ## [13,] -0.013476829 ## [14,] 0.011677981 ## [15,] -0.011416528 ## [16,] -0.009643222 ## [17,] 0.013992953 ## [18,] -0.010335582 ## [19,] 0.006080048 ## [20,] -0.006988010 Question 1.2 lm2.bt=function(x,i) { temp<-lm(x[,2]~x[,1]) y.star<-temp$fitted + temp$residual[i] temp1<-lm(y.star~x[,1])$coef ratio<-temp1[2]/temp1[1] return(ratio) } set.seed(5); boot2<-boot(data=z, statistic=lm2.bt, R=5); boot2 2 ## ## ORDINARY NONPARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = z, statistic = lm2.bt, R = 5) ## ## ## Bootstrap Statistics : ## original bias std. error ## t1* -0.1850722 -0.004011592 0.01066213 print(boot2) ## ## ORDINARY NONPARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = z, statistic = lm2.bt, R = 5) ## ## ## Bootstrap Statistics : ## original bias std. error ## t1* -0.1850722 -0.004011592 0.01066213 boot.ci(boot2, type = "norm") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 5 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = boot2, type = "norm") ## ## Intervals : ## Level Normal ## 95% (-0.2020, -0.1602 ) ## Calculations and Intervals on Original Scale Question 1.3 lm3.bt=function(x) { temp<-lm(x[,2]~x[,1])$coef; ratio<-temp[2]/temp[1]; ratio } lm.gen=function(x, mle) { n<-nrow(x); err<-rnorm(n, mean=0, sd=mle$sigma) y<-mle$beta0+x[,1]*mle$beta1+err return(cbind(x[,1],y)) } 3 temp<-lm(z[,2]~z[,1]) mle.list<-list(beta0=temp$coef[1], beta1=temp$coef[2]) mle.list$sigma<-sqrt(sum(temp$residˆ2)/temp$df.resid) set.seed(5); boot3<-boot(data=z, stat=lm3.bt, R=20, sim="parametric", ran.gen=lm.gen, mle=mle.list) print(boot3) ## ## PARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = z, statistic = lm3.bt, R = 20, sim = "parametric", ## ran.gen = lm.gen, mle = mle.list) ## ## ## Bootstrap Statistics : ## original bias std. error ## t1* -0.1850722 0.0001882346 0.006385126 boot.ci(boot3, type = "norm") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 20 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = boot3, type = "norm") ## ## Intervals : ## Level Normal ## 95% (-0.1978, -0.1727 ) ## Calculations and Intervals on Original Scale Question 1.4 cor.bt=function(x,i) { cor(x[i,1],x[i,2]) } set.seed(5); boot4<-boot(data=z, stat=cor.bt, R=5000); boot4 ## ## ORDINARY NONPARAMETRIC BOOTSTRAP ## ## ## Call: ## boot(data = z, statistic = cor.bt, R = 5000) ## ## ## Bootstrap Statistics : ## original bias std. error ## t1* -0.9847435 -5.809659e-05 0.006875357 4 bc<-boot.ci(boot4, type = "norm") plot(density(boot4$t), ylim=c(0,65), lwd=2) hist(boot4$t, breaks=30, freq=F, add=T) abline(v = boot4$t0, col = "red", lwd = 2) abline(v = bc$normal[,2], col = "green", lwd = 2) abline(v = bc$normal[,3], col = "green", lwd = 2) −1.00 −0.98 −0.96 −0.94 −0.92 −0.90 0 1 0 2 0 3 0 4 0 5 0 6 0 density.default(x = boot4$t) N = 5000 Bandwidth = 0.00104 D e n si ty 5 Libaries and data needed for Question 2 rm(list=ls()) library(e1071) set.seed(114) f1 <-1 f2 <-1 m1 <- cbind(rnorm(10) + f1, rnorm(10)) m2 <- cbind(rnorm(10), rnorm(10) + f2) sds <- sqrt(0.2) n <- 100 xx1 <- xx2 <- matrix(0, n, 2) set.seed(5) t1<-t2<-rep(0,n) for (i in 1:n) { t1[i] <- sample(1:10, 1) xx1[i, ] <- c(rnorm(1, m1[t1[i], 1], sds), rnorm(1, m1[t1[i], 2], sds)) t2[i] <- sample(1:10, 1) xx2[i, ] <- c(rnorm(1, m2[t2[i], 1], sds), rnorm(1, m2[t2[i], 2], sds)) } xl1<-min(xx1[,1],xx2[,1]) xl2<-max(xx1[,1],xx2[,1]) yl1<-min(xx1[,2],xx2[,2]) yl2<-max(xx1[,2],xx2[,2]) plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri 6 −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] gnum <- 200 x1 <- seq(xl1, xl2, length.out = gnum) y1 <- seq(yl1, yl2, length.out = gnum) x <- rep(x1, gnum) #rep(c(1,2,3),10) y <- rep(y1, each = gnum) Question 2.1 # par(mfrow=c(3,1)); plot(dis1); plot(dis2); plot(dis1-dis2) dis1 <- dis2 <- numeric(length(y)) for (j in 1:10) { dis1 <- dis1 + dnorm(x, m1[j, 1], sds) * dnorm(y, m1[j, 2], sds) dis2 <- dis2 + dnorm(x, m2[j, 1], sds) * dnorm(y, m2[j, 2], sds) } z <- matrix(dis1 - dis2, gnum, gnum) contour(x1, y1, z, levels = 0, xlim = c(xl1, xl2), ylim = c(yl1, yl2), lwd = 2, lty = 2, col = "black") points(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri 7 0 −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 8 Question 2.2 plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri X <- rbind(xx1, xx2) Y <- factor(c(rep(-1, 100), rep(1, 100))) svmfit <- svm(X, Y, cost = 1.5, kernel = "polynomial", degree = 1, coef = 1) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] 9 Question 2.3 plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1.5, kernel = "polynomial", degree = 3, coef = 1) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] 10 Question 2.4 plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1.5, kernel = "radial", gamma = 1) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] 11 Question 2.5 polynomial with increased cost plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1500, kernel = "polynomial", degree = 3, coef = 1) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] 12 Question 2.5 radial with cost versus gamma plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1500, kernel = "radial", gamma = 1) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1.5, kernel = "radial", gamma = 50) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") 13 −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] plot(xx1[, 1], xx1[, 2], col = "dodgerblue3", pch = 3,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #plus points(xx2[, 1], xx2[, 2], col = "forestgreen", pch =24,xlim=c(xl1,xl2),ylim=c(yl1,yl2)) #tri svmfit <- svm(X, Y, cost = 1500, kernel = "radial", gamma = 50) zpred <- matrix(attr(predict(svmfit, cbind(x, y), decision.values = TRUE),"decision.values")[, 1], gnum, gnum) mycontour <- contourLines(x = x1, y = y1, zpred, levels = 0) lines(mycontour[[1]]$x, mycontour[[1]]$y, lwd = 2, lty = 1, col = "firebrick") 14 −3 −2 −1 0 1 2 3 4 − 3 − 2 − 1 0 1 2 3 xx1[, 1] xx 1 [, 2 ] A small cost (C) creates a soft margin and allows more misclassifications. On the other hand, a large cost creates a hard margin and permits fewer misclassifications. A large C gives you low bias and high variance, the bias is low because you have penalized the cost of missclasification a lot. For SVM with radial basis kernel, a high value of gamma leads to more accuracy in terms of classification. If gamma is high, the impact of C will become negligible. If gamma is weak, affect of C becomes significant 15 Libaries and least square estimate needed for Question 1 Question 1.1 Question 1.2 Question 1.3 Question 1.4 Libaries and data needed for Question 2 Question 2.1 Question 2.2 Question 2.3 Question 2.4 Question 2.5 polynomial with increased cost Question 2.5 radial with cost versus gamma A small cost (C) creates a soft margin and allows more misclassifications. On the other hand, a large cost creates a hard margin and permits fewer misclassifications. A large C gives you low bias and high variance, the bias is low because you have penalized the cost of missclasification a lot. For SVM with radial basis kernel, a high value of gamma leads to more accuracy in terms of classification. If gamma is high, the impact of C will become negligible. If gamma is weak, affect of C becomes significant