vroomパッケージを試してみたが……

追記あり

Twitterで話題になっていたcsvファイルなどを高速に読み込むR用パッケージ vroom を試してみました。


github.com


結論から言うと自分の環境ではdata.table::fread()の方が速かったです。


id:songcunyouzai さんの記事のデータを使わせていただきます。
y-mattu.hatenablog.com

library(microbenchmark)
library(ggplot2)
file <- "Sales.csv"

compare <- microbenchmark('readr::read_csv' = {gc(); gc(); readr::read_csv(file)},
                          'data.table::fread' = {gc(); gc(); data.table::fread(file)},
                          'vroom::vroom' = {gc(); gc(); vroom::vroom(file, progress = FALSE)},
                          times = 100
                          )
compare
autoplot(compare)

devtools::session_info()
結果
expr min lq mean median uq max neval cld
readr::read_csv 1024.9577 1033.3797 1085.3734 1037.4678 1155.311 1310.8491 100 b
data.table::fread 808.9454 816.2527 864.3883 830.7628 936.749 979.5104 100 a
vroom::vroom 1141.2515 1197.7813 1247.7234 1232.6876 1318.041 1361.9606 100 c


f:id:bob3:20190518101809p:plain


うーん、vroom全然速くない……


なんでこうなるんだろう?


ちょっと気になったのは id:y-mattu さんはMacで私はWindowsという点。

OSによる挙動の違いだったりするのだろうか?

session_info()
  • Session info --------------------------------------------------------------------------------------------------
setting value version R version 3.6.0 (2019-04-26) os Windows 10 x64 system x86_64, mingw32 ui RStudio language (EN) collate Japanese_Japan.932 ctype Japanese_Japan.932 tz Asia/Tokyo date 2019-05-18
  • Packages ------------------------------------------------------------------------------------------------------
package * version date lib source assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) backports 1.1.4 2019-04-10 [1] CRAN (R 3.6.0) callr 3.2.0 2019-03-15 [1] CRAN (R 3.6.0) cli 1.1.0 2019-03-19 [1] CRAN (R 3.6.0) codetools 0.2-16 2018-12-24 [1] CRAN (R 3.6.0) colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0) crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) data.table 1.12.2 2019-04-07 [1] CRAN (R 3.6.0) desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.0) devtools 2.0.2 2019-04-08 [1] CRAN (R 3.6.0) digest 0.6.18 2018-10-10 [1] CRAN (R 3.6.0) dplyr 0.8.1 2019-05-14 [1] CRAN (R 3.6.0) fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) ggplot2 * 3.1.1 2019-04-07 [1] CRAN (R 3.6.0) glue 1.3.1 2019-03-12 [1] CRAN (R 3.6.0) gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0) highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) hms 0.4.2 2018-03-10 [1] CRAN (R 3.6.0) knitr 1.22 2019-03-08 [1] CRAN (R 3.6.0) lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.0) lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0) magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0) MASS 7.3-51.4 2019-03-31 [1] CRAN (R 3.6.0) Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.6.0) memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0) microbenchmark * 1.4-6 2018-10-18 [1] CRAN (R 3.6.0) multcomp 1.4-10 2019-03-05 [1] CRAN (R 3.6.0) munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0) mvtnorm 1.0-10 2019-03-05 [1] CRAN (R 3.6.0) pillar 1.4.0 2019-05-11 [1] CRAN (R 3.6.0) pkgbuild 1.0.3 2019-03-20 [1] CRAN (R 3.6.0) pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.6.0) pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0) plyr 1.8.4 2016-06-08 [1] CRAN (R 3.6.0) prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.6.0) processx 3.3.1 2019-05-08 [1] CRAN (R 3.6.0) ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0) purrr 0.3.2 2019-03-15 [1] CRAN (R 3.6.0) R6 2.4.0 2019-02-14 [1] CRAN (R 3.6.0) Rcpp 1.0.1 2019-03-17 [1] CRAN (R 3.6.0) readr 1.3.1 2018-12-21 [1] CRAN (R 3.6.0) remotes 2.0.4 2019-04-10 [1] CRAN (R 3.6.0) rlang 0.3.4 2019-04-07 [1] CRAN (R 3.6.0) rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0) rstudioapi 0.10 2019-03-19 [1] CRAN (R 3.6.0) sandwich 2.5-1 2019-04-06 [1] CRAN (R 3.6.0) scales 1.0.0 2018-08-09 [1] CRAN (R 3.6.0) sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) survival 2.44-1.1 2019-04-01 [1] CRAN (R 3.6.0) testthat 2.1.1 2019-04-23 [1] CRAN (R 3.6.0) TH.data 1.0-10 2019-01-21 [1] CRAN (R 3.6.0) tibble 2.1.1 2019-03-16 [1] CRAN (R 3.6.0) tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.6.0) usethis 1.5.0 2019-04-07 [1] CRAN (R 3.6.0) vroom 1.0.1 2019-05-14 [1] CRAN (R 3.6.0) withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) xfun 0.7 2019-05-14 [1] CRAN (R 3.6.0) zoo 1.8-5 2019-03-21 [1] CRAN (R 3.6.0)
Vrooom Vrooom

