library(coda)
## Warning: package 'coda' was built under R version 4.0.3
library(MASS)
## Warning: package 'MASS' was built under R version 4.0.3
library(mcsm)
## Warning: package 'mcsm' was built under R version 4.0.3
demo(Chapter.1)
##
##
## demo(Chapter.1)
## ---- ~~~~~~~~~
##
## > # Chapter 1 R commands
## >
## > # Section 1.3.2
## >
## > x=matrix(1:4,ncol=3)
## Warning in matrix(1:4, ncol = 3): la longitud de los datos [4] no es un
## submúltiplo o múltiplo del número de columnas [3] en la matriz
##
## > print(x[x>5])
## integer(0)
##
## > print(x[1.])
## [1] 1
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > # Section 1.3.3
## >
## > x = list(a = 1:10, beta = exp(-3:3),
## + logic = c(TRUE,FALSE,FALSE,TRUE))
##
## > print(lapply(x,mean))
## $a
## [1] 5.5
##
## $beta
## [1] 4.535125
##
## $logic
## [1] 0.5
##
##
## > print(sapply(x,mean))
## a beta logic
## 5.500000 4.535125 0.500000
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > # Section 1.4
## >
## > x=rnorm(20)
##
## > y=3*x+5+rnorm(20,sd=0.3)
##
## > reslm=lm(y~x)
##
## > print(summary(reslm))
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36944 -0.12842 -0.00643 0.14187 0.42563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.02797 0.05656 88.90 <2e-16 ***
## x 3.00450 0.05426 55.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2447 on 18 degrees of freedom
## Multiple R-squared: 0.9942, Adjusted R-squared: 0.9938
## F-statistic: 3066 on 1 and 18 DF, p-value: < 2.2e-16
##
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > # Section 1.5
## >
## > b=matrix(1:9,ncol=3)
##
## > print(var(b))
## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 1 1 1
## [3,] 1 1 1
##
## > print(sd(b)^2)
## [1] 7.5
##
## > x=rnorm(25) #produces a N(0,1) sample of size 25
##
## > out=t.test(x)
##
## > print(names(out))
## [1] "statistic" "parameter" "p.value" "conf.int" "estimate"
## [6] "null.value" "stderr" "alternative" "method" "data.name"
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > attach(faithful) #resident dataset
##
## > print(cor.test(faithful[,1],faithful[,2]))
##
## Pearson's product-moment correlation
##
## data: faithful[, 1] and faithful[, 2]
## t = 34.089, df = 270, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8756964 0.9210652
## sample estimates:
## cor
## 0.9008112
##
##
## > print(ks.test(jitter(faithful[,1]),pnorm))
##
## One-sample Kolmogorov-Smirnov test
##
## data: jitter(faithful[, 1])
## D = 0.94858, p-value < 2.2e-16
## alternative hypothesis: two-sided
##
##
## > print(shapiro.test(faithful[,2]))
##
## Shapiro-Wilk normality test
##
## data: faithful[, 2]
## W = 0.92215, p-value = 1.015e-10
##
##
## > print(wilcox.test(faithful[,1]))
##
## Wilcoxon signed rank test with continuity correction
##
## data: faithful[, 1]
## V = 37128, p-value < 2.2e-16
## alternative hypothesis: true location is not equal to 0
##
##
## > #S=readline(prompt="Type <Return> to continue : ")
## >
## > x=seq(-3,3,le=5) # equidispersed regressor
##
## > y=2+4*x+rnorm(5) # simulated variable
##
## > print(lm(y~x))
##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 2.748 4.156
##
##
## > print(summary(lm(y~x)))
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## 1 2 3 4 5
## -0.7342 0.4557 1.0177 -0.4659 -0.2734
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7476 0.3719 7.389 0.005126 **
## x 4.1562 0.1753 23.710 0.000164 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8315 on 3 degrees of freedom
## Multiple R-squared: 0.9947, Adjusted R-squared: 0.9929
## F-statistic: 562.2 on 1 and 3 DF, p-value: 0.0001644
##
##
## > out=lm(y~x)
##
## > print(sqrt(sum(out$res^2)/out$df))
## [1] 0.8314914
##
## > #S=readline(prompt="Type <Return> to continue : ")
## >
## > print(summary(lm(weight ~ feed, data = chickwts)))
##
## Call:
## lm(formula = weight ~ feed, data = chickwts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -123.909 -34.413 1.571 38.170 103.091
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 323.583 15.834 20.436 < 2e-16 ***
## feedhorsebean -163.383 23.485 -6.957 2.07e-09 ***
## feedlinseed -104.833 22.393 -4.682 1.49e-05 ***
## feedmeatmeal -46.674 22.896 -2.039 0.045567 *
## feedsoybean -77.155 21.578 -3.576 0.000665 ***
## feedsunflower 5.333 22.393 0.238 0.812495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 54.85 on 65 degrees of freedom
## Multiple R-squared: 0.5417, Adjusted R-squared: 0.5064
## F-statistic: 15.36 on 5 and 65 DF, p-value: 5.936e-10
##
##
## > print(anova(lm(weight ~ feed, data = chickwts)))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## feed 5 231129 46226 15.365 5.936e-10 ***
## Residuals 65 195556 3009
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > library(MASS)
##
## > print(glm(formula = type ~ bmi + age, family = "binomial", data = Pima.tr))
##
## Call: glm(formula = type ~ bmi + age, family = "binomial", data = Pima.tr)
##
## Coefficients:
## (Intercept) bmi age
## -6.49870 0.10519 0.07104
##
## Degrees of Freedom: 199 Total (i.e. Null); 197 Residual
## Null Deviance: 256.4
## Residual Deviance: 215.9 AIC: 221.9
##
## > print(arima(diff(EuStockMarkets[,1]),order=c(0,0,5)))
##
## Call:
## arima(x = diff(EuStockMarkets[, 1]), order = c(0, 0, 5))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 intercept
## 0.0054 -0.0130 -0.0110 -0.0041 -0.0486 2.0692
## s.e. 0.0234 0.0233 0.0221 0.0236 0.0235 0.6990
##
## sigma^2 estimated as 1053: log likelihood = -9106.23, aic = 18226.45
##
## > print(acf(ldeaths, plot=F))
##
## Autocorrelations of series 'ldeaths', by lag
##
## 0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333
## 1.000 0.755 0.397 0.019 -0.356 -0.609 -0.681 -0.608 -0.378 -0.013 0.383
## 0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000
## 0.650 0.723 0.638 0.372 0.009 -0.294 -0.497 -0.586
##
## > print(acf(ldeaths))
##
## Autocorrelations of series 'ldeaths', by lag
##
## 0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333
## 1.000 0.755 0.397 0.019 -0.356 -0.609 -0.681 -0.608 -0.378 -0.013 0.383
## 0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000
## 0.650 0.723 0.638 0.372 0.009 -0.294 -0.497 -0.586
##
## > print(acf(ldeaths,type="partial"))
##
## Partial autocorrelations of series 'ldeaths', by lag
##
## 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 0.9167
## 0.755 -0.403 -0.269 -0.372 -0.192 -0.116 -0.153 0.005 0.202 0.295 0.121
## 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000
## 0.022 0.138 -0.027 -0.061 0.094 0.134 0.051
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > # Section 1.6
## >
## > plot(faithful)
##
## > jpeg(file="faith.jpg")
##
## > par(mfrow=c(1,2),mar=c(4,2,2,1))
##
## > hist(faithful[,1],nclass=21,col="grey",main="",xlab=names(faithful)[1])
##
## > hist(faithful[,2],nclass=21,col="wheat",main="",xlab=names(faithful)[2])
##
## > dev.off()
## png
## 2
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > plot(as.vector(time(mdeaths)),as.vector(mdeaths),cex=.6,pch=19,xlab="",ylab="Monthly deaths from bronchitis in the UK")
##
## > lines(spline(mdeaths),lwd=2,col="chocolate",lty=3)
##
## > ar=arima(mdeaths,order=c(1,0,0))$coef
##
## > lines(as.vector(time(mdeaths))[-1], ar[2]+ar[1]*(mdeaths[-length(mdeaths)]-ar[2]),col="grey",lwd=2,lty=2)
##
## > title("Splines versus AR(1) predictive")
##
## > ar=arima(mdeaths,order=c(1,0,0),seasonal=list(order=c(1,0,0),period=12))$coef
##
## > lines(as.vector(time(mdeaths))[-(1:13)],ar[3]+ar[1]*(mdeaths[-c(1:12,72)]-ar[3])+ar[2]*(mdeaths[-(60:72)]-ar[3]),
## + lwd=2,col="steelblue",lty=2)
##
## > title("\n\nand SAR(1,12) predictive")
##
## > legend(1974,2800,legend=c("spline","AR(1)","SAR(1,12)"),col=c("chocolate","grey","steelblue"),
## + lty=c(3,2,2),lwd=rep(2,3),cex=.5)
##
## > S=readline(prompt="Type <Return> to continue : ")
## Type <Return> to continue :
##
## > # Section 1.7
## >
## > sqrnt=function(y) {
## +
## + x=y/2
## + while (abs(x*x-y) > 1e-10) x=(x+y/x)/2
## + x
## + }
##
## > sqrnt(9)
## [1] 3
##
## > sqrnt(2)
## [1] 1.414214