Rによる機械学習 最尤推定法その1
よくある例題を使って最尤推定法をやってみます
> training [1] 0.02841715 0.34097049 1.37103609 0.52697056 0.88287394 0.20673168 [7] -0.92268387 -1.41716486 -0.94624667 -0.56053292 0.21504500 > test [1] 0.2744391 0.2354888 1.0005627 0.7502364 0.5542613 0.1038660 [7] -0.7504576 -0.2694111 -0.9737203 -0.5597848 0.2300776
グラフはtrainingのデータをプロットしています
ある関数f(x)にノイズεが含まれていますことがわかっています
ノイズεは確率変数で何らかの確率分布に従っています
y=f(x)+ε
このデータからf(x)を推定することを考えます
こういった問題は最小二乗法や最尤推定法が用いられます
今回は最尤推定法を用いてf(x)を推定してみようと思います。
最尤推定法は以下のステップからなります
⓪ノイズの確率分布を仮定する。確率分布を定めるパラメータを導入する(平均、分散など)
①関数形を仮定する。関数を定めるパラメータを導入する(例えばg(x)=axのaがパラメータ)
②観測されたデータが出現する確率Pを計算する
③Pが最大になるようなパラメータを決める
最尤推定法は
「観測されたデータが最も出現確率が高いデータのはずだ」
という仮説のもとでパラメータを決定します
これは例えば表が確率q、裏が確率1-qででるコイン(pは未知)を5回投げて
表、表、表、裏、表
と出た時にq=4/5と推定するのが最も合理的だ、というのと同じ考え方です
ちなみにコインの例で最尤推定法を行うと
観測されたデータが出現する確率PはP=(1-q)q^4となってこれを最大化するようなpは
∂P/∂q=0よりq=4/5となります
TeXかなんかで数式を入れたかったのですが、何が一番良いのか模索中なので続きはまた次回
Rによる統計解析 日経平均先物その2
日経平均先物の分布について確認してみたいと思います
Black-Scholesモデルなど、株価のリターンの分布は対数正規分布と仮定されることが多いですが、実際はどのような分布になっているか調べてみます
今回は2010/6/21〜2016/6/21の日足データを使います
Macであれば楽天証券のMarketSpeedからデータをcsvで保存します
この中で、
> file <- file("NikkeiSakimono.csv", open="r", encoding="cp932") > data <- read.table(file, header=TRUE, sep=",")
ですがこのままだとデータがfactor型になってしまいます
factor型のままでは解析できません
> sapply(data,class) 日付 始値 高値 安値 終値 前日比 出来高 X5DMA "factor" "factor" "factor" "factor" "factor" "factor" "factor" "factor" X25DMA X75DMA X5DVMA X25DVMA X75DVMA "factor" "factor" "factor" "factor" "factor"
桁区切りを無くせばread.tableした時にnumeric型になってくれます
ということでExcelで桁区切りを一旦なくしてから再度読み込みました
あと日付のデータ型を日付型に変換します
> data$日付<-as.Date(data$日付)
株価のリターンとしては前日比を用いることにして、ヒストグラムをかきます
リターンは収益率にするべきでした
リターン=前日比/前日の終値
とします
> return <- data$前日比[-1]/data$終値[1:length(data$終値)-1]
[-1]で第1成分以外を指定し、[1:length(data$終値)-1]で最後の成分以外を指定します
リターンのヒストグラムは
> hist(return,main="Histogram of Return",xlab="return",freq=FALSE,breaks=50)
freq=FALSEとすることで縦軸を頻度ではなく確率密度にしています
breaksは分割数を指定しています
さて、パッと見正規分布に見えますがどうでしょう
> mu=mean(data$前日比) > sigma=sd(data$前日比) > curve(dnorm(x,mu,sigma),add=T)
dnorm(x,mu,sigma)で平均mu分散sigmaの正規分布を作ってヒストグラムに重ねてみます
これを見ると日経平均先物は正規分布と比べて端の方に分布が偏っているように見えます
実際、分散は0.0145くらいなので6σ程度まで分布が存在することになります
6σといえば1億分の1以下なので正規分布だとまあほぼほぼ出現しないはずですが、実際にはその辺まで現れます
暴騰暴落が起こりやすいということでその辺に気をつけないと死んでしまいそうです
Rによる統計解析 日経平均先物その1
Rで統計解析を行ってみます(参考文献 : Rによる統計解析)
とは言ってもHello World!から行きましょう
まずはRのインストール(Macでやってます)
The Comprehensive R Archive Network
インストールしたら
$ R R version 3.3.0 (2016-05-03) -- "Supposedly Educational" Copyright (C) 2016 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin13.4.0 (64-bit) R は、自由なソフトウェアであり、「完全に無保証」です。 一定の条件に従えば、自由にこれを再配布することができます。 配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 R は多くの貢献者による共同プロジェクトです。 詳しくは 'contributors()' と入力してください。 また、R や R のパッケージを出版物で引用する際の形式については 'citation()' と入力してください。 'demo()' と入力すればデモをみることができます。 'help()' とすればオンラインヘルプが出ます。 'help.start()' で HTML ブラウザによるヘルプがみられます。 'q()' と入力すれば R を終了します。 >
で起動できる(簡単)
では早速
> cat('Hello World!\n') Hello World!
無事Hello World!できました
具体的に何かのテーマを決めてやった方がわかりやすそうなので
日経平均先物のデータを使ってみたいと思います
データはこちらから
http://k-db.com/futures/
今回は2016年09月限の6/17の5分足データを扱ってみます
こいつをcsvでダウンロードしてきます。
このデータに関してですが、とりあえず6/17に取引された値段が5分ごとに記録されたものです。
Rにcsvファイルを読み込ませるにはread.tableを用います
> ( data <- read.table("futures_F101-1609_5m_20160617.csv", header=TRUE, sep=",") )
ここで
headerは1行目の"日付 時刻 始値 高値 安値 終値 出来高 売買代金"をヘッダとして扱うことを意味します
またカンマ区切りなので sep="," としてやります
Rはデータフレームというリスト構造を持っており、read.tableによってデータフレームdataにcsv内のデータを構造を保ったまま保存することができます
ここ大事
>> ( data <- read.table("futures_F101-1609_5m_20160617.csv", header=TRUE,sep=",") ) make.names(col.names, unique = TRUE) でエラー: '<93><fa><95>t' に不正なマルチバイト文字があります
とエラーになってしまいます
Rが想定している文字コードがMacではutf-8なんですがhttp://k-db.com/futures/のファイルの文字コードがcp932であることに起因します
そのため、encodeしてやる必要があります
> filename <- file("futures_F101-1609_5m_20160617.csv", open="r", encoding="cp932") > ( data <- read.table(filename, header=TRUE, sep=",") ) 日付 時刻 始値 高値 安値 終値 出来高 売買代金 1 2016-06-18 03:00 15580 15580 15580 15580 82 1277560 ...
これで準備ができました。
日付や高値などは変数として扱うことができます
データフレーム名$変数名
データフレーム名 [, 列番号]
データフレーム名 [列番号]
によって変数にアクセスすることができます
例えば
> data$安値 [1] 15580 NA 15560 15560 15550 15540 15540 15540 15540 15540 15530 15500 [13] 15500 15500 15490 15490 15510 15510 15510 15500 15500 15500 15500 15480 [25] 15470 15450 15460 15450 15440 15430 15410 15420 15440 15430 15420 15440 [37] 15440 15470 15460 15470 15460 15440 15440 15460 15490 15510 15500 15490 [49] 15480 15490 15500 15490 15480 15480 15480 15510 15500 15500 15500 15500 [61] 15510 15520 15520 15510 15510 15520 15520 15530 15540 15540 15530 15530 [73] 15530 15520 15510 15510 15520 15510 15500 15510 15500 15500 15520 15530 [85] 15530 15530 15530 15530 15530 15540 15550 NA 15540 15530 15520 15530 [97] 15530 15530 15550 15560 15560 15560 15550 15540 15530 15540 15540 15530 [109] 15530 15530 15510 15510 15520 15510 15490 15480 15470 15480 15460 15480 [121] 15480 15490 15540 15540 15530 15530 15540 15530 NA 15570 15540 15560 [133] 15550 15550 15570 15570 15550 15540 15550 15560 15550 15530 15540 15560 [145] 15570 15580 15580 15580 15580 15570 15560 15550 15560 15560 15580 15580 [157] 15610 15620 15610 15630 15640 15680 15680 15680 15670 15640 15640 15640 [169] 15640 15640 15630 15630 15640 15630 15620 15620 15620 15590 15600 15590 [181] 15560 15640 15630 15610 15610 15640 15650 15680 15700 15700 15690 15680 [193] 15660 15620 15610 15630 15630 15630 15620 15590 15620 15630 15610
いくつか NA となっていますが、これはデータの欠損を示しています
日経平均先物では取引が行われない時間帯があります
www.rakuten-sec.co.jp
実際、取引が行われない時間帯のデータが欠損していることがわかります
今は問題ないですが、実際に統計値を求めるときにはNAの扱いに注意が必要になります
次に読み込んだデータをプロットしてみようと思います
日付と時刻データが分離しているため、単純にプロットすると分かりづらくなりそうです(例えば横軸を時刻にすると日付の順がおかしくなるかもしれない)
そこで日付と時刻をくっつけてみます
Rでは日付型(年月日)とPOSIXct,POSIXlt,POSIXt(日付時分秒)があります
POSIXctは1970年元旦からの経過秒で、POSIXltは可読性が高くなるように経過秒を日付時分秒の文字列リストにしたものです
POSIXtはその両方を含んでいます
今回はpaste関数で日付と時刻を結合し、as.POSIXlt関数で日付時分秒の文字列リストに変換します
> data$時刻 <- as.POSIXlt(paste(data$日付,data$時刻))
これで日付データもいらないのでデータフレームから削除します
> data <- data[, colnames(data) != "日付"]
colnamesが"日付"でない部分集合をdataに再代入しています
> data 時刻 始値 高値 安値 終値 出来高 売買代金 1 2016-06-18 03:00:00 15580 15580 15580 15580 82 1277560 ...
こんな感じになります
ここまできたら終値をプロットできます
> plot(data$時刻, data$終値, type="l", xlab = "time", ylab = "closing price")
簡素なグラフが出来上がりました
最後にsummary関数を利用した簡単なデータ解析をしてみたいと思います
summary(データフレーム名[列番号の範囲])
としてやります
実際、日付や時刻の平均や最大値は無意味なので始値から売買代金までのデータをまとめてみます
> summary(data[2:7]) 始値 高値 安値 終値 Min. :15420 Min. :15440 Min. :15410 Min. :15420 1st Qu.:15510 1st Qu.:15520 1st Qu.:15500 1st Qu.:15510 Median :15540 Median :15550 Median :15530 Median :15540 Mean :15554 Mean :15564 Mean :15544 Mean :15554 3rd Qu.:15590 3rd Qu.:15610 3rd Qu.:15580 3rd Qu.:15590 Max. :15720 Max. :15740 Max. :15700 Max. :15730 NA's :3 NA's :3 NA's :3 NA's :3 出来高 売買代金 Min. : 0.0 Min. : 0 1st Qu.: 88.5 1st Qu.: 1372845 Median : 166.0 Median : 2579150 Mean : 356.0 Mean : 5549347 3rd Qu.: 403.5 3rd Qu.: 6292590 Max. :5224.0 Max. :81704200
Min.は最小値、1st Qu.は第1四分位数、Medianは中央値、Meanは平均値、3rd Qu.は第3四分位数、Max.は最大値となっています
また、NA'sはNAの出現回数を表しています
summary関数を用いることで簡単な統計量を求めることができました