Functionals Excercise

Ex1. 為何兩個lapply呼叫結果相同

set.seed(100)
trims <-  c(0,0.1, 0.2,0.5)
x <- rcauchy(100)
lapply(trims, function(trim) mean(x, trim = trim) )
lapply(trims, mean, x = x)

1.第一種情況逐一將 trims 元素傳入呼叫 自行定義函式 function(trim)

lapply(trims, function(trim) mean(x, trim = trim) )
## [[1]]
## [1] 0.7701293
## 
## [[2]]
## [1] 0.1492742
## 
## [[3]]
## [1] 0.1315753
## 
## [[4]]
## [1] 0.03109642
mean(x, trim = 0)
## [1] 0.7701293
mean(x, trim = 0.5)
## [1] 0.03109642

2.第二種方式 lapply … 會將額外參數傳入呼叫 mean，因為 x 與 mean參數名稱相同，trims 每個元素被視為trim 設值

lapply(trims, mean, x = x)
## [[1]]
## [1] 0.7701293
## 
## [[2]]
## [1] 0.1492742
## 
## [[3]]
## [1] 0.1315753
## 
## [[4]]
## [1] 0.03109642
mean(0, x=x)
## [1] 0.7701293
mean(0.5, x=x)
## [1] 0.03109642

Ex.2 自行定義函式將向量數值正規化到[0,1]區間

scale01 <- function(x) {
        rng <- range(x, na.rm =  TRUE)
        (x-rng[1])/(rng[2] - rng[1])
}

如何將函式套用到 data frame的所有欄位

scale.mtcars <- sapply(mtcars, scale01)
head(scale.mtcars)
##            mpg cyl      disp        hp      drat        wt      qsec vs am
## [1,] 0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.2830478 0.2333333  0  1
## [2,] 0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.3482485 0.3000000  0  1
## [3,] 0.5276596 0.0 0.0920429 0.1448763 0.5023041 0.2063411 0.4892857  1  1
## [4,] 0.4680851 0.5 0.4662010 0.2049470 0.1474654 0.4351828 0.5880952  1  0
## [5,] 0.3531915 1.0 0.7206286 0.4346290 0.1797235 0.4927129 0.3000000  0  0
## [6,] 0.3276596 0.5 0.3838863 0.1872792 0.0000000 0.4978266 0.6809524  1  0
##      gear      carb
## [1,]  0.5 0.4285714
## [2,]  0.5 0.4285714
## [3,]  0.5 0.0000000
## [4,]  0.0 0.0000000
## [5,]  0.0 0.1428571
## [6,]  0.0 0.0000000

如何將函式套用到 data frame的所有Numeric欄位

numeric.cols <-  which(sapply(mtcars, class) %in% "numeric")
scale.numeric.mtcars <- sapply(mtcars[numeric.cols], scale01)
head(scale.numeric.mtcars)
##            mpg cyl      disp        hp      drat        wt      qsec vs am
## [1,] 0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.2830478 0.2333333  0  1
## [2,] 0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.3482485 0.3000000  0  1
## [3,] 0.5276596 0.0 0.0920429 0.1448763 0.5023041 0.2063411 0.4892857  1  1
## [4,] 0.4680851 0.5 0.4662010 0.2049470 0.1474654 0.4351828 0.5880952  1  0
## [5,] 0.3531915 1.0 0.7206286 0.4346290 0.1797235 0.4927129 0.3000000  0  0
## [6,] 0.3276596 0.5 0.3838863 0.1872792 0.0000000 0.4978266 0.6809524  1  0
##      gear      carb
## [1,]  0.5 0.4285714
## [2,]  0.5 0.4285714
## [3,]  0.5 0.0000000
## [4,]  0.0 0.0000000
## [5,]  0.0 0.1428571
## [6,]  0.0 0.0000000

Ex.3 Fit model with different foumula

fitMods <- lapply(formulas, lm, data = mtcars)

formulas <- list(
        mpg ~ disp,
        mpg ~ I(1/disp),
        mpg ~ disp + wt,
        mpg ~ I(1/disp) + wt
)