余談ですが、King Crimson に Vrooom と Vrooom Vrooom という曲があります。
名曲なのでみんな聴いてね!

www.youtube.com

追記


というご指摘があったので型を指定して試してみました。

library(microbenchmark)
library(ggplot2)
file <- "Sales.csv"

compare <- microbenchmark('readr::read_csv' = {gc(); gc(); readr::read_csv(file)},
                          'data.table::fread' = {gc(); gc(); data.table::fread(file)},
                          'vroom::vroom' = {gc(); gc(); vroom::vroom(file, progress = FALSE)},
                          'readr::read_csv(col_types)' = {gc(); gc();
                            readr::read_csv(file,
                                            col_types = c(UserID = "c", ProductID = "i", Timestamp = "c"))},
                          'data.table::fread(colClasses)' = {gc(); gc();
                            data.table::fread(file,
                                              colClasses=c(UserID = "character", ProductID = "integer", Timestamp = "character"))},
                          'vroom::vroom(col_types)' = {gc(); gc();
                            vroom::vroom(file, progress = FALSE,
                                         col_types = c(UserID = "c", ProductID = "i", Timestamp = "c"))},
                          times = 100
                          )
compare
autoplot(compare)
expr min lq mean median uq max neval cld
readr::read_csv 1481.9142 1495.2920 1584.0018 1516.5437 1660.4270 2282.7862 100 d
data.table::fread 1141.0688 1150.7112 1183.2774 1162.5873 1199.1071 1414.1438 100 b
vroom::vroom 1436.1859 1487.1280 1506.6311 1498.5312 1517.7352 1628.3527 100 c
readr::read_csv(col_types) 1420.3024 1433.1942 1503.9445 1450.5208 1538.2399 1892.8654 100 c
data.table::fread(colClasses) 1137.5738 1151.8764 1184.1724 1173.6862 1205.0418 1293.7130 100 b
vroom::vroom(col_types) 654.5454 659.3483 691.8653 677.3639 712.9926 849.8384 100 a


f:id:bob3:20190518131720p:plain


ということで、型を指定するとvroomがグッと速くなります。


アドホックな分析の場合、型を調べて、指定して……というのは手間なのでfreadでよさそうですが、定型データを定期的に分析するような場合はvroomを活用するのがよさそうです。

lavaanで作ったモデルからパス図を描く関数

Rで構造方程式モデリング(共分散構造分析)を実行する場合、semパッケージかlavaanパッケージを使うことが多いです。

socialsciences.mcmaster.ca

lavaan.ugent.be

構造方程式モデリングを実施したらパス図が書きたくなります。
semパッケージやlavaanパッケージで作ったモデルからパス図を描くパッケージとしてsemPlotパッケージがありますが、なんかしっくりこない。
sachaepskamp.com

f:id:bob3:20190226183044p:plain
semPlot
細かく調整すればもうちょいましになるんでしょうが、面倒くさい……


そこで自分で lavaan のモデルからDiagrammeRを利用してキレイなパス図を描く関数を作りました。
rich-iannone.github.io

