Rの使い方
基本
- 変数名は英数字、ドット . 、先頭に数字は不可、大文字小文字は別
- アンダーバー _ は使えない(昔代入記号として使われていた)。今は使えなくもないが、混乱を避けるために使わないほうがいいかも。
- コメント: 先頭#,または
1if(0){}で囲う
- 代入
1<-
2x<- 6
3y<- x+2- べき乗
1^ あるいは **- 割り算と余り
1 %/% 整数除算
2 %% 剰余- 演算の優先順位
:の優先順位が一番高いので、
1:3-1と1:(3-1)は違うよ!!
1 a[1:4-1]とa[1:(4-1)]の違いに注意!!- まいなす記号
< <= :: -記号をつけるときは、代入と間違われないように、空白をきちんとつける
- 論理演算
1!
2& && 論理積
3| || 論理和
4xor()
5TRUE, FALSE (T,F)- 分岐と繰り返し
1if(条件式){式1}else{式2}
2switch(条件式,ケース)
3repeat{式}
4while(条件){式}
5for(i in M){式}
6next :
7break- 関数
1max(),min()
2abs
3sqrt()
4round(),floor(),ceiling() 丸め、切り捨て、切り上げ
5sum(), prod() ベクトルの要素の総和、要素の総積
6log(), log10(), log2()
7exp()
8sin(),cos(),tan()
9asin(),acos(),atan()
10pi
11Inf,-Inf
12
13mode()
14attributes(),attr(), structure() 属性
15typeof()
16
17ls()で登録されている名前一覧表示
18
19rm(hoge)で名前hoge削除help
1help(hoge)初期化
1rm(list=ls(all=TRUE)) #初期化
2setwd("/media/takiyama/m2sata/docs/Research/AUT_ENG/imep_ana")型
- character:文字列
- numeric:実数
- integer:整数
- logical:論理値
データ構造
ベクトル:c()
関数c(): combine,あるいはconcatenateの略
1x<-c(4,1,5)とすると処理だけ
2(x<-c(4,1,5))とかっこで囲むと表示までしてくれる中身が空のベクトル
1#直接大きさ指定はできないので以下のようにする
2a<-c()
3a[1:100]<-NA差の変化が規則正しいベクトル
1seq(2, -5, by=-0.5)で2,1.5,1,...-4.5,-5が得られる
2
3#rep()である一定パターンの繰り返し
4rep(c(3,2,1),3)で、3,2,1,3,2,1,3,2,1が得られる。行列:matrix
二次元 :defaultは列方向に値が設定される
1matrix(c(要素),nrow=nr,byrow=TRUE)#こうしたら行(row)方向に値が設定中身が空の行列:
1matrix()::中身,大きさが空
2matrix(nrow=2,ncol=3),行数(nrow)列数(ncol)指定
3matrix(要素,nrow=2,ncol=3)大きさ指定したいならば要素必要
1mm<-matrix(0,2,3)#なら空行列できる
2mm<-matrix(NA,2,3)#なら空行列できる
3行列zをベクトルとすると
4dim(z)<-c(2,3)で2行3列の行列要素を行方向に設定
1a<-c(1,2,3,4,5,6)
2matrix(a, ncol=3)で(2行)3列のmatrix
3matrix(a,nrow=3)で3行(2列)のmatrixができる。
4
5matrix(a,nrow=2,byrow=T) :: byrow=Tで左から右へ行方向に要素が埋められ、
6指定がない時は、上から下へ列方向に要素が設定
7
8matrix(data, nrow,ncol,byrow=FALSE)
9
10t(a)で転置行列単位行列は
1diag(4)配列
1array()多次元,arrayは行列の拡張版
m行,n列の行列を重ねたもの
1array(1:20,c(4,5))#で4行5列の行列、列から1,2..と埋めていく
2
31 5 ... 17
42 6 ... 18
53 7 ... 19
64 8 ... 20となるデータができる。
以上の構造中の要素はすべて同じ型。 数値、文字などを混在できない。 自動で変換されてしまう。
配列の大きさを調べる
1length(hoge)#で長さ
2dim(hoge)#で行列の大きさ配列のアクセス
- 要素の抽出
1a[2,3]
2a[2,]---二行目
3a[,3]---三列目- 1:10 で1から10まで。
:の優先順位が一番高いので、 1:3-1と1:(3-1)は違うよ!!
data.frame()
列ごとに異なる型を持たせることができる
https://htsuda.net/stats/data-basics.html
つまり1行方向はあるものの特性
1id name age sex married
21 Alice 10 Female FALSE
32 Bob 15 Male FALSE
43 Charlotte 20 Female FALSE
54 Dick 35 Male TRUEtable(x)
配列に含まれる値の出現数をカウント
https://www.marketechlabo.com/r-time-series-analysis/
ts
ts: Rの基本の時系列オブジェクト。ほとんどの時系列ライブラリはこの形式の時系列データを扱う
xts
xts: データフレームとtsの間に位置づけられる、時系列データを便利に扱えるようにした形式
zoo
zoo: データフレームとxtsの間の形式
少し高度な型
- factor(因子型)
- factorは因子型:カテゴリカルなデータを扱うときに使う
https://htsuda.net/stats/data-basics.html
- factor()はベクトルを因子にする。
- gl()は水準パターンを与えて因子を作る
http://www.okadajp.org/RWiki/?因子Tips大全
データ変換
一次元データを二次元に変換
1#dをnrow行,byrow=Tで左から右へ行方向に要素が埋められる
2#d(col*DATA_NUM)をs(DATA.NUM行,col列)に変換
3s<-matrix(d,nrow=DATA.NUM,byrow=T)
4s<-matrix(d,ncol=col,byrow=T)#列方向指定しても同じ型とデータ構造の変換
型変換:as.*****()
1as.numeric()
2as.character()
3as.vector()
4as.matrix()
5as.data.frame()問い合わせ
1is.*****()問い合わせ::ifなどの判断で用いる
(=, ! < >などの比較演算子以外に利用可)
データ型の確認
1class(hoge)配列の加工と操作::データ構造の演算
大きさを指定した配列の変数に代入という考え方はできない。
すべて、ベクトルの結合(リスト)、という考え方。
結合後に操作する。
コマンド
ベクトルx,yの結合
1c(x,y)#xの後ろにyを結合
2append(x,y)#xの後ろにyを結合
3append(x,y,after=2)#xの2番目以降にyを挿入置換
1replace(x,y,z)
2#x:置換対象
3#y:置換位置
4#z:置換内容1x<-c(1,2,3,4,5,6)
2y<-c(2,4)
3z<-c(10,20)
4w<-replace(x,y,z)
5w:(1,10,3,20,5,6)sort(x)昇順で並べ替え
1sort(x,decreasing)#降順で並べ替えorderで,大きさ順がわかる
1s<-c(1,10,2)
2order(s)
31 3 2
4rev(order(s))ベクトルの結合
- rbind:row(行)要素としてbind
- cbind:col(列)要素としてbind
1> a<-c(1,2,3)
2> b<-c(5,6,7)
3> d<-cbind(a,b)
4> d
5 a b
6[1,] 1 5
7[2,] 2 6
8[3,] 3 7
9> e<-rbind(a,b)
10> e
11 [,1] [,2] [,3]
12a 1 2 3
13b 5 6 7- data.frame(元,追加):行データを追加
- list:大きさが異なっても結合できる
クラスを保持したまま結合するのはlistがおすすめ by: https://www.karada-good.net/analyticsr/r-317/
データフレームに行追加をしないほうがいい::型が変わってしまう https://www.jaysong.net/RBook/datastructure.html
ベクトルに名前つける
1names(x)<-c("a","v")
2names(x)=NULLデータフレームの列方向毎に名前をつける
1FU<-1
2IMEP<-2
3
4colnames(sd)<-c(
5"FU",
6"IMEP")#列変数名1列目にFU,2列目にIMEPと名前がついた
名前を用いた演算
変数名$データ名で取り出せる
1IMEP.mean<-mean(sd$IMEP)
2ymax<-max(sd$IMEP)
3
4j<-1
5#IMEP.fil<-array(0,dim=c(DATA.NUM))
6IMEP.fil<-c()#大きさを指定しない
7for (i in 1:DATA.NUM){
8 if(sd$IMEP[i] < IMEP.mean+150){
9 IMEP.fil[j] <- sd$IMEP[i]#ベクトルに追加
10 j<-j+1
11 }
12}ハイテクニックな結合
which(names(df2)=="c")
で"c"というラベルの列を参照できる。
https://rstudio-pubs-static.s3.amazonaws.com/573963_94a58958207d44d884ab52d2f8a8ffac.html
問題:a,b,eを縦方向に結合する
1df1 <- data.frame(a = 1:3, b = c("a", "b", "c"), e = c(1, 3, 5))
2df2 <- data.frame(a = 7:9, b = c("t", "t", "t"), c = 4:6, d = c("q","w", "e"), e = c(6, 6, 6))
3
4df2a <- c(which(names(df2) == "c"), which(names(df2) == "d"))
5df2b <- df2[ , -df2a]
6rbind(df1, df2b)- df2からcとdを削除::
- df2aにcとdを集めて、それをdf2から引く->df2bとする
- その後df1とdf2bを結合
データ読込
1 read.table("ファイル名",header=ヘッダ行有無,sep="区切り文字")
2
3 read.csv("hoge") はread.table("hoge",header=T,sep=",")
4 read.delim("hoge")はread.table("hoge",sep="\t")
5
6# headerつきがdefault
7# 大規模なデータはscanのほうが効率いいらしい
8
9 scan("hoge",list(No,,l,l,l),skip=1,sep=",")
10
11# skipはヘッダ行を読み飛ばしと、開始行を指定。CSVデータ読み込み
1 d<-read.csv("ファイル名",header=F,skip=")バイナリデータ読み込み
1DATA.NUM <- 2500 #行数
2col <- 9 #列数
3
4con <- file("../data/2023_0712/w_IMEP/AF_const/data_2023_0712_143906","rb")
5d <- readBin(con, what=double(), n=col*DATA.NUM, size=8)
6close (con)
7#配列に放り込むoptionなし,すべて1次元データとして読み込み
8s<-matrix(d,ncol=col,byrow=T)#col列の二次元データに変換
9sd<-as.data.frame(s)#データフレームに変換データ保存
大きさが異なる変数はwrite.{table,csv}で保存できない
1library("dplyr")
2bind_rows()NAで補完してサイズを合わせてくれる
https://a-habakiri.hateblo.jp/entry/2019/11/07/203512
コマンド
1getwd() #現在ディレクトリの設定
2setwd() #WorhDirの設定
3#デリミタは \\ または /
4source("hoge.R")#でソースの読み込みと実行
5dump("hoge.R")#でRのオブジェクトをダンプ
6
7save.image("h")でデータ、関数の保存
8load("h")でデータ、関数の読み込みデータフレーム
data.frame() フレームでいろいろ変えられる
read.tableのようにデータフレームに変換したいときは
1as.data.frame(オブジェクト名)
2
3a<-read.csv("hoge")
4data.entry(a)
5
6#または
7
8fix(a)とすると、読み込んだデータが表形式で現れる。
(データエディタといい、編集もできるそうである)
名前を用いた参照
1 year abc
2Japan 2006, 2
3USA 2000, 3とすると
1a<-read.table("hoge")
2a$yearで列ごと
1a["Japan",]で行ごと 参照できる。
1attach(a)とすると
1year
2abcだけで参照できる
1colnames(a)<- c("hoge","fuga")で追加しても良いし、
1 x1<-c(1,2)
2 x2<-c(2,3)
3 xx<-data.frame(hoge=x1,fuga=x2)とすると名前hoge,fugaの列ができる! -列要素の足し算で、データが構成できる。!!
列要素の加工
1 a1<-2
2 a2<-7
3 a<-a[,-c((a1-1),(a1+1):(a2-1),(a2+1))]とすると aから,1,3:6,8行を取り除く行列の出来上がり!
1a[1:3]#とすると1列から3列
2a[1:3,]#と,をつけないと、行指定の選択にならない
3#:の優先順位が一番高いので、()をつけること!!図を描く
1plot(sd$IMEP,xlim=c(0,DATA.NUM),ylim=c(0,ymax),col='black',ann=F)
2par(new=T)#上書きするおまじない
3#plot(IMEP.fil,type="l",xlim=c(0,DATA.NUM),ylim=c(0,ymax),col='blue')
4plot(IMEP.fil,type="o",xlim=c(0,DATA.NUM),ylim=c(0,ymax),col='blue')
5#type::l:線,o:点と線,p(点,default),b,cなどありpar(…)
グラフィックパラメータの設定
help(par)で参照 http://takenaka-akio.org/doc/r_auto/chapter_05.html
1lheight(line height)
2lwd(line width)
3lty(line type)
4xlog新しい画面を開く
1dev.new()画面を分割する
1layout(matrix(1:2,2,1,byrow=TRUE))#2行1列,行方向に並ぶ
2par(mfcol=c(2,1))layout(hoge)の引数は分割サイズと図の配置位置も指定できる。
1 #1行目に図を1つ、2行目2列に分けて、右側に図を書く
2 nf <- layout(matrix(c(1,1,0,2), 2, 2, byrow = TRUE), respect = TRUE)
3 layout.show(nf)
4
5 #0は図を書かない、幅3:1,高さ1:3,
6 nf<-layout(matrix(c(2,0,1,3),2,2,byrow=TRUE),c(3,1),c(1,3),TRUE)
7 layout.show(nf)
8 #http://eau.uijin.com/advgraphs/layout.html画面を閉じる
1dev.off()#今見ている画面
2dev.off(2)#指定画面
3graphics.off()#全画面文字を画面に出力する
1cat("IMEP fil mean = " , IMEP.fil.mean ,"\n")
2cat("IMEP fil sd = ", IMEP.fil.sd, "\n")
3
4#printで出力できるのは1つだけ
5#または
6#sprintf("IMEP fil mean = %f " , IMEP.fil.mean): 使えないlmを使う
lmの返すオブジェクト
Rstat.pdf,p132.
1plot(ans$fitted.values,df$AirPassengers)#回帰データと元データをplot
2xlab="x axis name",
3ylab="y axis name",
4xlim=(0, 600), ylim=(0, 600),#描画範囲
5asp=1)
6abline(0,1)#傾き1の直線,fitしてたら一致するはず
7abline(ans)#回帰直線
8
9layout(matrix(1:4, 2, 2, byrow=TRUE))#一画面に複数の図を出す時のおまじない
10par(mar=c(4, 4.5, 1.5, 0.5), mgp=c(2.1, 0.8, 0))
11plot(ans, cex=0.7)
12plot(ans) # 回帰診断グラフp.135(149) 図7.4(p.136) は,以下のようにして描くことができる。
1> plot(x, y, pch=19, cex=1.5) # 元のデータの散布図
2> x1 <- seq(0, 10, by=0.05) # 新しいデータ
3> new.data <- data.frame(x=x1) # データフレームにする
4> y1 <- predict(ans1, new.data) # 新しいデータに対する予測値
5> y3 <- predict(ans3, new.data)
6> y5 <- predict(ans5, new.data)
7> matlines(x1, cbind(y1, y3, y5), # 3 種類の予測値を描画
8+ lty=3:1, col="black")
9> legend("bottomright", # 凡例を付ける
10+ c("1 次式", "3 次式", "5 次式"), lty=3:1)
11==低水準関数:abline
1abline(h=0,v=0)#座標軸を表示
2R-tips p.131(132)
3abline(切片,傾き,h=y,v=x)#ベクトルで与えたy座標の水平線、x座標の垂直線を引く
4abline(result)result は関数lm() で直線回帰を行った結果が入ったオブジェクト で,回帰直線が描かれる.
lm()が返すオブジェクトobj
+ obj<-lm()とすると summary(obj),AIC(obj),coefficient(obj)などで値を取り出すことができる r-tips p.188(190) + obj$coeffciant, obj$residuals, obj$fitted.valuesなどもあり (Rstat, p,123(137)にしか説明はなし -->helpにfited.valuesあり。グラフ表示もできた! ) + R-intro-1.7.0,p.53(59)にはcoeffients(obj)などのみの説明 ($coeffientsの説明はなし) + formula(obj)#モデル式を抽出
coefficients(obj)でも、obj$coefficientsでも同じ結果が表示されるけど、 lmの場合だけできる特典。coefficients(obj)が一般的
多変数回帰
1Y . X1 + X2 + X3 + X4 # モデル式
2
3Df Sum of Sq RSS AIC # 独立変数を除いたモデルのAIC
4- X4 1 3.233 72.763 27.846 # X4 を除くとAIC が小さくなる
5< none > 69.531 29.392 # 何も取り除かないモデルのAIC
6- X1 1 41.797 111.328 32.099
7- X2 1 43.966 113.497 32.292
8- X3 1 165.279 234.810 39.562
9
10Step: AIC=27.85 # このステップでのAIC
11Y . X1 + X2 + X3 # モデル式
12Df Sum of Sq RSS AIC # 独立変数を除いたモデルのAIC
13< none > 72.763 27.846 # 今のモデルのAIC が最小
14- X1 1 39.909 112.673 30.219 # 以下のどれを取り除いても
15- X2 1 84.007 156.770 33.522 # AIC は大きくなる
16- X3 1 168.982 241.746 37.853step で最小のAIC を持つ重回帰モデルlm(formula = Y 〜 X1 + X2 + X3, data
;;Rstat.pdf,p.128(142)
stepの挙動;
- そのバラメータを取り除いたときのAICを計算し、 小さい順に並べる
- 何も取り除かない時 < none > のAICと比較し、それよりもAICが小さい時は そのパラメータを取り除くことで、AICがよくなる(AICが小さくなる)ことを示しているので、そのパラメータを取り除くことにする。
- 一番小さいAICから順に、1つずつ(2)を繰り返す。
- <none>の時のAICが一番小さい値になるまで繰り返す。
- こうして残ったパラメータが有効なパラメータである。これについてPrなどを計算している
時系列解析
Rの標準パッケージに含まれる時系列データの型としてts(time series)がある
https://qiita.com/YM_DSKR/items/2528548913378bfbf9bc
https://logics-of-blue.com/時系列解析_理論編
AIC()
AR
1y(t)=a1*y(t-1)+a2*y(t-2)+...+ap*y(t-p)+epsMA
1y(t)=mu+b1*eps(t-1)+b2*eps(t-2)+...+bq*eps(t-q)+eps(t)ARMA
1y(t)=a1*y(t-1)+a2*y(t-2)+...+ap*y(t-p)+b1*eps(t-1)+...Autoregressive Moving Avarage, yとepsで回帰
ARIMA
Autoregressive integrated Moving Avarage
1y(t)-y(t-d)=c+a1*y(t-1)+a2*y(t-2)+...+ap*y(t-p)データの差分を取ってからARMAを適用したモデルです。
コレログラム
https://qiita.com/nash_efp/items/a35c1796e249277ec43f
acf(IMEP.fil)#自己相関.1,2,3
pacf(IMEP.fil)#偏自己相関:注目している時以外を無視した自己相関係数
1acf
2pacfコレログラムとはそれぞれのラグについての自己相関係数を表した図です。
日ごとのある株の収益率データを使って見ていきましょう。
ラグ0は現在の自分と現在の自分の相関を表しているので当然1.0です。
右にいくほど下がっている、つまり過去にいくほど現在との相関が弱くなって いることがわかります。
varモデル
1library(vars)VARモデル
1今年のサンマ = 去年のサンマ×傾きA + 去年のプランクトン×傾きB +切片C
2今年のプランクトン= 去年のプランクトン×傾きD + 去年のサンマ×傾きE +切片Fhttps://logics-of-blue.com/varモデル
arxと同等???
Rの移動平均は
statsパッケージに含まれるfilterを利用
dplyrとかぶるので、stats::filterとする
1library(tidyverse)
2https://datasciencehenomiti.com/programing/r/dplyr_moving_average/
3mutate(
4stats::filter(data,rep(1,n)/n,sides=1)
5)filterという関数を用いて書く
http://tips-r.blogspot.com/2015/01/r_1.html
1filter(x,c(1,1,1))/n
2c:時点n-1,現時点n,時点n+1,の重み
3c(1,1,1)とrep(1,3)は意味同じ
4rep(1,n)は1をn回繰り返しdplyr
1RcppRoll:roll_sum(x,)
2RcppRoll:roll_mean(x,)zoo::rollapply()もあり https://qiita.com/Sickle-Sword/items/d93aff1d29e1176e25fe
時系列データを扱うパッケージ
zoo dplyr tidyr purrr tibble
などあり https://rstudio-pubs-static.s3.amazonaws.com/640203_28880d5e47b44029b9279d8e24385fba.html
1mutate()データセットに変数を追加する関数
https://kazutan.github.io/kazutanR/hands_on_170730/mutate.html
1mutate(iris,k=S.len*2)irisのS.lenを2倍したkをirisの最後列に追加 名前が同じなら、上書きする 行数が同じでないとだめ
data.frameの仕様と同じ!?!?
パッケージ
package(library)のインストール
tseries,statsはdebianのパッケージになし
suしてRを起動
Rのコンソールから
install.packages("tseries", dependencies=TRUE)
https://stats.biopapyrus.jp/r/basic/package.html
libの指定がないので
/usr/local/lib/R/site-library
にinstall
依存対象::TTR,curl,jsonlite,quadprog,quantmod,xts,zooも入った
library(tseries)
tseriesにadfある
library(gvlma)
1fit1=lm(sa~a+b+c)
2fit2=lm(sa~a+b)
3
4AIC(fit1,fit2)
5#https://htsuda.net/stats/regression.htmllibrary(MASS)
1#stepAIC(fit,direction="backward")
2
3#https://qiita.com/purple_jp/items/b5b36bc78a4ccfb42027
4
5dev.new()
6acf(IMEP.fil)#自己相関.1,2,3
7pacf(IMEP.fil)#偏自己相関:注目している時以外を無視した自己相関係数
8
9imep.ar <- ar(IMEP.fil)#係数表示,31次!!参考文献、URL
forecast
パッケージforecastは、自動SARIMA推定だけでなく、色々な便利な関数が付い ています
https://logics-of-blue.com/時系列解析_理論編
https://logics-of-blue.com/時系列分析_実践編
TSstudio
今回は、「MLmetrics」(精度評価用のライブラリー)と「forecast」(時系 列解析のライブラリー)も利用するので、インストールされていない方は、イ ンストールしておいてください。 https://www.salesanalytics.co.jp/datascience/datascience011/
lbrts.pdf
統計解析ソフトRのスクリプト集
https://www.educa.nagoya-u.ac.jp/~ishii-h/materials/Rscripts_ishii.pdf
R 制御 モデルでgoogleした結果
北川源四郎
Rによる 時系列モデリング入門 岩波
代表的な手法と応用へのポイントを、RパッケージTSSSを用いた分析により実 践的に学ぶ。
統計数理研究所で開発されたRパッケージTSSSの使用法と解析例を新たに多数 追加した。
時系列解析(6) –ARモデルの推定
東京大学 数理・情報教育研究センター http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/時系列解析(6)_2019.pdf
北川源四郎 時系列解析 http://www.mi.u-tokyo.ac.jp/mds-oudan/lecture_document_2019_math7/%E6%99%82%E7%B3%BB%E5%88%97%E8%A7%A3%E6%9E%90%EF%BC%88%EF%BC%96%EF%BC%89_2019.pdf