R 學習日誌: 2008

2008年11月8日星期六

時間數列 (二) – R 實作

R提供一個特殊的資料型態 ts, 作為Time Series基本類別，並提供相關方法, 作法與 matrix 類似, 這真是太厲害了，完全是OOP的標準作法，讓我對R的發明者更欽佩三分。

以下就是一個簡的例子
> library(MASS) # 匯入測試資料
> plot(huron,type="l") # 製作線圖
> setwd("j:/r/my") # 設定工作目錄
在繪圖視窗按滑鼠右鍵, 選【儲存BMP】，再貼到【小畫家】中存檔
> y2=huron$mean.height*2 # 產生應變數數列
> y=ts(y2) # 將data frame資料轉為time series資料型態
> x=huron$mean.height # 產生自變數數列
> x=ts(x)
> x2=lag(x,-1) # 產生第二個自變數數列
> d=ts.union(y,x,x2) # 組成一個ts資料來源
> lm(y~x+x2, data=d) #產生線性模型
lm(formula = y ~ x + x2, data = d)
Coefficients:
(Intercept) x x2
2.572e-13 2.000e+00 -4.273e-17

時間數列 (一) – 基本理論

一般迴歸分析(Regression)是假設資料誤差是完全隨機的，但現實社會中的現象，卻往往不是如此，例如股價，它的波動會受到前一日或更前幾天的影響，因此，要建構這類的模型，我們就可以應用時間數列(Time Series)，來分析或預測未來的走勢。

在談到時間數列前，我們最好先了解迴歸模型，以簡單線性迴歸為例:

y=a+bx+ε

其中ε為誤差, 假設他服從常態分數，平均數為0，變異數為σ²，且ε_i及ε_j無關聯

之後我們用誤差的最小平方法(Least Square)，求出a及b，進而得到

總平方和(SST) = 殘差平方和(SSE) + 迴歸平方和(SSR)
R² = SSR / SST，R²是檢定x與y的相關係數
檢定一元迴歸可用自由度=n-2的t分配，多元迴歸則用F分配
F = (SSR / (k -1)) / (SSE / (n – k))
其中k-1為自變數個數

時間數列與迴歸分析不同點就是假設誤差不是完全隨機的, 而是自我相關的，即前期的誤差會影響後期的誤差

ε_t= ρε_{t -1} + a_t

其中ρ絕對值<1，a_t稱為干擾項(Disturbance Term)，服從常態分數，平均數為0，變異數為σ²

檢定ρ是否為0，一階(誤差只與前一期有關)可用Dubin-Watson統計量檢定，多階可用ARIMA統計量。
ARIMA： Auto Regression Integrated Moving Average，為Box/Jenkins提出

此外，若考慮不是σ²固定常數，也就是異質變異數(Heteroscendasticity)，可用加權最小平方法處理，這就是ARCH或GARCH模型。

註:此文主要是節錄自林茂文先生所注之時間數列之分析與預測一書。

2008年11月1日星期六

迴歸(Regression)

重拾統計, 真是令人緬懷當年學生生活.

就從迴歸作起吧.

載入測試資料, 也可以自己輸入
> load("j:/r/programs/testdata/usingr.rdata")
使用 moths 變數, 先輸入moths, 查看資料
以moths$meters 為應變數, moths$A、moths$P 為自變數
> lm(moths$meters ~ moths$A+moths$P)
得到簡單迴歸的截距(Interception)及斜率(Slope)
Coefficients:(Intercept) moths$A moths$P
45.1017 0.3712 -0.6005

R 的編輯器

R提供較簡易的介面, 方便使用者輸入資料, 試試看下列方式吧.

b=de(a) # de 為 data.entry的簡稱, b 資料型態為 list
b=fix(a) # fix 為指令編輯視窗
b=edit(a) # 可輸入數學式子, 也可呼叫外部編輯器

執行外部程式

R可以下指令執行外部程式, 執行DOS指令時須加command /c, 用法如2

system("c:/windows/notepad.exe")
system(paste(Sys.getenv("COMSPEC"),"/c", "dir"))

相反的, 外部程式要執行R的程式也有兩種

R xxx.r # xxx.r 為R的程式碼
使用com介面, 須加裝 RSrv250_pl1.exe, 內附範例可參考

R 變數類型

R的變數類型有下列幾種:

向量: 為預設類型, 即一維陣列
Data frame: 具有欄位名稱的二維陣列
List: 可儲存不同類型的物件集合, 例如summary(), 要轉成向量, 可用unlist()
S3: 可具有class的list, 如果list被函數呼叫時, R會使用class的方法取代標準的方法, 請看下例:
# 取三個常態分配隨機亂數
h=list(a=rnorm(3), b= "it cannot be printed ")
# 設定h變數的class名稱為MyClass
class(h)= "MyClass "
# 實作MyClass函數
print.MyClass=function(x){cat("A = ", x$a, "n ")}
# 列印h, 發現R以MyClass取代print
print(h)
S4: 類似S3, 有點複雜, 以後再說明
日期變數: 要使用日期函數處理, 必須先把變數轉為日期
l 如為字串變數且格式為yyyy-mm-dd, 可用as.Date轉換
l 如為數字變數且格式為yyyymmdd, 須先用as.character轉為文字再使用as.Date轉換
> a=20081028
> as.Date(as.character(a), format="%Y%m%d")
[1] "2008-10-28"
> as.Date(as.character(a), format="%Y%m%d")+1
[1] "2008-10-29"

日期格式設定可輸入help(“strtime”)查閱

檔案處理

自外部讀入一個文字檔非常容易, 只要一個指令

mydata=read.table("J:\\R\\Programs\\testdata.dat",header=TRUE, sep=",")