fitMods <- lapply(formulas, lm, data = mtcars)
fitMods[[4]]
## 
## Call:
## FUN(formula = X[[i]], data = ..1)
## 
## Coefficients:
## (Intercept)    I(1/disp)           wt  
##      19.024     1142.560       -1.798
fitMod4 <- lm(data=mtcars, formula= mpg ~ I(1/disp) + wt)
fitMod4
## 
## Call:
## lm(formula = mpg ~ I(1/disp) + wt, data = mtcars)
## 
## Coefficients:
## (Intercept)    I(1/disp)           wt  
##      19.024     1142.560       -1.798
rsq <- function(mod) summary(mod)$r.squared

lapply(fitMods, rsq)
## [[1]]
## [1] 0.7183433
## 
## [[2]]
## [1] 0.8596865
## 
## [[3]]
## [1] 0.7809306
## 
## [[4]]
## [1] 0.8838038
summary(fitMod4)$r.squared
## [1] 0.8838038

Ex.4 Fit model with different sample data frame

fitMods <- lapply(bootstraps, lm, formula = mpg ~ disp)

bootstraps <- lapply(1:10, function(i){
        rows <- sample(1:nrow(mtcars), rep=TRUE)
        mtcars[rows,]
})


fitMods <- lapply(bootstraps, lm, formula = mpg ~ disp)
fitMods[[10]]
## 
## Call:
## FUN(formula = ..1, data = X[[i]])
## 
## Coefficients:
## (Intercept)         disp  
##    26.96260     -0.03117
fitMod10 <- lm(data=bootstraps[[10]], formula= mpg ~ disp)
fitMod10
## 
## Call:
## lm(formula = mpg ~ disp, data = bootstraps[[10]])
## 
## Coefficients:
## (Intercept)         disp  
##    26.96260     -0.03117
lapply(fitMods, rsq)
## [[1]]
## [1] 0.8096264
## 
## [[2]]
## [1] 0.8056564
## 
## [[3]]
## [1] 0.6835524
## 
## [[4]]
## [1] 0.7359403
## 
## [[5]]
## [1] 0.7204606
## 
## [[6]]
## [1] 0.6629101
## 
## [[7]]
## [1] 0.7562879
## 
## [[8]]
## [1] 0.7705363
## 
## [[9]]
## [1] 0.767219
## 
## [[10]]
## [1] 0.6171997
summary(fitMod10)$r.squared
## [1] 0.6171997

lapply friends Exercise

Ex1.a 用vapply 計算數值 data frame每個欄位的標準差

vapply(mtcars, sd, numeric(1))
##         mpg         cyl        disp          hp        drat          wt 
##   6.0269481   1.7859216 123.9386938  68.5628685   0.5346787   0.9784574 
##        qsec          vs          am        gear        carb 
##   1.7869432   0.5040161   0.4989909   0.7378041   1.6152000
vapply(iris, sd, numeric(1))
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
##    0.8280661    0.4358663    1.7652982    0.7622377    0.8192319

Ex1.b 用vapply 計算 data frame中數值欄位的標準差

vapply(iris, class, class(1))
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
##    "numeric"    "numeric"    "numeric"    "numeric"     "factor"
num.idx <- which(vapply(iris, class, class(1)) %in% "numeric")
vapply(iris[,num.idx], sd, numeric(1))
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    0.8280661    0.4358663    1.7652982    0.7622377

Ex.3 下列程式碼對 non-normal 資料進行 t-test 用 sapply跟匿名函數取得每個測試的p-value

trials <- replicate(100, 
                    t.test(rpois(10,10), rpois(7,10)), 
                    simplify =  FALSE
)
trail.p.1 <- sapply(trials, function(t) {t$p.value} )
tail(trail.p.1,3)
## [1] 0.6610243 0.7206040 0.1067571
trail.p.2 <- sapply(1:100, function(i) {trials[[i]]["p.value"]} )
tail(trail.p.2,3)
## $p.value
## [1] 0.6610243
## 
## $p.value
## [1] 0.720604
## 
## $p.value
## [1] 0.1067571
trail.p.3 <-sapply(trials, `[[`, 'p.value')
tail(trail.p.3,3)
## [1] 0.6610243 0.7206040 0.1067571
trials[[100]]["p.value"]
## $p.value
## [1] 0.1067571

