Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2014-08-24

やってみよう分析!Rで強化学習(Q-learning, ε-greedy行動選択)

| 06:57 | やってみよう分析!Rで強化学習(Q-learning, ε-greedy行動選択)を含むブックマーク

のR実装版。強化学習自体の解説は上の記事読んどいたらいい。めんどいのでとりあえずε-greedyのみやった。計算結果は

> Qlearning()
     [,1]     [,2]
[1,]    5 20.00000
[2,]   10  4.99903

となって、上のLINK先と同じなのでコードは間違ってなさそうだ。

以下、全コード。めんどいのでノーコメント。後で昔の俺が未来の俺の首を絞める可能性大だ。

#行動選択
selectAction <- function(s, numA, Qt)
{
  qa <- Qt[s,]
  max_indexes <- (1:length(qa))[qa==max(qa)]
  if(length(max_indexes) == 1){
    max_indexes
  }else{
    max_indexes[sample(length(max_indexes), 1)]
  }
}
epsilonGreedy <- function(epsilon, numA, s, Qt)
{
  if(runif(1) < epsilon)
  {
    sample.int(numA, 1)
  }else{
    selectAction(s, numA, Qt)
  }
}
nextState <- function(s, a)
{  
  ifelse(a==1, ifelse(s==1, 2, 1), s)
}
calcReward <- function(s, a, net)
{
  ifelse(a!=1 && s==1, net, 0)
}
Qlearning <- function()
{
  numA <- 2
  numS <- 2
  alpha <- 0.3
  gamma <- 0.5
  epsilon <- 0.1
  trialMax <- 10000
  stepMax <- 10
  net <- 10
  s = 1
  Qt <- matrix(0, numS, numA)
  
  for(i in 1:trialMax)
  {
    for(j in 1:stepMax)
    {
      a <- epsilonGreedy(epsilon, numA, s, Qt)
      sd <- nextState(s, a)
      reward <- calcReward(s, a, net)
      Qt[s, a] <- (1 - alpha)*Qt[s, a] + alpha*(reward + gamma*max(Qt[sd,]))

      if(reward > 0){
        s <- 1
        break
      } else{
        s <- sd        
      }
    }
  }
  Qt
}

2014-08-20

殺伐とした空気を吹き飛ばすアスキーアートなテラモナギ

| 21:22 | 殺伐とした空気を吹き飛ばすアスキーアートなテラモナギを含むブックマーク

なんか殺伐とした空気を吹き飛ばせるんじゃないか、そう強く感じさせるパッケージを2つみつけた。

まず、Githubからのインストール

library(devtools) 
install_github("wrathematics/Rfiglet")
install_github("sckott/cowsay")

Rfigletから使ってみる。

> library(Rfiglet)
> figlet("teramonagi")
 _                                                   _ 
| |_ ___ _ __ __ _ _ __ ___   ___  _ __   __ _  __ _(_)
| __/ _ \ '__/ _` | '_ ` _ \ / _ \| '_ \ / _` |/ _` | |
| ||  __/ | | (_| | | | | | | (_) | | | | (_| | (_| | |
 \__\___|_|  \__,_|_| |_| |_|\___/|_| |_|\__,_|\__, |_|
                                               |___/    
> figlet("teramonagi", font="broadway")
                                                                        
8888888 8888888888 8 8888888888   8 888888888o.            .8.          
      8 8888       8 8888         8 8888    `88.          .888.         
      8 8888       8 8888         8 8888     `88         :88888.        
      8 8888       8 8888         8 8888     ,88        . `88888.       
      8 8888       8 888888888888 8 8888.   ,88'       .8. `88888.      
      8 8888       8 8888         8 888888888P'       .8`8. `88888.     
      8 8888       8 8888         8 8888`8b          .8' `8. `88888.    
      8 8888       8 8888         8 8888 `8b.       .8'   `8. `88888.   
      8 8888       8 8888         8 8888   `8b.    .888888888. `88888.  
      8 8888       8 888888888888 8 8888     `88. .8'       `8. `88888. 
          .         .                                            
         ,8.       ,8.           ,o888888o.     b.             8 
        ,888.     ,888.       . 8888     `88.   888o.          8 
       .`8888.   .`8888.     ,8 8888       `8b  Y88888o.       8 
      ,8.`8888. ,8.`8888.    88 8888        `8b .`Y888888o.    8 
     ,8'8.`8888,8^8.`8888.   88 8888         88 8o. `Y888888o. 8 
    ,8' `8.`8888' `8.`8888.  88 8888         88 8`Y8o. `Y88888o8 
   ,8'   `8.`88'   `8.`8888. 88 8888        ,8P 8   `Y8o. `Y8888 
  ,8'     `8.`'     `8.`8888.`8 8888       ,8P  8      `Y8o. `Y8 
 ,8'       `8        `8.`8888.` 8888     ,88'   8         `Y8o.` 
,8'         `         `8.`8888.  `8888888P'     8            `Yo 
                                            
         .8.           ,o888888o.    8 8888 
        .888.         8888     `88.  8 8888 
       :88888.     ,8 8888       `8. 8 8888 
      . `88888.    88 8888           8 8888 
     .8. `88888.   88 8888           8 8888 
    .8`8. `88888.  88 8888           8 8888 
   .8' `8. `88888. 88 8888   8888888 8 8888 
  .8'   `8. `88888.`8 8888       .8' 8 8888 
 .888888888. `88888.  8888     ,88'  8 8888 
.8'       `8. `88888.  `8888888P'    8 8888  