my_lavaan_plot <- function(fit, standardized = TRUE) {
  require(lavaan)
  require(stringr)
  require(dplyr)
  require(DiagrammeR)
  
  #描画に必要なパラメータをdfとして取り出す
  pars <- parameterEstimates(fit,
                             standardized = standardized)
  # エッジ(パス)の設定、"=~"を抜き出しpaseteして因子負荷量
  temp01 <- dplyr::filter(pars, op == "=~")
  temp02 <- str_c('  \"', temp01$lhs,'\"', ' -> ','\"', temp01$rhs, '\"',
                  '[label=\"', round(temp01$std.all,2),
                  '\" color=blue penwidth=1.001];')
  # "~"を抜き出しpaseteして因子間回帰
  temp31 <- dplyr::filter(pars, op == "~")
  temp32 <- str_c('  \"', temp31$rhs,'\"', ' -> ','\"', temp31$lhs, '\"',
                  '[label=\"', round(temp31$std.all,2),
                  '\" color=black penwidth=1.001];')
  # "~~"抜き出す
  temp11 <- dplyr::filter(pars, op == "~~")
  temp12 <- str_c('  \"', temp11$rhs,'\"', ' -> ','\"', temp11$lhs, '\"',
                  '[label=\"', round(temp11$std.all,2),
                  '\" color=glay penwidth=1.001, dir = both];')
  # ノードの設定 観測変数と潜在因子
  temp00 <- inspect(fit)$lambda
  test.f <- colnames(temp00) #潜在因子
  test.v <- rownames(temp00) #観測変数
  temp21 <- str_c('  \"', test.f,'\" [fontname=\"Helvetica\" fontsize=14 fillcolor = green shape=ellipse style=filled];')
  temp22 <- str_c('  \"', test.v,'\" [fontname=\"Helvetica\" fontsize=14 fillcolor=\"transparent\" shape=box style=filled];')
  tempALL <- c('digraph \"res.s\" {',
               '  layout = dot;', #dot, fdp
               '  rankdir=LR;',
               '  size=\"8,8\";',
               temp21,
               temp22,
               temp02,
               temp32,
               '  center=1;',
               temp12,
               '}')
  
  grViz(tempALL)
}

使い方はlavaanで作ったモデルを食わせるだけ。

# テスト用のモデル
# ヘッドスタート計画の相関係数行列
HS <- matrix(c(1.000, 0.484, 0.224, 0.268, 0.230, 0.265,
               0.484, 1.000, 0.342, 0.215, 0.215, 0.295,
               0.224, 0.342, 1.000, 0.387, 0.196, 0.234,
               0.268, 0.215, 0.387, 1.000, 0.115, 0.162,
               0.230, 0.215, 0.196, 0.115, 1.000, 0.635,
               0.265, 0.295, 0.234, 0.162, 0.635, 1.000), 6)
colnames(HS) <- rownames(HS) <- paste("x", 1:6, sep="")

HS.mdl.2 <- '
  # 測定方程式
  f1 =~ x1 + x2
  f2 =~ x3 + x4
  f3 =~ x5 + x6
  # 構造方程式
  f3 ~  f1 + f2 '
fit <- sem(HS.mdl.2, sample.cov=HS, sample.nobs=1800)
summary(fit, fit.measures=TRUE, standardized=TRUE)

my_lavaan_plot(fit)

f:id:bob3:20190226183425p:plain
自作関数による描画

Tunnels and Trolls における TARO、DAROの出目の確率

昔、トンネルズ&トロールズ(Tunnels & Trolls, T&T)というTRPGがあって、最近「完全版」も出たんだけど、この中にTARO, DARO と呼ばれるルールがある。

トンネルズ & トロールズ 完全版

トンネルズ & トロールズ 完全版

TAROは"Tripples Add and Roll Over"の略で「能力値決め等の際に6面体のサイコロを3個振ってその合計値を使うが、3個のサイコロの出目がすべて等しい場合(ゾロ目の場合)、もう一度3個のサイコロを振り、先の出た出目も合計値に足す。2回目もゾロ目だった場合はまた同様に振りなおして加える。ゾロ目が出なくなるまで繰り返す。」というもの。

  1. 最初の出目が「4,6,3」だった場合は、そのまま足し合わせて13。
  2. 最初の出目が「3,3,3」のゾロ目だった場合、もう一度振り直し、出目が「5,6,3」だった場合、2回分を合計して23。