Ex.4 replicate 函式用途？他減少了那種類型的For loop?

函式參數與其他與lapply系列有何不同

?replicate

replicate 依照輸入次數，重複執行輸入運算式，通常搭配亂數產生不同運算結果

replicate is a wrapper for the common use of sapply for repeated evaluation of > an expression (which will usually involve random number generation).
沒傳入物件，所以應該是減少 數值指標類型的 for 迴圈 Three basic ways to loop over a vector
1. loop over the elements: for (x in xs)
2. loop over the numeric indices: for (i in seq_along(xs))
3. loop over the names: for (nm in names(xs))

lappply 依照傳入元素項，逐一元素應用傳入函式

lapply returns a list of the same length as X, each element of which is the > > result of applying FUN to the corresponding element of X.

Ex.5

Ex.6

試著實作 Map

matrix and data frame Excercise

Ex1. apply 函式如何安排結果？讀文件並做一些實驗

依照MARGIN 產算函式運算結果 - 1:by Rows - 2:by Column

apply(mtcars, 1, mean)
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##            29.90727            29.98136            23.59818 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##            38.73955            53.66455            35.04909 
##          Duster 360           Merc 240D            Merc 230 
##            59.72000            24.63455            27.23364 
##            Merc 280           Merc 280C          Merc 450SE 
##            31.86000            31.78727            46.43091 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##            46.50000            46.35000            66.23273 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##            66.05855            65.97227            19.44091 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##            17.74227            18.81409            24.88864 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##            47.24091            46.00773            58.75273 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##            57.37955            18.92864            24.77909 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##            24.88027            60.97182            34.50818 
##       Maserati Bora          Volvo 142E 
##            63.15545            26.26273
rowMeans(mtcars)
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##            29.90727            29.98136            23.59818 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##            38.73955            53.66455            35.04909 
##          Duster 360           Merc 240D            Merc 230 
##            59.72000            24.63455            27.23364 
##            Merc 280           Merc 280C          Merc 450SE 
##            31.86000            31.78727            46.43091 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##            46.50000            46.35000            66.23273 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##            66.05855            65.97227            19.44091 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##            17.74227            18.81409            24.88864 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##            47.24091            46.00773            58.75273 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##            57.37955            18.92864            24.77909 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##            24.88027            60.97182            34.50818 
##       Maserati Bora          Volvo 142E 
##            63.15545            26.26273
apply(mtcars, 2, mean)
##        mpg        cyl       disp         hp       drat         wt 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250 
##       qsec         vs         am       gear       carb 
##  17.848750   0.437500   0.406250   3.687500   2.812500
colMeans(mtcars)
##        mpg        cyl       disp         hp       drat         wt 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250 
##       qsec         vs         am       gear       carb 
##  17.848750   0.437500   0.406250   3.687500   2.812500

Ex.2 split() + vapply()，有需要這樣的整合應嗎？

lapply(split(mtcars, mtcars$cyl), colSums)
## $`4`
##      mpg      cyl     disp       hp     drat       wt     qsec       vs 
##  293.300   44.000 1156.500  909.000   44.780   25.143  210.510   10.000 
##       am     gear     carb 
##    8.000   45.000   17.000 
## 
## $`6`
##     mpg     cyl    disp      hp    drat      wt    qsec      vs      am 
##  138.20   42.00 1283.20  856.00   25.10   21.82  125.84    4.00    3.00 
##    gear    carb 
##   27.00   24.00 
## 
## $`8`
##      mpg      cyl     disp       hp     drat       wt     qsec       vs 
##  211.400  112.000 4943.400 2929.000   45.210   55.989  234.810    0.000 
##       am     gear     carb 
##    2.000   46.000   49.000
vapply(split(mtcars, mtcars$cyl), colSums, numeric(11))
##             4       6        8
## mpg   293.300  138.20  211.400
## cyl    44.000   42.00  112.000
## disp 1156.500 1283.20 4943.400
## hp    909.000  856.00 2929.000
## drat   44.780   25.10   45.210
## wt     25.143   21.82   55.989
## qsec  210.510  125.84  234.810
## vs     10.000    4.00    0.000
## am      8.000    3.00    2.000
## gear   45.000   27.00   46.000
## carb   17.000   24.00   49.000

