EX01

Exchange <- function(){
  money <- readline(prompt = "Enter NT$: ")
  cat("Got ", as.numeric(money)/29.31, " US$\n", sep = "")
}
Exchange()
## Enter NT$: 
## Got NA US$

EX02

dta02 <- ChickWeight
sapply(split(dta02, dta02$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

EX03

t2z <- function(df){
  curve(dnorm(x), -4, 4, col = 2, ylab = "dnorm(x)", lwd = 2)
  for(i in 1:length(df)) curve(dt(x, df[[i]]), col = 3, lty = 2, add = TRUE)
}
t2z(df = 1:32)

EX04

library(pacman)
pacman::p_load(MASS, tidyverse)

## method 1, 簡單的dataframe形式
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
## method 2, wide format
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
## method 3, long format
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
## method 4, 修改了原本的data
Cushings %>%
 group_by(Type) %>%
 summarize( t_m = mean(Tetrahydrocortisone), p_m = mean(Pregnanetriol))
## # A tibble: 4 x 3
##   Type    t_m   p_m
##   <fct> <dbl> <dbl>
## 1 a      2.97  2.44
## 2 b      8.18  1.12
## 3 c     19.7   5.50
## 4 u     14.0   1.20
## method 5, 在縮減資料的同時保留了原本的data，提供資料屬性。
Cushings %>%
 nest(-Type) %>%
 mutate(avg = map(data, ~ apply(., 2, mean)), 
        res_1 = map_dbl(avg, "Tetrahydrocortisone"), 
        res_2 = map_dbl(avg, "Pregnanetriol"))
## # 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

EX05

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()
}

# run
lawLN(4000, 100, 10)

EX06

dta06 <- read.table("cstat.txt", header=T)
# function
c.stat <- function(data, n = length(data)){
  cden <- 1-(sum(diff(dta06[1:n,1])^2)/(2*(n-1)*var(dta06[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))
}

c.stat(dta06, 42) 
## $c
## [1] 0.6450652
## 
## $z
## [1] 4.282524
## 
## $pvalue
## [1] 9.239272e-06

EX07

library(rgl)
knitr::knit_hooks$set(webgl = hook_webgl)

# data
dta07 <- read.table("hs0.txt", header=TRUE)   
math_score <- dta07$math

# define ssq and vssq functions
ssq <- function(mu, sigma, y) {sum(((y - mu) / sigma)^2)}
vssq <- Vectorize(ssq, c("mu", "sigma"))

# find minimun of the ssq and the sigma
math_mu = sum(math_score)/length(math_score)
math_sigma = sqrt(sum((math_score - math_mu)^2) / (length(math_score) - 1))
math_ssq = ssq(math_mu, math_sigma, math_score)

# calculate ssq for several mu and sigma
my_mu <- seq(math_mu - 5, math_mu + 5,by=0.2)
my_sigma <- math_sigma
vssq_result <- vssq(my_mu, my_sigma, math_score)

# visualized
open3d()
## wgl 
##   1
plot3d(my_mu, my_sigma, vssq_result,
       xlab = "mu", ylab = "sigma", zlab = "ssq",
       col = rainbow(length(my_mu)))

# log
cat('The ideal minimum of ssq is ', math_ssq, '.\n', sep='')
## The ideal minimum of ssq is 199.
cat('The minimum in vssq is ', min(vssq_result), '.', sep='')
## The minimum in vssq is 199.

You must enable Javascript to view this page properly.