同様にDAROは"Doubles Add and Roll Over"の略で「セービング・ロールなどの際に6面体のサイコロを2個振ってその合計値を使うが、2個のサイコロの出目が等しい場合(ゾロ目の場合)、もう一度2個のサイコロを振り、先の出目も合計値に足す。2回目もゾロ目だった場合はまた同様に振りなおして加える。ゾロ目が出なくなるまで繰り返す。」というもの。

  1. 最初の出目が「1,3」だった場合は、そのまま足し合わせて4。
  2. 最初の出目が「3,3」のゾロ目だった場合、もう一度振り直し、出目が「5,6」だった場合、2回分を合計して17。

で、普通にサイコロを2個振ったり3個振ったりする場合の出目の確率は単純に算出できるわけだが、TARO、DAROの場合はちょっとややこしくなる。

ちゃんと計算すれば出せるんだろうけど、ちょっと面倒くさかったのでシミュレーションで出してみた。
TARO,DAROともに100万回試行の結果から算出してます。
実用的にその出目以上が出る確率(累積確率)でまとめて表にしてあります。
T&TでGMする際の参考になれば幸いです

ちなみにTAROについてはほかにもやってる人がいました。
rpg.stackexchange.com

素の2d6の場合

期待値は7。

出目 この出目以上の確率
2 100.0%
3 97.2%
4 91.7%
5 83.3%
6 72.2%
7 58.3%
8 41.7%
9 27.8%
10 16.7%
11 8.3%
12 2.8%

DAROの場合

期待値は8.4。
セービング・ロールの場合は「出目=目標値-キャラレベル-能力値」で考えればOK。

出目 この出目以上の確率
3 100.0%
4 94.4%
5 88.9%
6 77.6%
7 66.3%
8 49.2%
9 37.7%
10 25.6%
11 19.2%
12 12.4%
13 11.4%
14 9.9%
15 8.9%
16 7.5%
17 6.4%
18 5.0%
19 4.1%
20 3.0%
21 2.4%
22 1.8%
23 1.5%
24 1.1%
25 1.0%

素の3d6の場合

期待値は10.5。

出目 この出目以上の確率
3 100.0%
4 99.5%
5 98.1%
6 95.4%
7 90.7%
8 83.8%
9 74.1%
10 62.5%
11 50.0%
12 37.5%
13 25.9%
14 16.2%
15 9.3%
16 4.6%
17 1.9%
18 0.5%

TAROの場合

期待値は10.8。

出目 この出目以上の確率
4 100.0%
5 98.6%
6 95.8%
7 91.7%
8 84.7%
9 75.0%
10 63.9%
11 51.4%
12 38.8%
13 27.6%
14 17.8%
15 10.7%
16 6.4%
17 3.5%
18 1.9%
19 1.8%
20 1.7%
21 1.5%
22 1.4%
23 1.2%
24 1.0%

何かのお役に立ては幸いです。

Rで表情解析

 Rで表情解析をします。

Microsoft の Emotion API をRから動かす形になります。

azure.microsoft.com

Emotion API は画像の中の人物の表情を入力として取り、Face API を使って画像の中の顔それぞれについて一連の感情の信頼度と、顔の境界ボックスを返します。ユーザーがすでに Face API を呼び出している場合は、オプション入力として顔矩形を指定することができます。
検出される感情は、怒り、軽蔑、嫌悪感、恐怖、喜び、中立、悲しみ、驚きです。これらの感情は、文化が異なっても特定の表情を伴って広く交わされると理解されています。

ということで、画像から人物の顔を抽出して、その表情を8つの指標で評価してくれます。
対応している画像は、4MB未満のJPEGPNG、GIF、BMPとのこと。
Emotion API を利用するにはAPIキーが必要になるので、前掲のページから取得しておいてください。

ではやってみましょう。

# 必要なライブラリを読み込む
library(tidyverse)
library(httr)
library(EBImage)
library(grid)