The Split-Apply-Combine Strategy for Data Analysis

library(plyr)
aq<-airquality
aq$source<-"qiuworld"
aq1<-aq
aq1$source<-c("gaziou")
set.seed(123)
x1<-runif(nrow(aq1),-0.5,0.5)
aq1$Temp<-aq1$Temp+x1
aq<-rbind(aq,aq1)

sapply(split(aq, list(aq$source, aq$Month)), 
       function(x) {mean(x$Temp)})
##   gaziou.5 qiuworld.5   gaziou.6 qiuworld.6   gaziou.7 qiuworld.7 
##   65.63339   65.54839   79.02871   79.10000   83.90313   83.90323 
##   gaziou.8 qiuworld.8   gaziou.9 qiuworld.9 
##   83.98611   83.96774   76.89319   76.90000
#use ~ source + Month define variables 
avgTemp<-ddply(aq, ~ source + Month,function(a) mean(a$Temp))
avgTemp
##      source Month       V1
## 1    gaziou     5 65.63339
## 2    gaziou     6 79.02871
## 3    gaziou     7 83.90313
## 4    gaziou     8 83.98611
## 5    gaziou     9 76.89319
## 6  qiuworld     5 65.54839
## 7  qiuworld     6 79.10000
## 8  qiuworld     7 83.90323
## 9  qiuworld     8 83.96774
## 10 qiuworld     9 76.90000
#use .(source,Month) define variables 
avgWind<-ddply(aq, .(source,Month),function(a) mean(a$Wind))
avgWind
##      source Month        V1
## 1    gaziou     5 11.622581
## 2    gaziou     6 10.266667
## 3    gaziou     7  8.941935
## 4    gaziou     8  8.793548
## 5    gaziou     9 10.180000
## 6  qiuworld     5 11.622581
## 7  qiuworld     6 10.266667
## 8  qiuworld     7  8.941935
## 9  qiuworld     8  8.793548
## 10 qiuworld     9 10.180000
suppressMessages(library(dplyr)) 
aq.wind.temp <- aq %>% 
        group_by(source, Month) %>%
        summarise(mean.wind = mean(Wind),
                  mean.Temp = mean(Temp))
aq.wind.temp
## Source: local data frame [10 x 4]
## Groups: source
## 
##      source Month mean.wind mean.Temp
## 1    gaziou     5 11.622581  65.63339
## 2    gaziou     6 10.266667  79.02871
## 3    gaziou     7  8.941935  83.90313
## 4    gaziou     8  8.793548  83.98611
## 5    gaziou     9 10.180000  76.89319
## 6  qiuworld     5 11.622581  65.54839
## 7  qiuworld     6 10.266667  79.10000
## 8  qiuworld     7  8.941935  83.90323
## 9  qiuworld     8  8.793548  83.96774
## 10 qiuworld     9 10.180000  76.90000

Ex3. 用Ｒ實作 Split 提示用 unique + subseting

split2 <- function(x, f) {
        groups <- unique(f)
        groups <- groups[order(groups)]
        out <- lapply(groups, function(var) { 
                x[which(f == var),]})
        names(out) <- levels(groups)
        return (out)
}
split2(mtcars, mtcars$gear)
## [[1]]
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## 
## [[2]]
##                 mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4      21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710     22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230       22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280       19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C      17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Fiat 128       32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic    30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla 33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Fiat X1-9      27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Volvo 142E     21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2
## 
## [[3]]
##                 mpg cyl  disp  hp drat    wt qsec vs am gear carb
## Porsche 914-2  26.0   4 120.3  91 4.43 2.140 16.7  0  1    5    2
## Lotus Europa   30.4   4  95.1 113 3.77 1.513 16.9  1  1    5    2
## Ford Pantera L 15.8   8 351.0 264 4.22 3.170 14.5  0  1    5    4
## Ferrari Dino   19.7   6 145.0 175 3.62 2.770 15.5  0  1    5    6
## Maserati Bora  15.0   8 301.0 335 3.54 3.570 14.6  0  1    5    8

