Rを使ってExcelファイルを読み込んでcsvファイルを出力するまで
那須野さんの記事が興味深かったので、Rでもやってみた。
takuminasuno.com
library(tidyverse) #データ操作用 library(readxl) #Excelファイルの読み込み用 library(stringi) #文字列操作用 # 列名の指定。空列は"blank"とする。 nameList <- c("age", "blank1", "blank2", "total_201905", "male_201905", "female_201905", "blank3", "total_201812", "male_201812", "female_201812", "total-ja_201812", "male-ja_201812", "female-ja_201812") # Excelファイルの読み込み。読み込む個所はrangeで指定。 # index(age)の補正。空白の除去と英数の半角化。「歳」の追加。 # 数値の桁数の補正(単位を揃える)。2019年は1万倍、2018年は千倍。 dat <- read_excel("05k2-1.xls", range = "A14:M34", col_names = nameList) %>% select(-starts_with("blank")) %>% mutate(age = age %>% str_replace_all("[[:space:]]", "") %>% stri_trans_nfkc() %>% if_else(str_detect(., pattern="歳"), ., str_c(., "歳")) ) %>% mutate_at(vars(contains("2019")), ~.*10000) %>% mutate_at(vars(contains("2018")), ~.* 1000) # 横持ちでcsv出力 write_excel_csv(dat, "output_pivotted.csv") # 縦持ちのデータに変換 dat2 <- dat %>% gather(key="sex_yyyymm", value = "value", -age) %>% separate(sex_yyyymm, into = c("sex", "yyyymm"), sep = "_") # 縦持ちでcsv出力 write_excel_csv(dat2, "output_unpivotted.csv")
まあ、RでもPythonでも慣れてる方で書けばいいのかな。
vroomパッケージを試してみたが……
追記あり。
Twitterで話題になっていたcsvファイルなどを高速に読み込むR用パッケージ vroom を試してみました。
結論から言うと自分の環境では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 |
うーん、vroom全然速くない……
なんでこうなるんだろう?
ちょっと気になったのは id:y-mattu さんはMacで私はWindowsという点。
OSによる挙動の違いだったりするのだろうか?
session_info()
- Session info --------------------------------------------------------------------------------------------------
- Packages ------------------------------------------------------------------------------------------------------
追記
もしかしたら col_types を指定すれば速くなるとか
— Haruhiko Okumura (@h_okumura) May 18, 2019
というご指摘があったので型を指定して試してみました。
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 |
ということで、型を指定するとvroomがグッと速くなります。
アドホックな分析の場合、型を調べて、指定して……というのは手間なのでfreadでよさそうですが、定型データを定期的に分析するような場合はvroomを活用するのがよさそうです。
lavaanで作ったモデルからパス図を描く関数
Rで構造方程式モデリング(共分散構造分析)を実行する場合、semパッケージかlavaanパッケージを使うことが多いです。
構造方程式モデリングを実施したらパス図が書きたくなります。
semパッケージやlavaanパッケージで作ったモデルからパス図を描くパッケージとしてsemPlotパッケージがありますが、なんかしっくりこない。
sachaepskamp.com
細かく調整すればもうちょいましになるんでしょうが、面倒くさい……
そこで自分で 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)
Tunnels and Trolls における TARO、DAROの出目の確率
昔、トンネルズ&トロールズ(Tunnels & Trolls, T&T)というTRPGがあって、最近「完全版」も出たんだけど、この中にTARO, DARO と呼ばれるルールがある。
- 出版社/メーカー: cosaic
- 発売日: 2016/09/24
- メディア: おもちゃ&ホビー
- この商品を含むブログ (7件) を見る
TAROは"Tripples Add and Roll Over"の略で「能力値決め等の際に6面体のサイコロを3個振ってその合計値を使うが、3個のサイコロの出目がすべて等しい場合(ゾロ目の場合)、もう一度3個のサイコロを振り、先の出た出目も合計値に足す。2回目もゾロ目だった場合はまた同様に振りなおして加える。ゾロ目が出なくなるまで繰り返す。」というもの。
- 最初の出目が「4,6,3」だった場合は、そのまま足し合わせて13。
- 最初の出目が「3,3,3」のゾロ目だった場合、もう一度振り直し、出目が「5,6,3」だった場合、2回分を合計して23。
同様にDAROは"Doubles Add and Roll Over"の略で「セービング・ロールなどの際に6面体のサイコロを2個振ってその合計値を使うが、2個のサイコロの出目が等しい場合(ゾロ目の場合)、もう一度2個のサイコロを振り、先の出目も合計値に足す。2回目もゾロ目だった場合はまた同様に振りなおして加える。ゾロ目が出なくなるまで繰り返す。」というもの。
- 最初の出目が「1,3」だった場合は、そのまま足し合わせて4。
- 最初の出目が「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から動かす形になります。
Emotion API は画像の中の人物の表情を入力として取り、Face API を使って画像の中の顔それぞれについて一連の感情の信頼度と、顔の境界ボックスを返します。ユーザーがすでに Face API を呼び出している場合は、オプション入力として顔矩形を指定することができます。
検出される感情は、怒り、軽蔑、嫌悪感、恐怖、喜び、中立、悲しみ、驚きです。これらの感情は、文化が異なっても特定の表情を伴って広く交わされると理解されています。
ということで、画像から人物の顔を抽出して、その表情を8つの指標で評価してくれます。
対応している画像は、4MB未満のJPEG、PNG、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
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)
抜き出された顔がこちら。
顔ですね。
では、モナ・リザさんのご機嫌は?
# 感情指標抽出 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)
普通に幸せってことでまさに「微笑」ですね。
続けてほかの画像の結果も。
マリリンモンローのびっくり顔。
File:Gentlemen Prefer Blondes Movie Trailer Screenshot (16-2).jpg - Wikimedia Commons
びっくりしてますね。
映画『ベン・ハー』より。
File:Jack Hawkins in Ben Hur trailer.jpg - Wikimedia Commons
複雑な表情です。
映画『風と共に去りぬ』より
File:De Havilland-Melanie.jpg - Wikimedia Commons
幸せそうです。
同上。
File:Clark Gable as Rhett Butler in Gone With the Wind trailer.jpg - Wikimedia Commons
ちょっとスケベっぽいけど。
以上です。
enjoy!
TUBE-01Jと真空管(6J1, 6AK5, 5654)
真空管アンプ(プリアンプ)TUBE-01Jを購入しまして、秋葉原で互換球を探して来ました。
FX-AUDIO- TUBE-01J【シルバー】真空管プリアンプ(ラインアンプ)NFJ オリジナルモデル
- 出版社/メーカー: North Flat Japan
- メディア: エレクトロニクス
- この商品を含むブログを見る
写真の左から、
日立 5654 :60年代製のヴィンテージ。
東芝 6AK5(通測用Hi-S):年代不明だが70年代か?
GE JAN 5654W:軍事用の白箱。米国産。
80年代製と思われる。
R.T.C. 5654RT:こちらも軍事用。フランスの企業だがオランダ産でPhilipsのOEMらしい。
曙光 6J1:TUBE-01Jの標準球。中国産。
テスト
テスト
library(MASS) truehist(iris[,1])