マレーシア航空旅客機の衛星写真を画像解析
マレーシア航空のいろんな話題が飛び交っている中、Twitterで「【助けて】全世界のインターネットユーザーに協力を呼びかけ！ この写真から「消息を絶ったマレーシア航空の旅客機」を見つけてください」という記事を見つけました。
衛星写真の画像が大量にあるので、人海戦術で破片を見つけよう、という趣旨です。私は分析屋なので、データサイエンスを使って手助けできないか？と思い、少し分析してみました。何かの一助になればと思い、ブログで公開します。
データサイエンスを利用するには、まずデータが必要です。衛星画像はこのサイトにあります。
そのサイトから、まずは特徴的な画像を拾ってきて、分析してみます。サイトを見てみると、画像はだいたい３つのパターンに分かれているようでした。
- 雲の画像
- 海の画像
- 何かの物体の画像
3つ目の”何かの物体”を画像から判定できれば、それが旅客機かもしれません。分析のロジックを作るために、ある船が含まれる画像を使って分析してみました。その画像が次のものです。サイト上ではこのページ（23975番）の画像です。
衛星画像をいくつか適当に眺めて、たまたま見つけたものです。まずは特徴的なこの画像を使って、分析ロジックを検討します。船とヘリポートが写っており、これを画像解析によって見つけられるか試します（画像データのダウンロードはできなかったので、Web画面のスクリーンショットを保存しました）。
分析には、いつものように統計解析ソフトであるRを利用します。Rでの画像処理に関しては、@wdkzさんの「Rでウォーリーを探してみた」の記事が大変参考になりました。この場を借りてお礼申し上げます。
まず、飛行機が写っているとしたら画像では白っぽくなるだろうと思い、画像から”白っぽい物体”を抽出するロジックを考えました。いろいろ試行錯誤しましたが、結果的に、次のような順で行うことで”白っぽい物体”を抽出できました。
１．白っぽい部分以外を１色に塗りつぶす
２．白っぽい部分も１色にする
３．色を少しぼやけさせる
４．ぼやけた部分をくっきりさせる
５．ノイズを除去する
最終的に出来た画像が、５の処理で行ったものです。赤く協調されているものが”白っぽいもの”です。もともとあった船とヘリポートがしっかりと認識されています。ただ、他にもいくつか赤い点があります。もとの画像をみると、恐らく海の波か何かなんじゃないかと思いますが、ちょっと分かりません。ただ、１つ目の画像処理の時点で、何か浮き出てますね。。何だろ、、、もしかしたら目的のものの一部？
これで”白っぽいものを認識する”ロジックはできたはずなので、確認のために他の画像にも当てはめてみます。先ほどの海域の左にある、雲の画像に当てはめてみたものが、次の結果です。
もとの画像
解析後の画像
２枚目が解析後の画像ですが、雲が１つの物体と認識されてないので明らかにおかしいです。多分、判定ロジックのどこかがおかしく、まだ改良しなくてはならなさそうです。これがうまくいったら次は、「雲かどうか」も何かしらのロジックで判別しなくてはならないですね。
とりあえず今思い付く方法はこれくらいです。これらをまとめると、次の３点です。
- 画像の白色判定はうまくいったので、飛行機を見つけるときの補助ができる（人が目で見る分には、処理１か２だけで良いかもしれない）
- 【課題】雲を見分ける必要がある（雲が無い衛星写真があればベスト）
- 【課題】サイトから画像が全部ダウンロードできないので、自動で取得してくる必要がある
今できるのはこれくらいですが、改良すればもっと実用的になるかもしれません。何かコメントがあればTwitterの方でどうぞ。全画像が手に入れば、ビッグデータ解析になります。
なるべく早く事案が解決することを願っています。
以下、Rのコードです。ご参考下さい。
# ImageMagickなどのインストール（ImageMagickは結局使ってない） library(installr) install.ImageMagick() install.packages(c("abline", "biOps")) source("http://www.bioconductor.org/biocLite.R") biocLite("EBImage") source("http://rimagebook.googlecode.com/svn/installRImageBook.R") installRImageBook() require(EBImage) require(biOps) require(RImageBook) # setwd直下に、画像番号をファイル名にしたキャプチャpng画像を置いておく # setwd直下に画像の番号を付けたフォルダを作る必要がある WhiteExtract <- function(number){ MalayAir <- readImage(file=paste("MalayAir", number, ".png", sep="")) # display(MalayAir) # gtkとImageMagickの設定が難しい。writeImageで代用する。 # 白色抽出 writeImage(MalayAir, paste(number, "/MalayAirTest0.png", sep="")) MalayAir.white <- MalayAir MalayAir.white[MalayAir.white > 0.1] <- 1 writeImage(MalayAir.white, paste(number, "/MalayAirTest1.png", sep="")) MalayAir.white2 <- thresh(MalayAir.white, w=50, h=50, offset=0.1) writeImage(MalayAir.white2, paste(number, "/MalayAirTest2.png", sep="")) MalayAir.white3 <- gblur(MalayAir.white2, sigma=1) writeImage(MalayAir.white3, paste(number, "/MalayAirTest3.png", sep="")) MalayAir.white4 <- (MalayAir.white3 > 0.1) writeImage(MalayAir.white4, paste(number, "/MalayAirTest4.png", sep="")) kern <- makeBrush(7, shape="diamond") MalayAir.white5 <- closing(opening(MalayAir.white4, kern), kern) writeImage(MalayAir.white5, paste(number, "/MalayAirTest5.png", sep="")) writeImage(MalayAir.white5, paste("MalayAir", number, "_判定後.png", sep="")) MalayAir.white5.label <- bwlabel(MalayAir.white5) print(max(MalayAir.white5.label)) } WhiteExtract(21387) WhiteExtract(24889)
RでDeep Learningの一種をやってみる
このブログのTips052で、RでDeep Learningをやっているのを紹介してもらったので、自分も試してみました。
「Deep Learningすげぇ！！」という話は良く聞くのですが、亜種がいっぱいあるみたいで、まだあまり調査しきれてません。また時間があれば調査してまとめられると良いのですが。
以下、RでDA(Denoising Autoencoders)を実行するプログラムです。
sigmoid <- function(x){
return (1 / (1 + exp(-x)))
}
dA <- setRefClass(
Class="dA",
fields=list(input="matrix", n_visible="numeric",
n_hidden="numeric", W="matrix",
W.prime="matrix", hbias="vector",
vbias="vector", rng="numeric"),
methods=list(
get_corrupted_input=function(input, corruption_level){
stopifnot(corruption_level < 1) # アサーション
return(rbinom(size=1, n=ncol(input) * nrow(input),
prob=1 - corruption_level) * input) # size=1とすることで２項分布からベルヌーイ分布を得る
},
# Encode
get_hidden_values=function(input){
# print("-------------------get_hidden_values()")
# print("input=")
# print(input)
# print(".self$W=")
# print(.self$W)
# print("-------------------")
return(sigmoid(input %*% .self$W + .self$hbias))
},
# Decode
get_reconstructed_input=function(hidden){
# print("-------------------get_reconstructed_input()")
# print("hidden=")
# print(hidden)
# print(".self$W.prime=")
# print(.self$W.prime)
# print("-------------------")
.self$W.prime <- t(.self$W)
return(sigmoid(hidden %*% .self$W.prime + .self$vbias))
},
train=function(lr=0.1, corruption_level=0.3, input=NULL){
if(is.null(input)==FALSE){
.self$input <- input
}
input <- .self$input
# print("-------------------train()")
# print("input=")
# print(input)
# print("-------------------")
tilde_x <- .self$get_corrupted_input(input, corruption_level)
y <- .self$get_hidden_values(tilde_x)
z <- .self$get_reconstructed_input(y)
L_h2 <- input - z
L_h1 <- (L_h2 %*% .self$W) * y * (1 - y)
L_vbias <- L_h2
L_hbias <- L_h1
L_W <- (t(tilde_x) %*% L_h1) + (t(L_h2) %*% y)
.self$W <- .self$W + lr * L_W
.self$hbias <- .self$hbias + lr * colMeans(L_hbias)
.self$vbias <- .self$vbias + lr * colMeans(L_vbias)
# print("=============================")
# print(L_hbias)
# print(L_vbias)
# print(colMeans(L_hbias))
# print(colMeans(L_vbias))
# print("=============================")
},
negative_log_likelihood=function(corruption_level=0.3){
tilde_x <- .self$get_corrupted_input(.self$input, corruption_level)
y <- .self$get_hidden_values(tilde_x)
z <- .self$get_reconstructed_input(y)
# print("=============================")
# print(tilde_x)
# print(y)
# print(z)
# print(log(z))
# print(log(1-z))
# print("=============================")
cross_entropy = - mean(rowSums(.self$input * log(z) + (1 - .self$input) * log(1 - z)))
return(cross_entropy)
},
reconstruct=function(x){
y <- .self$get_hidden_values(x)
z <- .self$get_reconstructed_input(y)
return(z)
},
# コンストラクタ
initialize=function(input=NULL, n_visible=NULL,
n_hidden=NULL, W=NULL, W.prime=NULL,
hbias=NULL, vbias=NULL, rng=NULL){
if(is.null(input) == TRUE){
input <<- matrix(NA)
}
else{
input <<- input
}
if(is.null(n_visible) == TRUE){
n_visible <- 2
n_visible <<- n_visible
}
else{
n_visible <<- n_visible
}
if(is.null(n_hidden) == TRUE){
n_hidden <- 3
n_hidden <<- n_hidden
}
else{
n_hidden <<- n_hidden
}
if(is.null(rng) == TRUE){
rng <- 1234
rng <<- rng
}
else{
rng <<- rng
}
if(is.null(W) == TRUE){
a <- 1 / n_visible
set.seed(rng)
W <- matrix(runif(n_visible * n_hidden, min=-a, max=a),
n_visible, n_hidden)
W <<- W
W.prime <<- t(W)
}
else{
W <<- W
W.prime <<- t(W)
}
if(is.null(hbias) == TRUE){
hbias <<- rep(0, n_hidden) # initialize h bias 0
}
else{
hbias <<- hbias
}
if(is.null(vbias)==TRUE){
vbias <<- rep(0, n_visible) # initialize v bias 0
}
else{
vbias <<- vbias
}
}
)
)
test_dA <- function(learning_rate=0.1, corruption_level=0.3,
training_epochs=50){
data <- rbind(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0))
rng <- 123 # seed
# construct dA
da <- dA$new(input=data, n_visible=20, n_hidden=5, rng=rng)
# train
for(epoch in seq(training_epochs)){
da$train(lr=learning_rate, corruption_level=corruption_level)
cost = da$negative_log_likelihood(corruption_level=corruption_level)
print(paste('Training epoch ', epoch, ' cost is ', cost))
learning_rate <- learning_rate * 0.95
}
# test
x <- rbind(c(1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0))
print(da$reconstruct(x))
}
test_dA()
とりあえず動くところまでは確認（@konsonsanありがとうございます！）。
今回はダミーデータに対してやってるので、またいつかコードレビューをしつつ、サンプルデータとか実データでもやってみたい。
【CROSS2013】『活躍するデータサイエンティストの人材像』のセッションメモ
『今日から始まるデータサイエンティスト』という２部構成のセッションのうち、前半の『活躍するデータサイエンティストの人材像』の進行を行なってきました。みなさま第一線で活躍されている方ばかりで、進行をさせて頂くのが恐縮でしたが。。準備も突貫で行ったのですが、みなさまのお陰で無事に実施できました！
進行を行なっていたのでメモを取れなかったのですが、とりいそぎ流れを記録しておきます（あとで覚えている範囲で追記していければ）。
【進行】
倉橋一成（iAnalysis合同会社代表・CAO）
【パネリスト】
草野隆史（株式会社ブレインパッド社長）
益村勝将（トランスコスモス・アナリティクス株式会社COO）
佐々木智之（株式会社gumi執行役員）
西郷彰（株式会社リクルートテクノロジーズ）
【内容】
・オーガナイザーの里洋平氏（株式会社ディー・エヌ・エー）より前座トーク
・各自自己紹介
・分析系の職種を選ばれた理由は？
・データサイエンティストの人物像・スキルセットは？
・データサイエンティストを育てていくためにはどうしたらよいか？
・データサイエンティストを採用するにあたって気を付けることは？
他にも「分析を効率よく進めるためのシステム環境は？」とか「データサイエンティストを企業利益に結び付けるための組織構成は？」とか聞きたかったのですが時間がなかったです。。また機会があれば良いです。
パネリストのみなさん表現方法は違うものの、全体的に人物像は同じようなものを想像されているように感じました。私の言葉で言えば、「学問知識」「ビジネス感」「エンジニアリング」の３軸です（ホームページでも書いています）。この３軸を磨いてスキルの”面積”を大きくしつつ、企画力・提案力を身につけていくことが大事ですね。
あと印象的だったのは、いろんな意味で「バランス力」を持っている人が大事ということ。偏っていると効果の高い分析は難しいのかもしれません。あくまでも本セッションは、マネジメントの視点からなので、そういった意見が多かったのだと思います。もちろんエンジニアリングや数理モデリングに特化した方も、非常に重要です。
また思い出したら追記していきたいと思います。
データサイエンティストを目指すなら知っておきたいRパッケージ10個＋α
元ネタのブログは「10 R packages every data scientist should know about」と「10 R packages I wish I knew about earlier」です。紹介されているパッケージはどれも良いのでメモしておきます。私が「取得した方がいいだろうなー」と思う順番に並べ替えてます。サンプルコードは後者の記事に載ってます。
- randomForest：超強力な汎用予測モデル
- RPostgreSQL, RMYSQL, RMongo, RODBC, RSQLite：各種データベースへの接続
- plyr：データ集約
- reshape2：データ加工
- forecast：時系列予測
- stringr：文字列操作
- lubridate：日付操作
- sqldf：SQLライクなデータ操作
- ggplot2：綺麗なプロットを描く
- qcc：品質管理
個人的には、下の３つは優先度低いです。理由は、sqldf：R使いっぽくない、ggplot2：指定の仕方が特殊なので、結局描きたい絵を描くのが難しい（エクセルやパワポを使った方が早いことも）、qcc：使い所があまり多くない。
プラス、下記に私のオススメパッケージを紹介しておきます。
- party：決定木が綺麗に描ける
- gbm：randomForestより汎用性の高い超強力な予測モデル
- survival：生存分析
- caTools, Epi：予測モデルの性能評価に必要なROC曲線が描ける、AUCを計算できる
- XLConnect：エクセルのデータを読み込める、Rオブジェクトをエクセルに保存できる
これらが全部使いこなせれば、データサイエンティストのR技術は充分のように思います。あとは必要になったものをCRANで調べながら分析を進めていくスキルが必要ですね。
【CROSS2013】『活躍するデータサイエンティストの人材像』のセッションメモ
『今日から始まるデータサイエンティスト』という２部構成のセッションのうち、前半の『活躍するデータサイエンティストの人材像』の進行を行なってきました。みなさま第一線で活躍されている方ばかりで、進行をさせて頂くのが恐縮でしたが。。準備も突貫で行ったのですが、みなさまのお陰で無事に実施できました！
進行を行なっていたのでメモを取れなかったのですが、とりいそぎ流れを記録しておきます（あとで覚えている範囲で追記していければ）。
【進行】
倉橋一成（iAnalysis合同会社代表・CAO）
【パネリスト】
草野隆史（株式会社ブレインパッド社長）
益村勝将（トランスコスモス・アナリティクス株式会社COO）
佐々木智之（株式会社gumi執行役員）
西郷彰（株式会社リクルートテクノロジーズ）
【内容】
・オーガナイザーの里洋平氏（株式会社ディー・エヌ・エー）より前座トーク
・各自自己紹介
・分析系の職種を選ばれた理由は？
・データサイエンティストの人物像・スキルセットは？
・データサイエンティストを育てていくためにはどうしたらよいか？
・データサイエンティストを採用するにあたって気を付けることは？
他にも「分析を効率よく進めるためのシステム環境は？」とか「データサイエンティストを企業利益に結び付けるための組織構成は？」とか聞きたかったのですが時間がなかったです。。また機会があれば良いです。
パネリストのみなさん表現方法は違うものの、全体的に人物像は同じようなものを想像されているように感じました。私の言葉で言えば、「学問知識」「ビジネス感」「エンジニアリング」の３軸です（ホームページでも書いています）。この３軸を磨いてスキルの”面積”を大きくしつつ、企画力・提案力を身につけていくことが大事ですね。
あと印象的だったのは、いろんな意味で「バランス力」を持っている人が大事ということ。偏っていると効果の高い分析は難しいのかもしれません。あくまでも本セッションは、マネジメントの視点からなので、そういった意見が多かったのだと思います。もちろんエンジニアリングや数理モデリングに特化した方も、非常に重要です。
また思い出したら追記していきたいと思います。
【CodeIQ】Rで解くデータサイエンティスト問題の解説（R Advent Calendar2012）
先日より、リクルート様のITエンジニアのための実務スキル評価サービス「CodeIQ」で、データサイエンティストに関する問題を出題させて頂いております（問題集はこちら）。先日12/12のおしゃスタ＠リクルートでも少し解説しましたが、Rでの解答例をお見せする時間がなかったので、この機会にブログで公開します（おしゃスタに関するCodeIQ様のブログはこちら）。去年に引き続き勢いだけで参加したR Advent Calendar 2012でしたが、ちゃんとネタが見つかって良かった！！！でも無計画に参加したらクリスマスイブの日に当たってしまったので、、、日付が変わるくらいにさっさと書いてしまいたいと思います！！！爆
【データサイエンティスト初級問題】
【前提】
とある転職サイトから、「とりあえずデータがあるんだけど、、、」と言われてデータを受け取りました。先方は何をして欲しいかまだはっきりと決まってない様子。
受け取ったデータは、
「応募した人の属性データ(oubo_zokusei.csv)」
「応募した時間の記録データ(oubo_kiroku.csv)」
の2種のCSVファイルです。
【課題】
幸いにもデータサイズはあまり大きくなく、Rで読み込める程度。まずはRで読み込んでどんなデータか把握してください。
解答はテキストファイル(.txt)で、
1. Rコード
2. データの概要：どんなデータかまとめた文章
3. 分析提案：もし応募数の予測モデルを作るとしたらあった方がよさそうな変数とその理由
の3つを、この順番で書いて下さい。
データサイエンティストによくある状況なんですが、「データがあるからとりあえず何かやってよ」と言われているパターンです。「仮説がはっきりしてないと分析できませんよ」と言いたくなる気持ちを抑えつつも、何か結果を出す必要があります。受け取ったときはデータの定義ファイルなんかもありません。データを読み込んで始めて分かる状況です。解答例としては、例えば以下のようになります。
【解答例】
1. Rコード
# データチェック
setwd("......")
ID <- read.csv(file("oubo_zokusei.csv", encoding="cp932"), as.is=T, header=F)
oubo <- read.csv("oubo_kiroku.csv", as.is=T, header=F)
head(ID)
head(oubo)
dim(ID)
dim(oubo)
names(ID) <- c("ID", "join_time", "age", "gakureki", "shokureki")
names(oubo) <- c("ID", "time")
oubo$month <- months(as.Date(oubo$time))
hist(ID$age)
summary(ID$age)
table(ID$gakureki, useNA="always")
table(ID$shokureki, useNA="always")
table(oubo$month)
# 全期間での応募数
ouboNum <- aggregate(oubo$ID, list(oubo$ID), length)
names(ouboNum) <- c("ID", "oubo_num")
# 各月の応募数
ouboNumMonth <- aggregate(oubo$ID, list(oubo$ID, oubo$month), length)
names(ouboNumMonth) <- c("ID", "month", "oubo_num")
ouboNumMonth$month <- ifelse(ouboNumMonth$month == "5月", "mon5",
ifelse(ouboNumMonth$month == "6月", "mon6", "mon7"))
library(reshape)
ouboNumMonth1 <- cast(ouboNumMonth, ID~month, value="oubo_num")
ouboNumMonth1 <- as.data.frame(ouboNumMonth1)
ouboNumMonth1[is.na(ouboNumMonth1)] <- 0
# 入会データとマージ
ouboID <- merge(ID, ouboNum, by="ID", all.x=T)
ouboID <- merge(ouboID, ouboNumMonth1, by="ID", all.x=T)
ouboID[is.na(ouboID)] <- 0
# 可視化
cor(ouboID[, 7:9])
plot(ouboID[, 7:9])
# モデル
ouboLM <- lm(mon7~age+gakureki+shokureki+mon5+mon6, data=ouboID)
summary(ouboLM)
cor(predict(ouboLM), ouboID$mon7)
plot(predict(ouboLM), ouboID$mon7, xlim=c(0, 100), ylim=c(0, 100))
abline(a=0, b=1)
library(randomForest)
ouboRF <- randomForest(mon7~mon5+mon6, data=ouboID)
summary(ouboRF)
cor(predict(ouboRF), ouboID$mon7)
plot(predict(ouboRF), ouboID$mon7, xlim=c(0, 100), ylim=c(0, 100))
abline(a=0, b=1)
2. データの概要レポート
zokuseiデータ：恐らくID・登録日・年齢・学歴・勤務状況が記録されている
ouboデータ：恐らく誰が・いつ応募したのか記録されている
20代〜30代、大卒が中心に利用。離職中が在職中の約半分。
とりいそぎ過去の応募数から次月の応募数を予測するモデルを作ってみたところ、単純な線形回帰モデルでR二乗が約0.33であった。ランダムフォレストも同じくらいの精度であった。次月の応募数のうち約３割が説明可能という解釈ができるが、プロットを確認すると応募数が多い外れ値に引きずられているようにも見える。
3. 分析提案
詳細な分析を行うためには、属性データの詳細が必要。またサイトのアクセス履歴も分析することで、ユーザーの行動特性が見つかるかもしれない（ただ、アクセス解析はRではなく専用ソフトが向いているかも）。
将来の応募数を予測することで、営業支援や売上予測につなけることも可能である。そのためには現在のモデルより精度を上げなければならないため、新たな変数が必要であるし、モデルも吟味する必要がある。
【解説】
だいたいここまで分析時間と文章に起こす時間で、１〜２時間といったところです（エラーが出たりしたらもっとかかりますが）。むしろここまでは割とさくっとできるのですが、ここから先が大変です。仮説をヒアリングしたり、この結果を丁寧に説明したり。丁寧に説明しなくても分かる人が増えてくれれば良いんですがね。マッキンゼーが「データサイエンティストをマネジメントする人材が150万人足りなくなる」って言っているのはこのことだと感じてます（詳しくはこの記事など）。あと実務でこういう事をやろうとすると、データを作るまでに結構な時間がかかってしまったりします。分析に慣れている方がDBを操作していると、簡単に作ってくれたりすることもあります。
しかし、出題した後に、結構レベルが高かったかな〜と思ったり。。まだみなさんの回答を見れてないのですが、むしろこんな出題の仕方で回答してもらえて恐縮でした。。。初級の次に中級・上級を予定してたのですが、それは断念しました。
その反省点を活かして、次の問題からは課題をはっきりとさせました。データサイエンティストで必要な学問知識のうち重要な「統計学」と「機械学習」に分けて出題しています。Codingして結果を出せるように、Rを使って回答してもらう、という形式にしています。以下問題と、私の解答例を紹介します。
データサイエンティスト〜統計学編1〜
【問題】
問1. Rを使い、5人分の身長に関し、下記の6つの統計量を求めてください。
# 5人分の身長データ
height1 <- c(168, 173, 152, 181, 175)
(1) 5人分の身長データの平均値
(2) 5人分の身長データの中央値
(3) 5人分の身長データの標本分散
(4) 5人分の身長データの不偏分散
(5) 5人分の身長データの標準偏差（不偏分散を使う）
(6) 5人分の身長データの標準誤差（不偏分散を使う）
問2. 下記のrnorm関数を使うと5人分の身長のデータを擬似生成できます。擬似生成データを使った以下の質問に答えてください。
# 5人分の身長を擬似生成するコード
height2 <- 170 + 10*rnorm(5)
(1) 擬似生成した5人分の身長データを使って
- 標本分散
- 不偏分散
- 標準偏差（不偏分散を使う）
- 標準誤差（不偏分散を使う）
の4つの統計量を計算してください。
計算に使ったRのコードと計算結果を提出してください。(2) 「100人分の身長データ」、「1,000人分の身長データ」を擬似生成し、上記4つの統計量（標本分散、不偏分散、標準偏差、標準誤差）を計算してみてください。人数（データ件数）が増えたとき、この4つの関係性はどうなるか議論してください。
【解答例】
# --1.
height1 <- c(168, 173, 152, 181, 175)
mean(height1)
median(height1)
sum( (height1 - mean(height1)) ** 2)/length(height1)
var(height1)
sd(height1)
sd(height1)/sqrt(length(height1))
# --2.
stat <- function(vec){
a <- sum( (vec - mean(vec))**2)/length(vec)
b <- var(vec)
c <- sd(vec)
d <- sd(vec)/sqrt(length(vec))
c(a, b, c, d)
}
height2 <- 170 + 10*rnorm(5)
stat(height2)
height3 <- 170 + 10*rnorm(100)
stat(height3)
height4 <- 170 + 10*rnorm(1000)
stat(height4)
# 標本分散と不偏分散はほとんど等しくなっていく
# 標準偏差は10付近をばらついている
# 標準誤差はゼロに近づいていく
データサイエンティスト〜統計学編2〜
【問題】
問1. DataScience_stat2.csvのデータは100万人分の身長を模擬的に作成したものです。
この100万人分の身長データから、1000人分のデータをRを使ってサンプリングしてください。
問2. 問1でランダムサンプリングしたデータを使って以下の値を求めてください。
2-1. 平均値
2-2. 標準誤差（不偏分散を使う）
2-3. 平均値の95%信頼区間
問3. 問2-3の信頼区間は、もとの100万人の集団の平均身長（真値）を推測するものです。サンプリングを無限に繰り返すと、理論的には95%の割合で信頼区間の範囲に真値が含まれます。ランダムサンプリングを10000回繰り返し、この事を確認してください。計算に使ったRのコードと簡単な説明を提出してください。
※DataScience_stat2.csvは以下のRプログラムによって作成しています。
set.seed(1)
height <- 170 + 10*rnorm(1000000)
【解答例】
# --1. height <- read.csv("DataScience_stat2.csv", as.is=T) mean(height$x) sd(height$x) sd(height$x)/sqrt(length(height$x)) # --2. set.seed(1) height_sample <- sample(height$x, 1000) MEAN <- mean(height_sample) SE <- sd(height_sample)/sqrt(length(height_sample)) LCL <- MEAN - 1.96*SE UCL <- MEAN + 1.96*SE # --3. CI <- data.frame(LCL=NULL, UCL=NULL) for(i in 1:10000){ set.seed(i) height_sample <- sample(height$x, 1000) MEAN <- mean(height_sample) SE <- sd(height_sample)/sqrt(length(height_sample)) LCL <- MEAN - 1.96*SE UCL <- MEAN + 1.96*SE CI <- rbind(CI, c(LCL, UCL)) print(i) } names(CI) <- c("LCL", "UCL") # ----真値が信頼区間に含まれている割合 CI$BetweenCI <- ifelse(CI$LCL <= mean(height$x) & mean(height$x) <= CI$UCL, 1, 0) table(CI$BetweenCI) prop.table(table(CI$BetweenCI)) # ----グラフでチェックする plot(c(168, 172), c(1, 100), type="n") for(i in 1:100){ lines(c(CI$LCL[i], CI$UCL[i]), c(i, i), col=2-CI$BetweenCI[i]) } abline(v=mean(height$x), lty=2)
データサイエンティスト〜機械学習編1〜
【問題】
問1. Rを使い、DataScience_ML1.csvを読み込み以下の線形回帰モデルを作成してください。
y=x1+x2
問2. Rを使い、問1で作ったモデルに対して以下の回帰診断を行なってください。
(1) ローデータの散布図
(2) 調整済みR二乗
(3) 残差プロット
(4) キャリブレーションプロット
※DataScience_ML1.csvは以下のRプログラムによって作成しています。
set.seed(1)
x1 <- rnorm(10000)
set.seed(2)
x2 <- rnorm(10000)
set.seed(3)
y <- 2*x1 + x2**2 + rnorm(10000)
Data <- data.frame(x1 = x1, x2 = x2, y = y)
【解答例】
# --1. Data <- read.csv(DataScience_ML1.csv, as.is=T) LM1 <- lm(y ~ x1 + x2, data=Data) # --2. plot(Data) summary(LM1) plot(LM1) plot(predict(LM1), Data$y) # --3. Data2 <- data.frame(Data, x3 = Data$x2**2) LM2 <- lm(y ~ x1 + x2 + x3, data=Data2) summary(LM2) plot(predict(LM2), Data2$y)
データサイエンティスト〜機械学習編2〜
【問題】
問1. Rを使い、DataScience_ML1.csvを読み込み、以下のモデルを作成してください。結果変数はyとします。・線形カーネルのSVRモデル
・3次多項式カーネルのSVRモデル
・ガウシアンカーネルのSVRモデル（radial basis）
・シグモイドカーネルのSVRモデル
問2. 作ったモデルに対して予測診断を行うために、各モデルのキャリブレーションプロットとR二乗をチェックしてください。
【解答例】
# --1. library(e1071) Data <- read.csv(DataScience_ML1.csv, as.is=T) SVR1 <- svm(y~., data=Data, type="eps-regression", kernel="linear") SVR2 <- svm(y~., data=Data, type="eps-regression", kernel="polynomial", degree=3) SVR3 <- svm(y~., data=Data, type="eps-regression", kernel="radial") SVR4 <- svm(y~., data=Data, type="eps-regression", kernel="sigmoid") # --2. plot(predict(SVR1), Data$y) plot(predict(SVR2), Data$y) plot(predict(SVR3), Data$y) plot(predict(SVR4), Data$y) cor(predict(SVR1), Data$y) ** 2 cor(predict(SVR2), Data$y) ** 2 cor(predict(SVR3), Data$y) ** 2 cor(predict(SVR4), Data$y) ** 2
以上が解答例になります。解説はまた今度の機会に行おうと思っています。CodeIQのサイトからデータのダウンロードはできなくなっていますが、統計学編と機械学習編は各自の環境でデータを作成できます。ぜひやってみて下さい。今回紹介した5問は分析では入り口ですが、これらができるだけでも実務ではとても役に立ちます。これくらい出来ていれば、他の手法も理解しやすくなると思います。おしゃスタでの発表資料でも問題の解説を少ししていますので、補助的に見て下さい。
今後も分析力を身に着けてレベルアップしていけるように、問題をステップアップできればと思っています。最終的には私の修士論文とか博士論文のネタに辿り着ければな〜と思ったり（ここに置いてます）。
それでは、今回の記事はこのへんで失礼致します。
みなさまよいクリスマスと年末年始をお過ごし下さい！！！！！！！！！！