cowsayパッケージにはAAがあるみたい。WinだとSys.setlocale(category = "LC_CTYPE", locale = "C")しておかないとうまく表示されなかった。

> library("cowsay")
> Sys.setlocale(category = "LC_CTYPE", locale = "C")
[1] "C"
> say('time')

 ----- 
 2014-08-20 21:20:06 
 ------ 
    \   ^__^ 
     \  (oo)\ ________ 
        (__)\         )\ /\ 
             ||------w|
             ||      ||

元ネタと例

2014-08-19

1:length(x) の代わりに seq_along(x) を使うと良いってごみ箱が言ってた

| 20:45 | 1:length(x) の代わりに seq_along(x) を使うと良いってごみ箱が言ってたを含むブックマーク

seq_alongの方が、空っぽのベクターに対しても安全に動作しますよっと。

> x <- 1:3
> 1:length(x)
[1] 1 2 3
> seq_along(x)
[1] 1 2 3
> x <- c()
> 1:length(x)
[1] 1 0
> seq_along(x)
integer(0)

2014-08-15

地図を描きたい俺はgoogleVisパッケージかggmapパッケージかどちらを使おうか迷うが、静的な地図を描きたい場合はggmapで良いみたい

| 19:39 | 地図を描きたい俺はgoogleVisパッケージかggmapパッケージかどちらを使おうか迷うが、静的な地図を描きたい場合はggmapで良いみたいを含むブックマーク

掲題の件、そういうことです。動的にインタラクティブに地図いじりたい!って人はgoogleVisパッケージでいいけど、静的な地図ならggmapなな印象だね!・・・と以下に書いてあった。

というわけで、今、静的な地図に興味がある俺はggmapをちょいとお勉強。

日本を中心として地図を描く

基本にして、ここで紹介するget_map, ggmap関数がggmapパッケージのメイン関数となる。

get_map関数Google Mapを筆頭とする4つのデータソースから地図を取得する関数になってて、source引数でデータソースを指定する。

また、それぞれにソースに対して直接アクセスする関数

- get_googlemap

- get_openstreetmap

- get_stamenmap

- get_cloudmademap

もある。あるが、俺はgoogle mapばっかり使うことになるんだろうなって。その他のデータソースは

らっぽい。

ggmap関数がget_mapで取得したデータの描画関数にあたる。

以下は大体日本を中心とした地図の描画方法。

library(ggmap)
map <- get_googlemap(center=c(139,35), scale=1, size=c(1260, 960), zoom=4)
ggmap(map) + 
  geom_point(aes(x=lon, y=lat), data=data.frame(lon=139, lat=35), size=10, colour='black')

f:id:teramonagi:20140815190945p:image

geocodeで地名から緯度・経度算出

北海道の緯度・経度。県庁所在地なのかな?

> geocode('Hokkaido')
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Hokkaido&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
       lon      lat
1 141.3468 43.06461

revgeocode関数で逆引き(緯度・経度⇒住所)できる。

> revgeocode(c(141.3468,43.06461))
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?latlng=43.06461,141.3468&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
[1] "6 Chome-1 Kita 3 Jonishi, Chuo Ward, Sapporo, Hokkaido, Japan"

