2017年11月2日 星期四

2017 R TAIWAN 研討會



會議名稱:2017 R TAIWAN 研討會

會議主軸:資料融合、分析與應用 (不限定R語言)
會議地址:台北市中正區貴陽街一段56號 (東吳大學城中校區)
會議教室:東吳大學城中校區 5211演講廳、5117教室、2123教室

宗旨:
Big Data 是近年來熱門的話題之一,根據 2016年 KDnuggets 的調查顯示,R仍然高居資料解析、資料探勘與資料科學主要程式語言的首位。但R絕對不是萬能無缺的!搭配其它資料處理與分析的工具,方能捕捉到 Big一字的意涵 – 從關鍵的數據中解析真正不同且重要的洞見。依稀記得補習班常用的一句廣告詞:好的老師帶你上天堂,平庸的老師讓你原地踏步。R 無疑是個關鍵且能造血活血的工具,包羅萬象的R套件可以持續提升我們的資料解析能力。除非妳(你)想安逸,否則會被它一步一步地往前推。人人都希望經常遇見好的老師,但好老師背後成功的原因是授之以漁,而非授之以魚。能活血造血的開源工具無疑可以成為我們的導師,所以強調跨域資料分析的資料科學家們一定要先慎選工具,再去投注寶貴的時間和金錢在艱辛的學習過程中。今年 2017 R Taiwan 研討會從雲端、產業根基與解析技術等開始談起,結合實務經驗豐富的業界專家與尖端研究的學界教授進行分享討論,期使 R軟體發揮槓桿作用,讓各行各業能從彈性的開源工具中取得資料加值的競爭優勢。

會議時間:2017/12/14 星期四 (研討會前導課程)2017/12/15 星期五 (研討會主議程)

會議網址:http://rtaiwan2017.fmhuang.net/
早鳥票報名網址: https://goo.gl/forms/IYUtY76EgRShd2Os1

研討會主議程報名費用:一般票1000元、早鳥票800元
研討會前導課程報名費用:一般票1500元、早鳥票1200元
線上報名時間:
早鳥票:2017/11/01 中午 12 點開始至 2017/11/13 中午 12 點止
一般票:2017/11/14 中午 12 點開始至 2017/12/04 中午 12 點止

主辦單位:東吳大學巨量資料管理學院、中華R軟體學會、臺灣資料科學與商業應用協會、國立臺北商業大學
合辦單位:國立臺北商業大學資訊與決策科學研究所
協辦單位:中華R軟體學會、臺灣資料科學與商業應用協會、臺北市商用數學發展協會、國立臺北商業大學,資訊與決策科學研究所、國立臺北商業大學,資料科學應用研究中心、臺北醫學大學,管理學院、臺北醫學大學,大數據研究中心、中華市場研究協會、中華資料採礦協會、國立台北科技大學,資訊與財金管理系、東吳大學,財務工程與精算數學系。

# end

2017年10月19日 星期四

by + scale 資料物件轉換為資料框並與原始資料合併




關鍵字:

  • scale 資料標準化
  • by 依群組計算
  • rbind 列合併
  • do.call 執行R函數
  • cbind 行合併

分析:

  • by 函數:提供不同群組資料執行計算. by 的結果為 list, 使用 as.data.frame 會有錯誤, 此時改用 data.frame(do.call("rbind" , x)) 即可解決此問題. 
  • scale 函數: 將資料值進行標準化轉換, (x - u)/s, u:平均值, s:標準差.
  • rbind 函數是上/下資料的列結合.
  • cbind 函數是左右資料的行結合.