List 操作 Exercises

Ex1. Why isn’t is.na() a predicate function? What base R function is closest to being a predicate version of is.na()?

Filter(x=iris, f= anyNA)
## data frame with 0 columns and 150 rows
Filter(x=iris, f= is.na)
## data frame with 0 columns and 150 rows
iris[151, 1] <- 5
iris[151, 4] <- 5
iris[151, 5] <- "virginica"

Ex2. Use Filter() and vapply() to create a function that applies a summary statistic to every numeric column in a data frame.

summaryNumeric <- function(data) {
        num.data <- Filter(x=data, f=is.numeric)
        na.data <- Filter(num.data, f=anyNA)
        v1 <- vapply(na.data, summary, vector(mode="numeric", length=7))   
        
        no.na.data <- Filter(num.data, f=Negate(anyNA))
        v2 <- vapply(no.na.data, summary, vector(mode="numeric", length=6)) 
        v2 <- rbind(v2, NA)
        cbind( v2, v1)
}
tail(iris, 3)
##     Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
## 149          6.2         3.4          5.4         2.3 virginica
## 150          5.9         3.0          5.1         1.8 virginica
## 151          5.0          NA           NA         5.0 virginica
summaryNumeric(iris)
##         Sepal.Length Petal.Width Sepal.Width Petal.Length
## Min.           4.300       0.100       2.000        1.000
## 1st Qu.        5.100       0.300       2.800        1.600
## Median         5.800       1.300       3.000        4.350
## Mean           5.838       1.225       3.057        3.758
## 3rd Qu.        6.400       1.800       3.300        5.100
## Max.           7.900       5.000       4.400        6.900
##                   NA          NA       1.000        1.000

Ex3. What’s the relationship between which() and Position()? What’s the relationship between where() and Filter()?

Ex4. Implement Any(), a function that takes a list and a predicate function, and returns TRUE if the predicate function returns TRUE for any of the inputs. Implement All() similarly.

Ex.5 Implement the span() function from Haskell: given a list x and a predicate function f, span returns the location of the longest sequential run of elements where the predicate is true. (Hint: you might find rle() helpful.)

練習實作 subset

subset2 <- function(x, condition){
        condition_call <- substitute(condition)
        r <- eval(condition_call, x)
        x[r, ]
}
subset2(mtcars, cyl==6)
##                 mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4      21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Valiant        18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Merc 280       19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C      17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Ferrari Dino   19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6

rollapply 參數說明

  • by: calculate FUN at every by-th time point rather than every point. by is only used if width is length 1 and either a plain scalar or a list. （每組計算後間隔幾筆資料再取值）

  • window: numeric vector or list. In the simplest case this is an integer specifying the window width (in numbers of observations) which is aligned to the original sample according to the align argument. Alternatively, width can be a list regarded as offsets compared to the current time, see below for details. (每次取幾筆資料計算)

  • align: 決定日期切齊方向

require(zoo)

產生 取得6個數值的 zoo 物件

z <- zoo(11:16, as.Date(31:36))
z
## 1970-02-01 1970-02-02 1970-02-03 1970-02-04 1970-02-05 1970-02-06 
##         11         12         13         14         15         16

每兩個取平均值，日期切齊左邊

rollapply(data = z, width = 2, mean, align = "left")
## 1970-02-01 1970-02-02 1970-02-03 1970-02-04 1970-02-05 
##       11.5       12.5       13.5       14.5       15.5

每兩個取平均值，每次間隔兩個日期，日期切齊左邊