これらのAPIに対する利用制限はgeocodeQueryCheckで取得できる。

> geocodeQueryCheck()
2498 geocoding queries remaining. 

get_mapする際の引数(maptype)を変えることで、いろんな種類の地図になる

以下はロードマップの例。他にもいろいろある。

> library(magrittr)
> get_map(maptype='roadmap') %>% ggmap(map, fullpage=TRUE)

f:id:teramonagi:20140815194134p:image

ggimage + hadleyで神降臨。

ありがたや。。。

> ggimage(hadley)

f:id:teramonagi:20140815190946p:image


2点間の経路を検索してその結果を地図に描画

マニュアルから引っ張って来たそのまんまですが、あぁ、こりゃすげぇや。。。

> from <- 'houson, texas'
> to <- 'waco, texas'
> route_df <- route(from, to, structure = 'route')
Information from URL : http://maps.googleapis.com/maps/api/directions/json?origin=houson,+texas&destination=waco,+texas&mode=driving&units=metric&alternatives=false&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
> qmap('college station, texas', zoom = 8) +
+     geom_path(
+         aes(x = lon, y = lat), colour = 'red', size = 1.5,
+         data = route_df, lineend = 'round'
+     )
Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=college+station,+texas&zoom=8&size=%20640x640&scale=%202&maptype=terrain&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=college+station,+texas&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms

f:id:teramonagi:20140815194135p:image

Open Street Map/Stamenからは取れなかった

Googleじゃないソースも使ってみようかなって思ったら、なんか怒られる。

> get_openstreetmap()
Error: map grabbing failed - see details in ?get_openstreetmap.
In addition: Warning message:
In download.file(url, destfile = destfile, quiet = !messaging, mode = "wb") :
  cannot open: HTTP status was '503 Service Unavailable'
> get_stamenmap()
Error in readPNG(destfile) : file is not in PNG format

cloudmademapはAPIキーが必要だし、googlemapかなぁ。

あとこれ画像張り付けるのめんどいからRPubs + Markdownでまとめるべき事案でしたわ。。。

2014-08-14

こういうことが言いたいだけなんじゃねぇの?

| 22:59 | こういうことが言いたいだけなんじゃねぇの?を含むブックマーク

他の部分に対する突っ込みはよくわからないが、とりあえずここはこうだろうと。

↑この本を買ったのだが、95ページ周辺にむちゃくちゃな事が書いてあって萎えた。

統計学では「仮説検定」という手法があります。

しかし、仮説検定は(中略)、人数が(サンプルサイズ)が多い場合には、ほとんどの結果が「統計的有意な差がある」という結果になります。

差がわずかしかない二つの集団を比べる場合は、そんなことはない。

本当?"標準偏差1"という値に比べて十分に小さい"0.1"という間隔を持つ2つの正規分布からなるサンプル(10^i個; i∈{1,2,3,4,5,6,7})に対して、各々、平均値の差を見る検定であるt検定をかけたところ、

library(magrittr)
r <- list()
for(i in 1:7){
  set.seed(100)
  r[[length(r)+1]] <- t.test(rnorm(10^i, mean=0.1), rnorm(10^i),var.equal=T)   
}
sapply(r, function(z)z$p.value) %>% plot(1:7, ., type="b", ylab="p-value", xlab="number of sample(10^i)")

f:id:teramonagi:20140814225713p:image

となって、サンプル数がでかくなると、p値が小さくなる傾向が見て取れる、従って書籍の言うように「サンプル数が大きくなると統計的有意であると判断する確率が高まる傾向がありそうだ」といっていいように思えるんだが。

あと、みんなmagrittrパッケージ使ったらいいと思う。

カジュアルにC++11 in R

| 21:33 | カジュアルにC++11 in Rを含むブックマーク

以下のプラグイン

plugins=c("cpp11")

をcppFunctionに指定するだけ。Winでもいけた。C++11なんで戻り値をdecltypeにしてやろうと思ったら、できなかった。windowsgccが古いせいかな?

library(Rcpp)
cppFunction('
  std::vector<double> twoTimes(std::vector<double> xs) 
  {
    for(auto &x : xs)
    {
      x *= 2.0;
    }
    return xs;
  }', plugins=c("cpp11"))

こんなかんじで。

> twoTimes(1:10)
 [1]  2  4  6  8 10 12 14 16 18 20

こうなる。