上例第一個參數為文字檔名, 第二個參數指明文字檔第一列為欄位名稱, 第三個參數指出欄位間分隔符號為逗點(,),其他可設定的參數非常多,可輸入指令help("read.table")查閱。

自R輸出一個文字檔只要改read為write即可

write.table(mydata,"J:\\R\\Programs\\testdata2.dat", sep=",",row.names = FALSE)

其中row.names = FALSE表示不輸出列號, 否則, R預設會輸出列號。

大檔案的處理

大檔一次讀入R, 會造成記憶體不足, 這時必須採取串流的方式, 一行行或一段段的讀, 才不會出問題, 指令如下

> con=file("j:/r/programs/testdata.dat","r") #開啟檔案

> readLines(con,1) #1表讀一列, 可一次讀取100列

> close(con) #關閉檔案

R 常用指令

變數指定
> a=1
變數預設資料型態為向量, 亦稱為集合(collection), 故使用c
> x = c(10.4, 5.6, 3.1, 6.4, 21.7)
特殊變數指定方式
> x = 1:12 # x=1,2,3,4, ..., 12
> a = seq(-5,5,by=.2) # -5~5間隔0.2產生一個值
計算結果
> a * 2
列出所有變數
> ls()
列出所有變數及資料型態(data type)
> ls.str()
註解(Comment): # 開頭
幫變數加註解
> comment(x1) = c("This data is from experiment #0234", "Jun 5, 1998")
移除變數 c, d
> rm(c,d)
移除所有變數
> rm(list=ls())
查看指令說明
> help(“rm”)
搜尋相關說明
> help.search("Time Series")
查看套件說明
> library(help="stats")
二維資料型態 data frame, 它具有欄位名稱, 如下例, 欄位名稱分別為f1、f2, 資料有三列
> a = data.frame(f1=c(1,2,3),f2=c(2,4,6))
取其中一欄的值
> a$f1
兩個向量合併為一data frame
> x=cbind(a,b) # 若 a, b 元素個數不等, R會自動重複元素, 已補齊至最大元素個數
> x=rbind(a,b) # a, b 元素已列方式合併
以特定欄位比對合併兩data frame為一data frame
#以兩data frame 的日期比對, 合併
> x2=merge(mydata,mydata2, by.x="Date", by.y="Date")
判斷資料型態是否為矩陣
> is.matrix(b)
顯示資料型態
> class(b)
矩陣轉置(transpose)
> t(b)
矩陣相乘, 須用%*%, 若用*, 系統會將相同位置的兩數相乘, 而非矩陣相乘
> t(b)%*%b
引入一段R程式, 假設你寫好一段公用函數, 你可以以下指令執行, 而不必每次重打, 或複製/貼上方式執行
source("c:/1.r")

R 的環境簡介

個人認為 R 的強大功能有兩項:

專為數學/統計設計的程式語言, 變數不必宣告, 直接使用, 函數庫(Library)不必指定, 只需載入即可, 向量、矩陣、微積分、統計分配函數及轉換樣樣不缺, 繪圖能力超強。
套件(Package)特多, 官方網站數量就超過百個, 也可以自製套件。

當然, 也不是沒有缺點, 套件說明稍嫌簡略, 不過這是Open Source的普遍現象吧.

在程式撰寫前, 應先了解開發環境:

不知道怎麼開始, 輸入help.start(), 即可連到各種使用手冊的首頁。
設定現行目錄: 執行選單【檔案】>【變更現行目錄】或直接下指令如下, 設定程式存放的目錄, 目錄分隔符號(\)須改為斜線(/)或使用雙倒斜線(\\)。
setwd("J:\\R\\Programs")
載入套件: 載入新套件前, 先下指令列出目前已載入的套件
library()
輸入下列指令載入套件
install.packages("graphics")
之後選擇下載的地點即可開始下載，並同時載入至目前環境中, 以上指令也可以執行選單, 達成相同效果。
儲存工作空間(Work Space): 輸入下列指令或執行選單【檔案】>【儲存工作空間】, 將之前的設定儲存到一個組態檔(*.rdata)中。
save.image("J:\\R\\Programs\\My.rdata")
載入工作空間(Work Space): 下次執行R時, 可先輸入下列指令或執行選單【檔案】>【載入工作空間】, 載入上次的設定，

R的通則:

指令大小寫有區別
函數一律要加()

R 的初次邂逅

初次看到R, 就被它的語言彈性及眾多的套件(Package)深深吸引了, 這麼好的軟體, 竟然是Open Source, 而且它在1980年代就被發明了, 可說是阿公級的產品, 歷久不衰, 在軟體技術不斷翻新的浪潮中, 更顯彌足珍貴, 只能說不用可惜啊 !

以下就簡單說明R的安裝與使用:

安裝: R的安裝很簡單, 可以在Windows或Linux安裝, 至 http://cran.r-project.org/ 下載並執行, 就搞定了。

啟動: 從Windows選單, 執行【開始】>【程式集】>【R】>【R 2.7.0】。

使用: 出現畫面如下, 在提示符號後輸入指令即可。例如
a=1
b=2
a+b

離開: 輸入q(), 按ENTER即可。

2008年11月8日星期六

時間數列 (二) – R 實作

時間數列 (一) – 基本理論

2008年11月1日星期六

迴歸(Regression)

R 的編輯器

執行外部程式

R 變數類型

檔案處理

R 常用指令

R 的環境簡介

R 的初次邂逅

R 學習日誌

關於我自己

網誌存檔

2008年11月8日 星期六

2008年11月1日 星期六

關於我自己

網誌存檔

2008年11月8日星期六

2008年11月1日星期六