# Microsoft's Azure Emotion API を使う準備
API_url <- "https://westus.api.cognitive.microsoft.com/emotion/v1.0/recognize"
KEY <- "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" #ここに各自のAPIキーを入力

最初のモデルはモナ・リザさん。
モナ・リザ - Wikipedia

f:id:bob3:20171217124134j:plain

res <- POST(url = API_url,
            add_headers(.headers = c("Ocp-Apim-Subscription-Key" = KEY)),
            body = '{"url":"https://upload.wikimedia.org/wikipedia/commons/thumb/6/6a/Mona_Lisa.jpg/317px-Mona_Lisa.jpg"}', #ここに画像のURL
            query = list(returnFaceAttributes = "emotion"),
            accept_json())
raw <- res %>%
       content()

ここで抽出されるのは検出された顔の領域と八つの表情の指標です。
これをもとにまずは顔を抜き出します。

# 顔の領域抽出と切り出し
faceRect <- raw %>%
            pluck(1, "faceRectangle")
Left      <- faceRect$left
LeftWidth <- faceRect$left + faceRect$width
Top       <- faceRect$top
TopHeight <- faceRect$top + faceRect$height

pic <- readImage("https://upload.wikimedia.org/wikipedia/commons/thumb/6/6a/Mona_Lisa.jpg/317px-Mona_Lisa.jpg") #ここも画像のURL
pic_face <- pic[Left:LeftWidth, Top:TopHeight,]
plot(pic_face)

抜き出された顔がこちら。
f:id:bob3:20171217125109j:plain:medium
顔ですね。

では、モナ・リザさんのご機嫌は?

# 感情指標抽出
scores <- raw %>%
          pluck(1, "scores") %>%   
          as.tibble() %>%
          t() %>%
          round(4)
scores <- tibble(emotion=c("怒り","侮蔑","嫌悪","恐れ","幸せ","普通","悲しみ","驚き"),
                 scores=scores[,1])

# 感情のチャート化
g <- ggplot(scores, aes(x=emotion, y=scores, fill=emotion)) +
     geom_bar(stat = "identity") +
     ylim(0, 1) +
     annotation_raster(pic_face, 4, 5, 0.875, 1,interpolate=TRUE)
plot(g)

f:id:bob3:20171217130257p:plain
普通に幸せってことでまさに「微笑」ですね。

続けてほかの画像の結果も。

マリリンモンローのびっくり顔。
File:Gentlemen Prefer Blondes Movie Trailer Screenshot (16-2).jpg - Wikimedia Commons
f:id:bob3:20171217141521j:plain:medium
f:id:bob3:20171217141555p:plain
びっくりしてますね。

映画『ベン・ハー』より。
File:Jack Hawkins in Ben Hur trailer.jpg - Wikimedia Commons
f:id:bob3:20160910132332j:plain:medium
f:id:bob3:20171217142759p:plain
複雑な表情です。

映画『風と共に去りぬ』より
File:De Havilland-Melanie.jpg - Wikimedia Commons
f:id:bob3:20171217144029j:plain:medium
f:id:bob3:20171217144303p:plain
幸せそうです。

同上。
File:Clark Gable as Rhett Butler in Gone With the Wind trailer.jpg - Wikimedia Commons
f:id:bob3:20170705140520j:plain:medium
f:id:bob3:20171217145625p:plain
ちょっとスケベっぽいけど。

以上です。

enjoy!

TUBE-01Jと真空管(6J1, 6AK5, 5654)

f:id:bob3:20161204234141j:plain
f:id:bob3:20161204234854j:plain
真空管アンプ(プリアンプ)TUBE-01Jを購入しまして、秋葉原で互換球を探して来ました。

写真の左から、
日立 5654 :60年代製のヴィンテージ。
東芝 6AK5(通測用Hi-S):年代不明だが70年代か?
GE JAN 5654W:軍事用の白箱。米国産。
80年代製と思われる。
R.T.C. 5654RT:こちらも軍事用。フランスの企業だがオランダ産でPhilipsのOEMらしい。
曙光 6J1:TUBE-01Jの標準球。中国産。