R程式解說:

  • [#1] 先將 by的結果儲存成資料物件x.
  • [#2] 使用 data.frame, do.call, rbind 將by結果合併.
  • [#3] 使用 cbind 將by 結果與原有資料物件進行行合併.
  • [#4] head 函數預設顯示前6筆資料.

R程式:

x <- by(iris[-5], iris$Species, scale) # list
tmp <- data.frame(do.call("rbind", x))
iris.scale <- cbind(iris, tmp)
head(iris.scale)
# end



2017年10月8日 星期日

蒙地卡羅法估計圓周率 (Monte Carlo Estimation of pi)

# 蒙地卡羅法
# 圓周率
# 模擬
# Monte Carlo
# pi
# Simulation
# Programming
# plotrix package

分析:

  • 本範例使用蒙地卡羅法以估計圓周率 pi。
  • 考慮正方形邊長為單位1, 面積為1平方單位. 中間包括圓形, 其半徑為0.5單位,圓形面積 = pi*(0.5^2) = pi/4.
  • 在正方形內隨機產生一個點(x, y), 點落在圓型內的機率大約等於扇型面積 / 正方形面積 = pi/4。
  • 考慮(x,y)落在圓形中的次數有numberInCircle次, 全部實驗次數為N, 則 pi 大約等於:
    (numberInCircle / N ) * 4.


R程式解說:

  • [#3] 首先載入 plotrix 套件, 準備繪製圓形圖.
  • [#4] 設定亂數種子, 準備隨機產生資料點(x, y)
  • [#5] mfrow=c(1,2) 設定繪圖結果為1列2行
  • [#6] 設定 type="n" 表示不繪圖. 因後續須客製化座標軸, axes=FALSE 不繪製座標軸, asp=1設定長寬比例為1. main 標題使用 expression 以標示數學符號pi.
  • [#7] 使用 draw.circle {plotrix} 繪製圓形圖.
  • [#8] 使用 lines 繪製正外形.
  • [#9] 加上4個角落位置的座標軸.
  • [#11] 建立 pi.simulation 函數
  • [#16-18] 如果點(x, y) 與 (0.5, 0.5) 之距離小於或等於0.5, 表示此點位於圓形之內, 此時 numberInCircle 個數加上1, 繪製紅色點.
  • [#21] 點(x, y)不位於圓形內, 繪製藍色點.
  • [#23] 計算估計pi值並儲存於 c 物件.
  • [#25] 回傳估計pi值物件c.
  • [#30] 取出最後一次計算之估計pi值.
  • [#39] 使用 par 函數將繪圖結果還原成1列,1行.

蒙地卡羅法估計圓周率圖:




R程式:



# title: Monte Carlo Estiamtion of pi
# date: 2017.10.9
library(plotrix)
set.seed(123)
op <- par(mfrow=c(1,2))
plot(0.5, 0.5, xlim=c(0, 1), ylim=c(0,1), type="n", axes=FALSE, asp=1, xlab="", ylab="", main=expression(paste("Monte Carlo for ", pi)))
draw.circle(0.5, 0.5, radius=0.5)
lines(x=c(0,1,1,0,0), y=c(0,0,1,1,0))
text(c(0.05,0.95,0.95,0.05), c(0.05,0.05,0.95,0.95), c("(0,0)", "(1,0)", "(1,1)", "(0,1)"))

pi.simulation <- function(samplesize) {
 c <- rep(0, samplesize)
 numberInCircle <- 0
 for (i in 1:samplesize) {
 x <- runif(2)
 if (sqrt((x[1]-0.5)^2 + (x[2]-0.5)^2) <= 0.5) {
 numberInCircle <- numberInCircle + 1
 points(x[1], x[2], col="red", pch=".")
 }
 else {
 points(x[1], x[2], col="blue", pch=".")
 }
 c[i] <- (numberInCircle / i) * 4
 }
 return(c)
}

size <- 10000
pi.sim <- pi.simulation(size)
estimation.pi <- pi.sim[size]
text(0.5, 0.5, paste0("sample size=", size), cex=1.5)
text(0.5, 0.4, paste("pi=",estimation.pi))
plot(pi.sim, type="l", 
 main="Monte Carlo method for pi",
 xlab="samples", ylab=expression(paste("Estimation of ", pi)), 
 ylim=c(0,6))
abline(h=estimation.pi, col="red", lty=3)
grid()
par(op)
# end

2017年9月22日 星期五

R等差數列(前後二個元素相減)與條件式計算


主題: 如何計算數列的前後二個元素相減產生的新等差數列與不同條件式數值計算

# diff
# c
# lapply
# <<

分析:
  • 感謝R友 Bic Ton提供此問題
  • 考慮 x <- c(2,5,6,1,3,8,4,5,6), 計算數列中後面減前面之結果, 例: 5-2=3, 6-5=1,...., 此時可以使用 diff(x).  diff 函數預設會計算數列中第2個元素減第1個元素的結果. 參閱線上說明 ?diff 
  • 如果希望計算前面減後面之結果, 例: 2-5=33, 5-6=-1,...., 此時可以使用 -diff(x)
  • 本例考慮3種不同解法:
    (1). 使用 for 迴圈判斷等差數列每個元素為正數或負數, 再加總結果.
    (2). 使用 for 迴圈判斷等差數列每個元素為正數或負數, 將結果用 c() 向量組合成新向量, 此方法可儲存正數與負數之結果. 注意: 本方法因為要儲存新等差數列的結果, 因此執行時間須較久.
    (3). 使用 lapply 函數以進行不同條件式數值計算, 本例使用全域變數指派符號 << 以利回傳計算結果.
  • 本例第2種方法 c() 分別儲存正/負數結果, 如果資料量較大時, 例: 100萬筆資料以上, 其執行時間須很久, 此時須考量其他儲存方式, 例: 使用 list串列 以加速程式之進行.
R程式碼:

x <- c(2,5,6,1,3,8,4,5,6)
x.diff <- -diff(x)
x.diff
# method 1 using for loop
gain1 <- 0
loss1 <- 0
for (i in 1:length(x.diff)) {
 if (x.diff[i] >= 0) gain1 <- gain1 + x.diff[i]
 else if (x.diff[i] < 0) loss1 <- loss1 + x.diff[i]
}
gain1
loss1

# method 2 using for loop with c()
gain2 <- c()
loss2 <- c()
for (i in 1:length(x.diff)) {
 if (x.diff[i] >= 0) gain2 <- c(gain2, x.diff[i])
 else if (x.diff[i] < 0) loss2 <- c(loss2, x.diff[i])
}
sum(gain2)
sum(loss2)

# method 3
gain3 <- 0
loss3 <- 0
lapply(x.diff, function(x) {
 ifelse(x > 0, gain3 <<- gain3 + x, loss3 <<- loss3 + x)
 #return(c(gain, loss))
 })
gain3
loss3
# end

2017年9月10日 星期日

R讀取中文檔案產生亂碼等錯誤問題

主題: R讀取中文檔案產生亂碼等錯誤問題

說明:

# read.table
# encoding="UTF-8-BOM"
# ANSI
  • 感謝R友-阿賢提供 encoding="UTF-8-BOM"解決亂碼問題.
  • 使用R讀取文字檔時, 有時會遇到資料匯入有錯誤訊息或中文亂碼問題.
  • 資料來源: https://data.gov.tw/dataset/35131, 匯入 open data 空氣品質監測小時值(一般污染物,每日更新) 所產生的問題與解決方式.
  • 匯入資料 read.table {utils} 常用參數:
    (1). fill = TRUE --> 使用時機: 錯誤訊息為 line x did not have xxx elements.
    (2). encoding --> 結果為亂碼.
    (3). fileEncoding  --> 結果為亂碼.
  • 考慮 Windows 執行環境, 如果有亂碼問題, 最簡單的解決方式之一是使用記事本開啟檔案, 另存新檔 畫面中, 編碼改為 ANSI.


  • 使用 Notepad++ [https://notepad-plus-plus.org/zh/] 開啟檔案, 在視窗右下角狀態列會有"UTF-8-BOM"編碼, 此時可加上 encoding 或 fileEncoding 參數, 本例使用 fileEncoding = "UTF-8-BOM" 參數即可完成匯入資料.


  • 另可參考: 資料集為CSV檔,打開來為亂碼,怎麼辦? https://data.gov.tw/node/18765
  • 結論: 使用另存ANSI編碼或加入fileEncoding="UTF-8-BOM"參數應該可以解決亂碼問題.
# 執行畫面:



# R程式碼:
x1 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",") # line 1 did not have 31 elements
x2 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",", encoding="UTF-8") # line 1 did not have 31 elements
x3 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",", encoding="UTF-8", fill=TRUE) # 欄位錯誤
x3[1:3, 1:6]
x4 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",", encoding="UTF-8-BOM") # line 1 did not have 31 elements
x5 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",", encoding="UTF-8-BOM", fill=TRUE) # 亂碼
x5[1:3, 1:6]
x6 <- read.table("ATM00626_20170910170405.csv", header=TRUE, sep=",", fileEncoding = "UTF-8-BOM", fill=TRUE) # OK
x6[1:3, 1:6]
x7 <- read.table("ATM00626_20170829185638-ansi.csv", header=TRUE, sep=",") # OK
x7[1:3, 1:6]
x8 <- read.csv("ATM00626_20170910170405.csv", header=TRUE) # 亂碼
x8[1:3, 1:6]
x9 <- read.csv("ATM00626_20170910170405.csv", header=TRUE, fileEncoding = "UTF-8") # invalid input found
x10 <- read.csv("ATM00626_20170910170405.csv", header=TRUE, encoding = "UTF-8-BOM") # 亂碼
x10[1:3, 1:6]
x11 <- read.csv("ATM00626_20170910170405.csv", header=TRUE, fileEncoding = "UTF-8-BOM") # OK
x11[1:3, 1:6]
# end

2017年9月9日 星期六

網路抓取 R CRAN 套件清單, 使用 ggplot2 套件繪圖, 建立第2個y軸座標.

主題: 網路抓取 R CRAN 套件清單, 使用 ggplot2 套件繪圖, 建立第2個y軸座標.
說明:

# ggplot2
# packages list
# XML
# geom_col
# goem_line
# geom_point
# scale_y_continuous



  • [#1-2] 首先載入 XML, ggplot2 套件.
  • [#3] 使用CRAN網站-依日期排列抓取現有1萬多個套件清單,
    例: http://cran.csie.ntu.edu.tw/web/packages/available_packages_by_date.html
  • [#4] 使用 readHTMLTable {XML} 函數以讀取網站中的套件清單表格, 將結果儲存為mydf資料物件.
  • [#7] 使用 trimws {base} 函數以刪除欄位名稱空白字元.
  • [#8] 原匯入第1欄 Date為字串資料型態, 使用 as.Date {base} 轉換為 日期(Date) 資料型態.
  • [#9] 使用 format {base}並取出套件更新年, 使用 table {base} 以計算各年套件個數.
  • [#10] 使用 cumsum {base}計算累計套件數並新增為 AccumulatedPAckages 欄位.
  • [#11] 使用 names {base} 設定前二欄名稱為 Year, Packages.
  • [#14] 使用 geom_col {ggplot2} 繪製"套件數(年)"長條圖, 另可使用 geom_bar {ggplot2} 繪製.
  • [#15-18] 設定主標題, x軸標題, y軸標題. theme {ggplot2} 可設定標題左右置中.
  • [#19] 使用 scale_y_continuous {ggplot2}可在繪圖區之右側建立y軸第2座標軸. 右側y軸對應長條圖的刻度. ggplot2 採用資料轉換概念, 因此左側y軸第1座標軸的刻度,對應至累計套件數, 其中最大值約11405, 右側最大值約4502, 11405/4502=2.5, 考慮以2倍計算, 即將左側刻度除以2, 轉換為右側刻度, 一般使用 trans = ~. /2 表示.
  • [#21] 使用 annotate {ggplot2} 可加上文字標題.
  • [#22-24] 使用 goem_line {ggplot2} 繪製累計套件數線圖.
  • [#25-27] 使用 goem_point {ggplot2} 繪製累計套件數點圖.

library(XML)
library(ggplot2)
urls <- "http://cran.csie.ntu.edu.tw/web/packages/available_packages_by_date.html"
mydf <- readHTMLTable(urls, stringsAsFactors = FALSE)
str(mydf)
mydf <- mydf[[1]]
names(mydf) <- trimws(names(mydf))
mydf$Date <- as.Date(mydf$Date)
mydf <- data.frame(table(format(mydf$Date, "%Y")), stringsAsFactors=FALSE)
mydf$AccumulatedPAckages <- cumsum(mydf$Freq)
names(mydf)[1:2] <- c("Year", "Packages")
mydf
ggplot(mydf, aes(x = Year)) +
geom_col(aes(y = Packages*2), fill="grey") + 
theme(plot.title = element_text(hjust = 0.5)) + 
xlab("年 (2005.1~2017.9)") + 
ylab("累計套件數") +
ggtitle("CRAN-套件統計圖 - by RWEPA") + 
scale_y_continuous(sec.axis = sec_axis(trans = ~. /2)) + 
theme(axis.title.y = element_text(colour = "blue")) + 
annotate("text", x=13, y=5500, label= "套件數(年)", angle = -90) + 
geom_line(aes(y = AccumulatedPAckages, group = 1), linetype = 2, color = "blue") +
theme(axis.text.x = element_text(face="bold", angle=45),
axis.text.y = element_text(face="bold", color="black", size=10)) +
geom_point(aes(y = AccumulatedPAckages, group = 1), color = "blue") +
geom_text(aes(y = AccumulatedPAckages, label = AccumulatedPAckages), 
vjust = -0.3, size = 4, color = "blue")
# end


2017年8月5日 星期六

圖例中顯示點線混合符號 legend

主題: 使用R繪製圖形時, 圖例採用點線混合同時繪圖

說明:

在 legend 函數中, pch 或 lty  等參數中設定為NA 即可混合使用點線等圖例.



# title: 圖例中顯示點線混合符號
# date: 2017.8.5
x1 <- c(1,2,5,4,3)
x2 <- c(4, 1.5, 1.9, 3, 7)
x3 <- c(5, 6, 9, 3, 1)
ymax <- max(c(x1,x2,x3)) + 1.5
plot(x1, type = "b", pch = 19, lty = 1, col = 1, 
 ylim = c(0, ymax), 
 main="圖例中顯示點線混合符號")
points(x2, pch = 17, col = 2)
lines(x3, lty = 2, col = 4)
legend("topleft", legend = c("x1", "x2", "x3"),
 pch = c(19, 17, NA), lty = c(1, NA, 2),
 col = c(1, 2, 4), text.col = c(1,2,4))
# end