4/9 HomeWork

HW1

# Make a convert function
NTtoUS <- function(){
  money <- readline(prompt = "Enter your NT$ here:")
  cat("You've got $", as.numeric(money)/29.32, " US dollars\n", sep = "")
}
NTtoUS()

## Enter your NT$ here:
## You've got $NA US dollars

HW2

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

HW3

# Transfer t to z
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)

# HW4

# method 1:比較單純的dataframe,輸出形式為list

m1 <- aggregate( . ~ Type, data = Cushings, mean)

# method 2:形式為矩陣，輸出形式為為數值，並且是wide format

m2 <- sapply(split(Cushings[,-3], Cushings$Type), function(x) apply(x, 2, mean))

# method 3:形式同樣為矩陣，但這次是long format
m3 <-do.call("rbind", as.list(
  by(Cushings, list(Cushings$Type), function(x) {
    y <- subset(x, select =  -Type)
    apply(y, 2, mean)
  }
)))


# method 4:以pipe作為分群與連接方式，修改了原本的data，輸出為list
m4 <-Cushings %>%
 group_by(Type) %>%
 summarize( t_m = mean(Tetrahydrocortisone), p_m = mean(Pregnanetriol))

# method 5:類似的pipe語法，但同時給了修改後與原本的資料

m5 <- 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

###

HW5

# set a function that can run a plot ramdonmly
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)

HW6

dta6 <- read.table("cstat.txt", header = TRUE)

# function
c.stat <- function(data, n = length(data)){
  cden <- 1-(sum(diff(dta6[1:n,1])^2)/(2*(n-1)*var(dta6[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(dta6, 42) 

## $c
## [1] 0.6450652
## 
## $z
## [1] 4.282524
## 
## $pvalue
## [1] 9.239272e-06