pacman::p_load(tidyr, broom, magrittr, dplyr, ggplot2,
nlme, kableExtra, furniture, purrr)
1
converter <- function( ) {
money <- readline(prompt="Enter the number of NT dollars you want to convert: ")
cat("You can convert your NT dollars to $", as.numeric(money)/29.32, " US dollars\n", sep = "")
}
converter()
## Enter the number of NT dollars you want to convert:
## You can convert your NT dollars to $NA US dollars
2
dta <- ChickWeight
sapply(split(dta, dta$Chick),
function(x) lm(weight ~ Time, data = x)$coef)
## 18 16 15 13 9 20 10
## (Intercept) 39 43.392857 46.83333 43.384359 52.094086 37.667826 38.695054
## Time -2 1.053571 1.89881 2.239601 2.663137 3.732718 4.066102
## 8 17 19 4 6 11
## (Intercept) 43.727273 43.030706 31.21222 32.86568 44.123431 47.921948
## Time 4.827273 4.531538 5.08743 6.08864 6.378006 7.510967
## 3 1 12 2 5 14
## (Intercept) 23.17955 24.465436 21.939797 24.724853 16.89563 20.52488
## Time 8.48737 7.987899 8.440629 8.719861 10.05536 11.98245
## 7 24 30 22 23 27
## (Intercept) 5.842535 53.067766 39.109666 40.082590 38.428074 29.858569
## Time 13.205264 1.207533 5.898351 5.877931 6.685978 7.379368
## 28 26 25 29 21 33
## (Intercept) 23.984874 20.70715 19.65119 5.882771 15.56330 45.830283
## Time 9.703676 10.10316 11.30676 12.453487 15.47512 5.855241
## 37 36 31 39 38 32
## (Intercept) 29.608834 25.85403 19.13099 17.03661 10.67282 13.69173
## Time 6.677053 9.99047 10.02617 10.73710 12.06051 13.18091
## 40 34 35 44 45 43
## (Intercept) 10.83830 5.081682 4.757979 44.909091 35.673121 52.185751
## Time 13.44229 15.000151 17.258811 6.354545 7.686432 8.318863
## 41 47 49 46 50 42
## (Intercept) 39.337922 36.489790 31.662986 27.771744 23.78218 19.86507
## Time 8.159885 8.374981 9.717894 9.738466 11.33293 11.83679
## 48
## (Intercept) 7.947663
## Time 13.714718
3
x <- seq(-4, 4, .05)
df <- c(2, 4, 8, 16, 32)
plot(x, dnorm(x), type = 'l', ylab = "dnorm(x)", col = "red")
for(i in 1:5){lines(x, dt(x, df[i]), lty = 2, col = "green")}
play <- function(){
dff <- readline(prompt = "Enter a number")
plot(x, dnorm(x), type = 'l', ylab = "dnorm(x)", col = "red")
lines(x, dt(x, as.numeric(dff)), lty = 2, col = "green")
legend("topright",
legend = c("Normal", "df = you typed"), lty = 2, col = c("red", "green"))
}
play()
## Enter a number
4
library(MASS)
## Warning: package 'MASS' was built under R version 3.4.4
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
str(Cushings)
## 'data.frame': 27 obs. of 3 variables:
## $ Tetrahydrocortisone: num 3.1 3 1.9 3.8 4.1 1.9 8.3 3.8 3.9 7.8 ...
## $ Pregnanetriol : num 11.7 1.3 0.1 0.04 1.1 0.4 1 0.2 0.6 1.2 ...
## $ Type : Factor w/ 4 levels "a","b","c","u": 1 1 1 1 1 1 2 2 2 2 ...
head(Cushings)
## Tetrahydrocortisone Pregnanetriol Type
## a1 3.1 11.70 a
## a2 3.0 1.30 a
## a3 1.9 0.10 a
## a4 3.8 0.04 a
## a5 4.1 1.10 a
## a6 1.9 0.40 a
# 1
#內容較精簡
aggregate( . ~ Type, data = Cushings, mean)
## Type Tetrahydrocortisone Pregnanetriol
## 1 a 2.966667 2.44
## 2 b 8.180000 1.12
## 3 c 19.720000 5.50
## 4 u 14.016667 1.20
# 2
# 排列方式和其他方法不同
sapply(split(Cushings[,-3], Cushings$Type), function(x) apply(x, 2, mean))
## a b c u
## Tetrahydrocortisone 2.966667 8.18 19.72 14.01667
## Pregnanetriol 2.440000 1.12 5.50 1.20000
# 3
# 相對複雜
do.call("rbind", as.list(
by(Cushings, list(Cushings$Type), function(x) {
y <- subset(x, select = -Type)
apply(y, 2, mean)
}
)))
## Tetrahydrocortisone Pregnanetriol
## a 2.966667 2.44
## b 8.180000 1.12
## c 19.720000 5.50
## u 14.016667 1.20
# 5
# 內容較難懂、提供資料屬性
Cushings %>%
nest(-Type) %>%
mutate(avg = map(data, ~ apply(., 2, mean)),
res_1 = map_dbl(avg, "Tetrahydrocortisone"),
res_2 = map_dbl(avg, "Pregnanetriol"))
## Warning: package 'bindrcpp' was built under R version 3.4.4
## # A tibble: 4 x 5
## Type data avg res_1 res_2
## <fct> <list> <list> <dbl> <dbl>
## 1 a <data.frame [6 x 2]> <dbl [2]> 2.97 2.44
## 2 b <data.frame [10 x 2]> <dbl [2]> 8.18 1.12
## 3 c <data.frame [5 x 2]> <dbl [2]> 19.7 5.50
## 4 u <data.frame [6 x 2]> <dbl [2]> 14.0 1.20
5
lawLN <- function(n, mu, s){
set.seed(0221)
random.sample <- rnorm(n, mu, s)
plot(x = 1:n, y = cumsum(random.sample)/1:n, type = "l", col = 3,
xlab = "Sample Size", ylab = "Running Average")
abline(h = mu, col = 2, lty = 2)
grid()
}
lawLN(4000, 100, 10)
6
dta <- read.table("/Users/user/Documents/cstat.txt", header = TRUE)
# function
c.stat <- function(data, n = length(data)){
cden <- 1-(sum(diff(dta[1:n,1])^2)/(2*(n-1)*var(dta[1:n,1])))
sc <- sqrt((n-2)/((n-1)*(n+1)))
pval <- 1-pnorm(cden/sc)
return(list(c = cden, z = cden/sc, pvalue = pval))
}
# run
c.stat(dta, 42)
## $c
## [1] 0.6450652
##
## $z
## [1] 4.282524
##
## $pvalue
## [1] 9.239272e-06
7
# data
dta <- read.table("/Users/user/Documents/hs0.txt", header = TRUE)
# ssq function in lecture
ssq <- function(mu, sigma, y) {sum(((y - mu) / sigma)^2)}
# vectorize function
vssq <- Vectorize(ssq, c("mu", "sigma"))
# minimal math vssq
m_vssq <- with(dta, vssq(mu = quantile(math, 1/4):quantile(math, 3/4),
sigma = var(math), y = math))
f <- function(data = NA){
z <- vssq(mu = quantile(data, 1/4):quantile(data, 3/4),
sigma = var(data),
y = data)
require(plot3D)
image3D(x = quantile(data, 1/4):quantile(data, 3/4),
y = var(data),
z = data,
xlab = "mu", ylab = "sigma", zlab = "score",
col = "steelblue")
}
f(data = dta$math)
## Loading required package: plot3D
## Warning: package 'plot3D' was built under R version 3.4.4