rollapply(data = z, width = 2, mean, align = "left", by = 2)
## 1970-02-01 1970-02-03 1970-02-05 
##       11.5       13.5       15.5

每兩個取平均值，日期切齊右邊

rollapply(data = z, width = 2, mean, align = "right")
## 1970-02-02 1970-02-03 1970-02-04 1970-02-05 1970-02-06 
##       11.5       12.5       13.5       14.5       15.5

每兩個取平均值，每次間隔兩個日期，日期切齊右邊

rollapply(data = z, width = 2, mean, align = "right", by = 2)
## 1970-02-02 1970-02-04 1970-02-06 
##       11.5       13.5       15.5

每3個取平均值，間隔4個日期，間隔四天後，只剩兩組資料，只有1筆結果，日期切齊左邊

rollapply(data = z, width = 3, mean, by = 4, align = "left")
## 1970-02-01 
##         12

每3個取平均值，間隔4個日期，間隔四天後，只剩兩組資料，只有1筆結果，日期切齊右邊

rollapply(data = z, width = 3, mean, by = 4, align = "right")
## 1970-02-03 
##         12

每2個取平均值，間隔4個日期，間隔四天後，足夠兩筆資料計算平均，會有2筆結果，日期切齊左邊

rollapply(data = z, width = 2, mean, by = 4, align = "left")
## 1970-02-01 1970-02-05 
##       11.5       15.5

每2個取平均值，間隔4個日期，間隔四天後，足夠兩筆資料計算平均，會有2筆結果，日期切齊右邊

rollapply(data = z, width = 2, mean, by = 4, align = "right")
## 1970-02-02 1970-02-06 
##       11.5       15.5

Homework 11/02

定義 Function 計算 歐式買權, 歐式賣權

require(fOptions)

[Black Sholes Equation](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_equation_ -TypeFlag 選擇權類型 - C︰買權 (call) - P︰賣權 (put)

-ncdf︰累加常態分配函數 (cumulative normal distribution function)

-S︰現貨價格 (spot price)

-K︰履約價格 (strike price)(exercise price)

-r︰無風險利率 (risk-free rate)，通常會用相同天期的公債殖利率

-deltaT︰離到期日還有多久，單位是年

-sigma(σ)︰波動率 (volatility)，單位是 %

MyBlackScholeOptionPrice <- function (TypeFlag, S, K, deltaT, r, sigma, ncdf = pnorm) {
        d1 = (log(S/K) + (r + 0.5 * (sigma ^2 ) * deltaT ))/(sigma * sqrt(deltaT))
        d2 = d1 - sigma * sqrt(deltaT)
        
        # Calculate Call Price
        if(grepl("C", TypeFlag, ignore.case = TRUE)) {
                price <- S * ncdf(d1) - 
                         K * exp(-r *deltaT) * ncdf(d2) 
        } else if(grepl("P", TypeFlag, ignore.case = TRUE)) {
                # Calculate Put Price                                                                     
                price = K * exp(-r * deltaT) * ncdf(-d2) - S * ncdf(-d1)
        } else {
                stop("Invalid TypeFlag, should be Call or Put")
        }
        return (price)

}
myCallPrice <- MyBlackScholeOptionPrice(TypeFlag = "c", S = 2.8, K = 2.9, deltaT= 0.083333333, r = 0.01, sigma = 0.2, ncdf = pnorm)
myCallPrice
## [1] 0.02729898
fOpCallPrice <- GBSOption(TypeFlag = "c", S=2.8, X=2.9, Time= 0.083333333, r = 0.01, sigma = 0.2, b= 0)

fOpCallPrice@price
## [1] 0.02737078
myPutPrice <- MyBlackScholeOptionPrice(TypeFlag = "p", S = 2.8, K = 2.9, deltaT= 0.083333333, r = 0.01, sigma = 0.2, ncdf = pnorm)
myPutPrice
## [1] 0.1248833
fOpPutPrice <- GBSOption(TypeFlag = "p", S=2.8, X=2.9, Time= 0.083333333, r = 0.01, sigma = 0.2, b= 0)
fOpPutPrice@price
## [1] 0.1272875