Rのspdepのパッケージを使って、コンビニのドミナント出店戦略を定量化する
こんにちは、集計野郎マクガイバーです。
コンビニの出店戦略ではセブンイレブンのドミナント戦略が有名で公式ホームページにも記載してあるくらいです。
しかし、その実態に関してはどの程度そうであるのかといったような厳密な検証情報はなく、特定都市における出店数を比較したようなものや、単に四国にセブンイレブンが出店していないという情報を根拠にしたものが多数見られます。
関西ではローソンがドミナントしているのはいうまでもないし、特定の地域を抜き出して一般的な結論を持ち出すのは問題に思われます。
幸いなことに去年集めたものですが、データが手元にあるので、今日はこのドミナント出店戦略を数値化、可視化してみます。
（以後、ですます調ではなくなります。）
ロウデータで見る各コンビニの出店戦略
ここでは、ロウデータを地図上にプロットして各コンビニチェーンの出店戦略がどのようになっているかを見てみる。
赤：セブンイレブン、青：ローソン、緑：ファミリーマート、オレンジ：サークルK、紫：その他
以下の図は東京都内の中で各チェーンの出店戦略を最も顕著に表していると考えられる箇所を抜粋したもである。
青のローソンが多く見えるが、これは色の見えやすさも要因であることに注意する。
ここから見てとれる出店戦略の仮説は以下の通りである。
セブンイレブン：東京の中心部からやや離れた地域に店舗が多く、電車の沿線ではないところにも出店している。
サークルK：出店地域は狭く傾向は読み取りにくいが、セブンイレブンと同様の出店傾向が見られる。
ローソン：都心を中心に幅広く出店している。沿線上の出店傾向が強い。
ファミリーマート：都心を重点的に出店している。
こうしてみてみると沿線上の出店傾向があるローソンもセブンイレブンと変わらないドミナント戦略のように見える。
ドミナンスを可視化する
続いて、各コンビニチェーンのドミナンスを可視化する。
ここでは、日本地図においてコンビニが存在している領域を切り出したバウンディングボックスを
一辺の最大セル数が300になるように正六角形で区切り、当該セル内における最も店舗数の多いコンビニチェーンを当該チェーンの支配している地域として色を塗りわける。
これを支配地域と呼ぶ。ただし、同一セル内で最も多い出店数が複数チェーンでタイであった場合はランダムにその中から選ぶものとする。
上記の図は後出しではあるが、各コンビニチェーンの出店戦略を最も表していると考えられる地図である。
中国地方から東海地方までを切り出している。色は下記の通りである。
赤：セブンイレブン、青：ローソン、緑：ファミリーマート、オレンジ：サークルK、紫：その他
各チェーンの特徴をまとめると以下のように読み解ける。
セブンイレブン：典型的なドミナント戦略、中国地方や滋賀などに重点的に出店しているが、他の地方には見られない。
サークルK：東海地方に重点的に出店しており、ほとんどの領域を支配している。
ローソン：足元である関西地方の店舗を支配している他に、比較的人口が少ないと思われる山間部や地方にも幅広く出店している。
ファミリーマート：ドミナントしている地域は少ないが、淡路島やその他人口が少ないと思われる地方に所々見られるように幅広く出店している。
またこの地図にはないが、東京の都心部、山手線圏内はファミリーマートが支配している（4セル）。
ちなみに日本全国で見ると以下のような感じになる。
ドミナント出店を定量化する。
可視化では仮説のレベルを出ないため、ドミナント出店をしているかしていないかを表す指標を作成する。
定量化するにあたって、その統計量として以下の3つの条件が必要である。
- ドミナント性・・・ある特定の地域を支配していること。これは先ほどの六角形グリッドのデータを用いることで満たす事が出来る。
- 空間的自己相関性・・・支配地域が空間的に隣接している数が多いこと。支配地域が隣接していること。
- 支配地域の数に対する非依存性・・・隣接数は支配地域が多いと高くなるので、店舗数の少ないチェーンには不利であるため。
この条件を満たす空間的自己相関の統計量として、今回はカテゴリカルな空間的自己相関のjoint-countの検定に用いるZスコアが使えそうなので、これを用いることにした。
（これ以外に上記の条件を満たし、ドミナント戦略を定量化できるよりよい指標は思い浮かばなかった。）
このスコアは帰無仮説（各チェーンにおける支配地域数を変えずにランダムにセルを配置した場合）におけるチェーンの支配地域の隣接数と実際のデータの隣接数を
標準偏差を踏まえた形で比較するものである。Rのspdepのjointcount.multiを用いることで算出できる。
結果
結果は以下の通りである。
Zスコアでは、サークルK・サンクスが最も高い数値を示している事がわかった。
実際にはサークルK・サンクスもセブンイレブンと変わらない程度の高いドミナント戦略をとっている事がわかった。
また、これと比べるとファミリーマートやローソンは幅広く出店する戦略を採択している。
関西を足元に持つローソンがドミナンス指標ではファミリーマートよりやや優位にでる。
また、店舗数比率/支配地域比率のインデックス（より少ない店舗で地域を支配しているという指標）と空間的自己相関性を表すZスコアが比例する形になった。
結論
セブンイレブンとサークルK・サンクスが最も高いドミナント戦略をしていたチェーンに群となる。
その後、ローソン、ファミリーマートと続いた。
最後に
地図はQJIS使って作っています。ドロネー分割したり、人口や沿線とその乗降客数の情報等も使うと他のチェーンの戦略がわかって面白そう。
ちなみに今回この分析をしているのは、コンビニの出店戦略に興味があったからではなく、別の技術的な検証を行うためでした。
それについてはまた時間があれば、ブログを更新したいと思います。
追記：すいません、結論内部で間違った記述をしていた箇所があるため、本文を修正しました。
arulesSequenceをWindowsRStudioから呼び出せるように修正するパッチ
わんばんこ、年末も近づいてきてアソシエーション真っ盛りの季節ですね！
さて、RのarulesSequenceパッケージのcspade関数をWindows7でRStudioから呼び出そうとしたところ、エラーになる事象が発生。
同様の事象は海外でも報告されている。（以下のリンクを参照）
http://grokbase.com/t/r/r-help/124gp26w8r/r-system2-and-cspade-on-windows
https://stat.ethz.ch/pipermail/r-help/2012-April/308802.html
特徴的なのは、RStudioから呼び出そうするとエラーになるが、普通のRConsoleから呼び出す時には動くこと。
で、原因調べてみるとcspade内でsystem2コマンドで実行した結果の標準出力の出力先をstdoutオプションで指定できるはずだが、stdout引数で指定された出力先にリダイレクトされていないのが原因の模様。このため、本来Rのテンプディレクトリに出力されているはずの標準出力結果をリダイレクトしたファイルが存在しないと怒られる。
> s<-cspade(zaki, parameter=list(support=0.5), control = list(verbose = TRUE)) parameter specification: support : 0.5 maxsize : 10 maxlen : 10 algorithmic control: bfstype : FALSE verbose : TRUE summary : FALSE preprocessing ... 1 partition(s), 0 MB [0.02s] mining transactions ... NA MB [0.01s] reading sequences ... 以下にエラー file(con, "r") : コネクションを開くことができません 追加情報: 警告メッセージ： In file(con, "r") : ファイル 'C:\Users\******\AppData\Local\Temp\Rtmp2VKdgt\cspade16503e135e49.out' を開くことができません: No such file or directory
でもって、さらに調べてみると、cspadeコードの上部で、この問題に対応できるようにWindowsの場合のsystem2関数を上書きしている。
ただし、GUIがRStudioの場合を考慮していない。このため、RGuiであれば正しく実行できるが、RStudioだと実行できないということになる。
そこで、この問題を解決するためのアドホックパッチを作成した。
下記のコードを実行すれば解決すると思われる。ただし、cspadeではなく、arulesSequencesのcspadeを呼び出すこと。
ポイントは#change modified function's environmentで修正した関数mod.cspadeに対して環境を設定しているところ。
これをしないとcspade内部で.formatSPというエクスポートされていない関数が呼び出せず落ちてしまうので要注意。
> s<-arulesSequences::cspade(zaki, parameter=list(support=0.5), control = list(verbose = TRUE)) parameter specification: 以下にエラー cat(.formatSP(parameter), sep = "\n") : 関数 ".formatSP" を見つけることができませんでした
パッチ
library(arulesSequences)
#modified cspade
mod.cspade<-function (data, parameter = NULL, control = NULL, tmpdir = tempdir())
{
#fix only here
if (.Platform$OS == "windows" && .Platform$GUI %in% c("Rgui", "RStudio"))
system2 <- function(command, args = character(), stdout = "",
...) {
if (is.character(stdout) && nzchar(stdout)) {
args <- c(args, ">", stdout)
stdout <- NULL
}
args <- c("/c", shQuote(command), args)
command <- Sys.getenv("COMSPEC")
if (!nzchar(command))
stop("environment variable 'COMSPEC' not set")
base::system2(command, args = args, stdout = stdout,
...)
}
if (!inherits(data, "transactions"))
stop("'data' not of class transactions")
if (!all(c("sequenceID", "eventID") %in% names(transactionInfo(data))))
stop("slot transactionInfo: missing 'sequenceID' or 'eventID'")
if (!all(dim(data)))
return(new("sequences"))
parameter <- as(parameter, "SPparameter")
control <- as(control, "SPcontrol")
if (control@verbose) {
t1 <- proc.time()
cat("\nparameter specification:\n")
cat(.formatSP(parameter), sep = "\n")
cat("\nalgorithmic control:\n")
cat(.formatSP(control), sep = "\n")
cat("\npreprocessing ...")
}
exe <- "bin"
if (.Platform$r_arch != "")
exe <- file.path(exe, .Platform$r_arch)
exe <- system.file(exe, package = "arulesSequences")
file <- tempfile(pattern = "cspade", tmpdir)
on.exit(unlink(paste(file, "*", sep = ".")))
opt <- ""
nop <- ceiling((dim(data)[1] + 2 * length(data@data@i)) *
.Machine$sizeof.long/4^10/5)
if (length(control@memsize)) {
opt <- paste("-m", control@memsize)
nop <- ceiling(nop * 32/control@memsize)
}
if (length(control@numpart)) {
if (control@numpart < nop)
warning("'numpart' less than recommended")
nop <- control@numpart
}
out <- paste(file, "stdout", sep = ".")
if (FALSE) {
asc <- paste(file, "asc", sep = ".")
write_cspade(data, con = asc)
if (system2(file.path(exe, "makebin"), args = c(asc,
paste(file, "data", sep = "."))) || system2(file.path(exe,
"getconf"), args = c("-i", file, "-o", file), stdout = out))
stop("system invocation failed")
file.append("summary.out", out)
}
else makebin(data, file)
if (system2(file.path(exe, "exttpose"), args = c("-i", file,
"-o", file, "-p", nop, opt, "-l -x -s", parameter@support),
stdout = out))
stop("system invocation failed")
file.append("summary.out", out)
if (length(parameter@maxsize))
opt <- paste(opt, "-Z", parameter@maxsize, collapse = "")
if (length(parameter@maxlen))
opt <- paste(opt, "-z", parameter@maxlen, collapse = "")
if (length(parameter@mingap))
opt <- paste(opt, "-l", parameter@mingap, collapse = "")
if (length(parameter@maxgap))
opt <- paste(opt, "-u", parameter@maxgap, collapse = "")
if (length(parameter@maxwin))
opt <- paste(opt, "-w", parameter@maxwin, collapse = "")
if (!length(control@bfstype) || !control@bfstype)
opt <- paste(opt, "-r", collapse = "")
if (control@verbose) {
t2 <- proc.time()
du <- sum(file.info(list.files(path = dirname(file),
pattern = basename(file), full.names = TRUE))$size)
cat(paste("", nop, "partition(s),", round(du/4^10, digits = 2),
"MB"))
cat(paste(" [", format((t2 - t1)[3], digits = 2, format = "f"),
"s]", sep = ""))
cat("\nmining transactions ...")
}
out <- paste(file, "out", sep = ".")
if (system2(file.path(exe, "spade"), args = c("-i", file,
"-s", parameter@support, opt, "-e", nop, "-o"), stdout = out))
stop("system invocation failed")
if (control@verbose) {
t3 <- proc.time()
du <- file.info(out)$size
cat(paste("", round(du/4^10, digits = 2), "MB"))
cat(paste(" [", format((t3 - t2)[3], digits = 2, format = "f"),
"s]", sep = ""))
cat("\nreading sequences ...")
}
out <- read_spade(con = out, labels = itemLabels(data))
out@info <- c(data = match.call()$data, ntransactions = length(data),
out@info, support = parameter@support)
if (control@verbose) {
t4 <- proc.time()
cat(paste(" [", format((t4 - t3)[3], digits = 2, format = "f"),
"s]", sep = ""))
cat("\n\ntotal elapsed time: ", (t4 - t1)[3], "s\n",
sep = "")
}
if (!control@summary)
unlink("summary.out")
out
}
#change modified function's environment
environment(mod.cspade)<-environment(cspade)
#unlock and assign
unlockBinding("cspade", getNamespace("arulesSequences"))
assign("cspade", mod.cspade, getNamespace("arulesSequences"))
lockBinding("cspade", getNamespace("arulesSequences"))
#doing cspade(note:cpade itself does not change.try to call cspade not arulesSequences::cspade)
data(zaki)
s<-arulesSequences::cspade(zaki, parameter=list(support=0.5), control = list(verbose = TRUE))
1次元ハイブリッドモンテカルロを実装してみた
1次元HMCのコードを書いたのでメモ
動機はNUTSに挑戦する前にパラメータと尤度、ハミルトニアンの振る舞いについて理解しておきたかったので
偉大なる先輩のblogはこちら
ハイブリッドモンテカルロを実装してみた
http://d.hatena.ne.jp/koh_ta/20121116/1353074113
########################################################
#HMC
#initial value
epsilon<-0.25
N<-1000
L<-4
p<-rep(NA, N)
q<-rep(NA, N)
s<-rep(NA, N)
h<-rep(NA, N)
p[1]<-rnorm(0, 1, n=1)
q[1]<-3
h[1]<-hamiltonian(q[1], p[1])
SD<-2
MU<-1
#q obey to N(3, 2.2) f:function itself, df:differentiation function
f<-function(x){dnorm(x, mean=MU, sd=SD)}
df<-function(x){dnorm(x, mean=MU, sd=SD) * (-1/(2*SD^2)*(2*x-2*MU) )}
hamiltonian<-function(q,p){log(exp(f(q))-(p^2)/2)}
#ploting function and differentiation function
plot(df, -1, 3)
par(new=T)
plot(f, -1, 3, col="red")
leapfrog<-function(q_var, p_var, epsilon){
p_var<-p_var-(epsilon*df(q_var)/2)
q_var<-q_var+epsilon*p_temp
p_var<-p_var-(epsilon*df(q_var)/2)
cat(q_var, p_var, "\r\n")
c(q_var, p_var)
}
for(i in 2:N){
p_temp<-p[i-1]
q_temp<-q[i-1]
for(j in 1:L){
qp_temp<-leapfrog(q_temp, p_temp, epsilon)
q_temp<-qp_temp[1]
p_temp<-qp_temp[2]
}
#accept by hamiltonian ratio
alpha<-min(1, hamiltonian(q_temp, p_temp)/hamiltonian(q[1], p[1]), na.rm=TRUE)
if(alpha>runif(1)){
q[i]<-q_temp
p[i]<- (-1)*p_temp
s[i]<-1
}else{
q[i]<-q[i-1]
p[i]<-p[i-1]
s[i]<-0
}
h[i]<-hamiltonian(q[i], p[i])
}
ylim<-c(-2, 4)
plot(p,type="l",ylim=ylim)
par(new=T)
plot(q, col="red",type="l",ylim=ylim)
par(new=T)
plot(s, col="blue", ylim=ylim)
par(new=T)
plot(h, col="green",type="l",ylim=ylim)
PostGISのよく使う機能をまとめた
PostGISの中でも頻繁に使用する機能について調べたので、参照用にSQLとその使用方法についてまとめておきます。
GISは空間参照系や楕円体のあたりなど技術的にもマニアックだが、とりあえず理解している範囲内で基本的なことをメモ。
PostGISについて
PosgreSQLには地理情報に関する機能を集めたPostGIS呼ばれる拡張がある。メリデメは以下の通り。
メリット
- 複雑な集計操作を豊富な関数で処理でき、GIS基盤のミドルウェアとして利用可能。
- データベースに格納できるので既存のデータとの結合やデータの共有が容易。
- 毎回必要に応じて集計すればよいので、粒度や集計単位の異なるshpファイルをたくさん持たなくてもよくなる。
- QGIS等とも相性がよく、手軽にSQLで集計して主題図を記述可能。
デメリット
- 毎回集計するのでメモリに乗るデータなら普通にGISで処理した方が早いかも。
- データベースへのロードやエクスポート作業が必要。
Geometry型
PostGISではGeometry型とよばれるデータ型のカラムを作成し、ジオメトリー操作用の関数(ST_XXXXXX)でこれらを操作し、地理的な検索や集計を行ったり、地図として表示させるポリゴンを取得します。
PostGISには、３つのジオメトリー型が存在する。点（位置、ランドマーク）を表すPOINT、線（道路や川）を表すLINESTRING、面（領域、町丁目）を表すPOLYGONの３つ。厳密にはそれらの複合型、集合型もある。
Geography型
Geometry型だけでなく、Geography型と呼ばれるデータ型もあるが、これは日本国内の地図を触っているだけなら必要なさそうなので無視することにした。
厳密な距離を算出したいような場合でも必要かもしれない。残念ながらGeography型は対応している関数が少ないのでGeometryで持っている方が無難。
必要な時だけGeometry型で持っているデータをGeographyに変換して操作するのがベター。
はまったポイント
PostGISを使う上で初心者である私がはまった注意点は以下の3つです。これから基盤を構築される型は気をつけましょう。
- SRIDを統一する。・・・SRIDとは、2次元座標（数値）を実際の位置と対応させる空間参照系のユニークなIDのことを指します。PostGISでは各ジオメトリー型カラムにSRIDを設定することができる。shp2pgsqlを使ってインポートする際に、SRIDを変換する機能があるのでそれを使ってあらかじめデータベース用の共通SRIDデータのみ含むようにしておくと無難。わからない人はとりあえずSRID＝4326に統一するようにするとよいと思われる。
- 文字コード・・・落としてきたshpファイルが大体Shift-Jisだったりするが、shp2pgsqlでDBへとインポートする際に文字コードを指定できるので、インポートする時点でDBと同じ文字コードに統一しておく。
- Oracleの情報でないか確認する・・・PostGISを触っていると色々と情報ない中、Googleで調べる羽目になるが、"Postgre"をキーワードに入れていても、"Oracle"の情報が入っていたりする。元々SQLが似ているので混同しないように気をつける。実際何度か間違えた。
インストール
インストールは少し古いが、こちらを参照すること。
Stack Builderで適当にやってたら入ると思います。万が一、はまった時は思いっきり悔しがりましょう。
データを作成する
ポリゴンデータのインポート・・・PostGISにはshpファイルをインポートするshp2pgsqlというコマンドがあります。
shp2pgsqlは、shpファイルを引数に取り、SQLを出力します。このSQLを実行することで、データをPostgreSQLに取り込むことが出来ます。
pgAdminでもインポートできるけど、文字コード周りはトラブルの種なので面倒がらずに、こちらを用いたほうが無難。
ここでは国勢調査のデータより、渋谷区の町丁目ポリゴン情報をインポートします。
- shpファイル用テーブルの作成
C:\Users\Test>shp2pgsql -s 4326 -W cp932 -c \Users\Test\Desktop\A00200521201
0DDSWC13113\h22ka13113.shp > createtable.sql
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
createtable.sqlの中身は以下の通り。これで"h22ka13113"にshpファイルと同様のテーブルが作成されます。
SET CLIENT_ENCODING TO UTF8; SET STANDARD_CONFORMING_STRINGS TO ON; BEGIN; CREATE TABLE "h22ka13113" (gid serial, "area" numeric, "perimeter" numeric, "h22ka07_" int4, "h22ka07_id" int4, "ken" varchar(2), "city" varchar(3), ... INSERT INTO "h22ka13113" ("area","perimeter","h22ka07_","h22ka07_id","ken","city",....B37EBD24140C10B47D5B2766140908F9FB0EBD24140'); COMMIT;
- 列の作成・・・テーブルschema.table_stに地理情報データ含むジオメトリ型の列を作成します。
ちなみに、今回のように一つのテーブルに複数の異なるタイプの地理情報型を持たせるのはあまり美しくないのでやめておきましょう。
（今回はテスト用のデータのために、複数のテーブルを作りたくなかったので一つにまとめた。）
SELECT AddGeometryColumn('schema','table_st','column_name_geom1',4326, 'POINT', 2); SELECT AddGeometryColumn('schema','table_st','column_name_geom2',4326, 'LINESTRING', 2); SELECT AddGeometryColumn('schema','table_st','column_name_geom2',4326, 'POLYGON', 2);
ポリゴンデータの作成・・・ポリゴンデータを作成します。
- 点の作成:渋谷駅
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (123, '渋谷駅', ST_GeomFromText('POINT(139.70066 35.65878)', 4326), NULL, NULL);
- 線の作成:LINESSTRING1
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (234, 'LINE-1', NULL,ST_GeomFromText('LINESTRING(139.66657 35.67380, 139.67651 35.67690, 139.68333 35.67981)',4326),NULL);
- 面の作成:POLYGON-1
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (345, 'PLOYGON-1', NULL, NULL, ST_GeomFromText('POLYGON((139.69962 35.67493,139.69615 35.65587,139.70990 35.65822,139.69962 35.67493))',4326));
- 索引の作成・・・ジオメトリーカラムに索引を作成します。索引を使うにはSQLの書き方で工夫することが必要だが、念のため作成しておく。
CREATE INDEX index_geom1 ON public.table_st USING gist (column_name_geom1);
- SRIDの設定・・・列の作成時に全て設定しているためここでは特に必要ないが、列のSRIDの設定に関するSQLは下記の通り。
SELECT UpdateGeometrySRID('schema','table_st','column_name_geom1',4326);
＊類似の関数でST_SetSRIDというのがあるが、これはテーブル定義で持っているSRIDの値を変えるだけで座標系を直接変更するものではないらしいので注意。
- ジオメトリー型のSRIDを変更した情報を取得する。
SELECT UpdateGeometrySRID('schema','table_st','column_name_geom1',4326);
上記のSQLで、作成した結果は以下の通りです。
ポリゴン情報の取得
- 重心の計測・・・ポリゴンから重心を取得します。
SELECT ST_Centroid(column_name_geom2) FROM table_st WHERE id = 234;
- 距離の計測・・・点と点の距離を測る。
SELECT ST_distance(yoyogi.column_name_geom1,shibuya.column_name_geom1) dist, ST_distance(yoyogi.column_name_geom1,shibuya.column_name_geom1, false) meter FROM table_st shibuya, table_st yoyogi WHERE shibuya.name = '渋谷駅' AND yoyogi.name = '代々木八幡駅';
| dist | meter |
| 0.016518376262419 | 1658.68010562005 |
- 面と面の距離を測る。
SELECT ST_distance(p2.column_name_geom3,p1.column_name_geom3) dist, ST_distance(p2.column_name_geom3,p1.column_name_geom3, false) meter FROM table_st p1, table_st p2 WHERE p1.name = 'POLYGON-1' AND p2.name = 'POLYGON-2';
| dist | meter |
| 0.00402964942162027 | 366.006648511634 |
＊値の結果より、面と面の距離は最短距離を計測しているようである。
- 線と面の距離を測る。
SELECT ST_distance(l1.column_name_geom2,p1.column_name_geom3) dist, ST_distance(l1.column_name_geom2,p1.column_name_geom3, false) meter FROM table_st p1, table_st l1 WHERE p1.name = 'POLYGON-1' and l1.name = 'LINE-1';
| dist | meter |
| 0.0170052491895865 | 1568.26832718958 |
＊値の結果より、線と面の距離は最短距離を計測しているようである。
- 単一の線の距離を測る。
SELECT ST_length(column_name_geom2) length, ST_length(column_name_geom2, false) meter FROM table_st WHERE id = 234;
| length | meter |
| 0.0178270688501295 | 1657.57701436593 |
- X、Yの座標値を取得する。
SELECT ST_X(column_name_geom1) X, ST_Y(column_name_geom1) Y FROM table_st WHERE name = '渋谷駅'
| x | y |
| 139.70066 | 35.65878 |
検索と集計
- 共有・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにジオメトリーが重なっているかどうかに基づいて検索をかけます。
POLYGON-1と重複する渋谷区の町丁目ポリゴンを取り出すのが下記のSQLです。ST_Intersectsを使用します。
SELECT moji FROM h22ka13113,table_st WHERE ST_Intersects(column_name_geom3, geom) AND table_st.name = 'POLYGON-1';
上記のSQLの結果であるPOLYGON-1と渋谷区の共有部分のある地域は以下のとおりです。
- 包含・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにジオメトリー間の包含関係に基づき検索を行います。
ここでは渋谷駅が含まれるポリゴンを検索します。ST_Withinを使用します。
geomにcolumn_name_geom1が含まれています。引数の順番に注意してください。
SELECT moji FROM h22ka13113,table_st WHERE ST_Within(column_name_geom1,geom) AND table_st.name = '渋谷駅';
- バッファ検索・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにバッファ検索を行います。
今度は渋谷駅から半径500m以内のバッファに含まれる町丁目を検索します。ST_DWithinを使用します。
column_name_geom1から500m以内の範囲にgeomが含まれています。引数の順番に注意してください。
＊何故か私の環境だと最後のuse_spheroidのbooleanに値を入れないとメートル単位の検索にならず、これはtrueでもfalseでも同様の結果です。
SELECT moji FROM h22ka13113,table_st WHERE ST_DWithin(column_name_geom1,geom, 500, true) AND table_st.name = '渋谷駅';
以下の図は上記の二つのSQLの結果です。
濃い紫が「道玄坂1丁目」で渋谷駅を含む包含関係のSQLの結果です。
ピンクが渋谷駅から500m以内のバッファ検索のSQLの結果です。
- ポリゴンの結合・・・ジオメトリーを結合します。市区町村のポリゴンを都道府県に変更したりする場合などに用います。
GROUP BYの単位で結合されます。下記のSQLで渋谷区全体のポリゴンが返されます。
SELECT st_union(geom) FROM h22ka13113;
- ポリゴンの簡素化・・・町丁目レベルで保有しているような大量のポリゴンを結合すると、一つのポリゴンのサイズが大きくなりすぎて複雑すぎる形になるので簡素化する必要があります。
ポリゴンの簡素化にはST_Simplify、ST_SimplifyPreserveTopologyを用います。前者だと実際に扱うような汚いポリゴンだと、不正なポリゴンが返ってくることがあったので
後者をお勧めします。なお、アルゴリズムにはDouglas–Peucker algorithmが使用されており、簡素化を行う閾値となるパラメータεを指定する必要があります。
パラメータεは大きければ大きいほど簡素化されます。下の図は渋谷区のポリゴンをST_Unionで合成した後に、ST_Simplifyでパラメータを変えて簡素化した際の例です。
SELECT st_union(geom) FROM h22ka13113;#赤：簡素化なし SELECT st_simplify(st_union(geom), 0.0008) FROM h22ka13113;#緑：ε=0.0008 SELECT st_simplify(st_union(geom), 0.002) FROM h22ka13113;#青：ε=0.002
管理
- PostGISのバージョンを問い合わせる
SELECT postgis_full_version();
- 地理情報カラムの情報を問い合わせる
SELECT * FROM geometry_columns;
| postgres | public | table_st | column_name_geom1 | 2 | 4326 | POINT |
| postgres | public | table_st | column_name_geom2 | 2 | 4326 | LINESTRING |
| postgres | public | table_st | column_name_geom3 | 2 | 4326 | POLYGON |
| postgres | public | h22ka13113 | geom | 2 | 4326 | MULTIPOLYGON |
- 空間参照系の情報を問い合わせる・・・ここにIDの存在するSRIDしかPostGISは認識できないので注意、必要なものがなければ追加する。
SELECT * FROM spatial_ref_sys;
| 4324 | EPSG | 4324 | GEOGCS["WGS 72BE",DATUM[... | +proj=longlat +ellps=WGS... |
| 4326 | EPSG | 4326 | GEOGCS["WGS 84",DATUM["]... | +proj=longlat +datum=WGS... |
| 4463 | EPSG | 4463 | GEOGCS["RGSPM06",DATUM["... | +proj=longlat +ellps=GRS... |
| 4470 | EPSG | 4470 | GEOGCS["RGM04",DATUM["Re... | +proj=longlat +ellps=GRS... |
DBからRへとポリゴンデータをダウンロードする
最後に、Rで主題図を書いたりする時に必要になる、DB上に保存されているジオメトリーデータをRに読み込む方法について解説します。
やり方は少し厄介ですが、以下の方法で実現します。
- RからRPostgreSQLを使ってデフォルトドライバーでPostgreSQLにアクセスします。
- PostGISのST_AsText関数でポリゴンを地理情報を文字列で表すWKTに変換してSQLを実行しRにcharacterとして読み込みます。
- これをRのrgeos::readWKTでspクラスに変換します。
このとき、ODBCなんかを使うと文字列長制限で大きなジオメトリーデータはダウンロードできないため、RPostgreSQLを使用しています。
とりあえずもっと賢い方法はあるんでしょうけど、今のところこれで特に不便もしていないのでこれを使用しています。
下記はPostgreSQLに入ってる渋谷区のデータをダウンロードして、人口に関する主題図を作成するRコードです。
library(DBI)
library(RPostgreSQL)
library(rgeos)
#PostgreSQLのドライバーを使用して、DB接続を行う。
conn <- dbConnect(dbDriver("PostgreSQL"),user="postgres",host="localhost", password="password",dbname="dbname")
#文字コードを変更する。環境によっては不要
postgresqlpqExec(conn, "SET CLIENT_ENCODING TO 'SJIS'")
#手もとの環境では、2億文字までは文字列で取得可能であった。20億文字でメモリ不足で落ちた。
#使用する際はあまり文字が大きくなり過ぎないようにST_Simplifyなどを使うように注意する。
nchar(fetch(dbSendQuery(conn, "SELECT repeat('yeah!', 200000000);"), n=-1))
#渋谷区のポリゴン情報をテキストとして抜き出す。
df<-fetch(dbSendQuery(conn, "SELECT ST_AsText(geom) geom, jinko from h22ka13113"), n=-1)
plot(readWKT(text=df$geom[1], id=rownames(df)[1]))
#rgeos::readWKTを使ってDBから取得しcharacter型のWKTをRのspクラスに変換する。
spps<-SpatialPolygons(sapply(1:nrow(df), function(i){readWKT(df[i, "geom"], id=rownames(df)[i])@polygons}))
spdf<-SpatialPolygonsDataFrame(Sr=spps, data=df)
#1000人ずつで7色に塗り分ける。
plot(spdf,col=heat.colors(7)[as.integer(cut(spdf$jinko, breaks=0:7*1000))])
title("渋谷区人口")
legend(x=spps@bbox[1, 1]-0.00005, y=spps@bbox[2, 1]+0.02, title="人口",
legend=levels(cut(spdf$jinko, breaks=0:7*1000)), fill=heat.colors(7), bg="white", cex=0.8)
enjoy!
お役立ちリンク
| タイトル | URL | 摘要 |
| PostGIS 2.0.0マニュアル日本語訳 | http://www.finds.jp/docs/pgisman/2.0.0/postgis.html | マニュアルです |
| PostGIS入門 | http://workshops.boundlessgeo.com/postgis-intro-jp/index.html | PostGIS使い方に関する包括的なチュートリアル、一度は目を通すべき |
| GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ | http://d.hatena.ne.jp/tmizu23/20091215/1260868350 | 測地系、EPSGコードに関するよい解説 |
| PostGISとは？ | http://lets.postgresql.jp/documents/tutorial/PostGIS/1/ | PostGISに関する基礎的なチュートリアル、少し情報は古いが、インストール周りの解説が丁寧で読みやすい |
第1回BUGS/Stan勉強会を開催しました
こんにちはStanエヴァンジェリストの駆け出しベイジアンです。掲題の通り、第1回BUGS/stan勉強会を開催したので報告します。
BUGS/Stanって何？
改めて、BUGS/Stanとは何ぞやということを説明しておくとベイズ推定による柔軟なモデリングを実現するためのDSL言語です。
これでモデルを記述すると、非エンジニアなデータ解析屋のみなさんでも（比較的）簡単にデータ分析が出来ます！
このあたりの技術は日本語で枯れていないのでこの勉強会を通して色々と技術を枯らして行きたいと思います。
特にベイズモデリングは恣意性が高く、細かいテクニックと経験が必要なのでそういった悩みのあたりも共有できたらと思っています。
会場はドリコムさんの100人くらいは入る広い会場をお借りしました。ドリコムさんに感謝！！
参加者は全員で10名程度でした。かなり高学歴な方ばかりが集まっていらっしゃったので大変恐縮でしたが、私からの発表でした。
1.「Introduction Of Stan（発表者 @TeitoNakagawa）」
Stanの基本的な使い方について説明するものでした。スライドが英語なのはとある外国の方からのお願いです。
発表内容はStanを使いたい動機から始まってインストール、文法、実行方法について解説しています。
技術的なことは資料を読んで頂くとして以下に要点だけまとめておきます。
BUGSやStanを用いてベイズ推定を行うモチベーション
- 仕事でデータ解析をしているが、お客様から提供されるデータや国勢調査のような統計データはデータが少なく欠損も多く、売上のようなデータは偏りが大きく個体差を検討せずにモデリングできない。
- 一方で、コンサルという仕事柄、機械学習のような方法を用いることは出来ずデータ間の関係性に関して詳細な言及を要求されるので、表現力豊かなモデリングを行いたい。
その中でもStanには以下のようなメリットがあるので私はStanがお勧めです。
Stanのメリット
- StanはBUGSと比べてもインストールも簡単。
- Rと親和性が高く、Rから簡単に使える。
- C++にコンパイルされて実行ファイルに変更されるので実行スピードが速い（ただし、コンパイル時間はかかる）
- HMCとNUTSにより収束も早い。
- データ型が多く、プログラミングが容易（人によって分かれるとは思います。）
2.「ベイズ統計学とサンプラー(発表者 @kos59125)」
モデルとは何か、ベイズ推定とは、そしてMCMCでどうしてベイズ推定が出来て、どのような手法でサンプリングできるのかGibbsサンプラー、HMCのさわりとDSL言語についてというこの分野の入門的な内容を全て網羅して頂きました。とってもスライドがわかりやすく勉強になりました。
印象に残っているのは、われわれはデータを入手できるがそのデータの出自であるモデルとは「神のみぞ知る世界」であり、「モデルは真実ではない」という言葉です。
しかしながら、我々としてはなんとしてもそのモデルに近づきたいので限られたデータから、ランダムサンプリングでそのモデルを推定する努力をせざるを得ないわけです。
ちなみにDSL言語とは、特定の領域に特化した言語のことでありBUGSやStanはこのうちベイズモデリングに特化した言語ということになります。
これらのモデリングに関して自らC++などでコードを書かれている方も多いようですが、そういったコードの書けない初心者はBUGS/Stanでお手軽にモデリングしちゃいましょう。
3.「StanでLotka-Vlotraモデルのパラメータ推定を試してみたけどうまくいかなかった話(発表者 @kos59125)」
Stanを使って、生態学で被食者捕食者のモデリングに用いられるロトカ・ボルテラ方程式の解を求めるというものでした。人工的に生成したデータにランダムノイズを加え、Stanでモデリングして微分方程式を解くというのをやってみたけど、うまく収束しなかったので何か問題あるでしょうかという趣旨の発表でしたが@berobero11さんから「加えているノイズが大きすぎるのではないか」とのご回答がありました。懇親会で微分方程式を解く場合はBUGSやStanを使うより、普通にプログラミングをしたほうが良いかもとの話もありました。
4．「BUGSを使うメリット (発表者 @berobero11)」
BUGSを使うメリット from berobero11
BUGS界のハルク・ホーガンこと@berobero11さんの発表です。大人のビジネスマンである@berobero11さん曰く、まずはメリットデメリットを語るのが社会人の基本とのことです。勉強になります。ちなみにBUGSのメリットは以下の通り、
メリット
- ドメイン知識がフル活用できる。
- 書くコードが少なくてすむ。
- 結果がRからとても使いやすい。
- パラメータの同時分布が得られる。
デメリット
- 時間がかかる
- モデルの書き方がわかりにくい
- エラーメッセージが意味不明で殺意がわく
だそうです。
BUGS/Stan/JAGSのうち、皆さんがどれをやるべきかのチャートが最後のほうにあるので皆さんも是非やってみてください。
ちなみに私はJAGSになりました。
が、質疑応答の中でJAGSの話になった時、「JAGSとStanは現時点で収束の早さに関して同程度のパフォーマンスのようですが、おそらくStan1.4.0でStanが圧勝するだろう」とのことだったのでしばらくはStanを使い続けようと思います。
５．「TrueSkillモデルの紹介 (発表者 @berobero11)」
大トリも、続いて@berobero11さんの発表でした。TrueSkillモデルというのは1対1対戦の競技などで、個人の持つ勝負ムラを反映して実力を評価するモデリング手法です。
これらの1VS1の勝負の世界ではレーティングメソッドが既に確立されていますが、強い人もはじめは低い初期レートに設定される、勝負ムラを考慮できていない、成績上位者と下位者の評価が落ち着かない等の欠点があるそうです。
そこでMicrosoftの研究者によって開発されたのがTrueSkillモデルです。例えば女子卓球を例に挙げると、福原愛選手と平野早矢香選手では、福原愛選手はムラが大きく格上にも勝つが、格下にも負けるそうです。平野早彩香選手は格上には勝つ事が少ないが格下に負けることも少ないそうです。単なる選手のパフォーマンスレーティングではなく、こうしたバラツキも含めて推定を行う事が出来るのがTrueSkillモデルです。
発表では簡単な例としてドラゴンボールにおけるキャラ同士の勝敗表からTrueSkillモデルを当てはめたユニークな例に始まり、将棋の勝敗データにこれを当てはめた実際の解析例も示して頂きました。質疑応答でも実際にモデリングするときの課題やTipsについて面白おかしく説明して頂き、とっても勉強になりました！
次回予定
懇親会の中で次回が年内にもう一度開かれることになりました。今後BUGS/Stanはマーケティングの主流になると思われそうな階層ベイズ、より複雑なモデリングを行う最高の手段だと思います。@berobero11さんのような素敵な野生のスモールデータサイエンティスト目指して、みなさんも是非参加されてみてはいかがでしょうか。
会場の様子
感想
正直、少し難しかったですがなんとか話にはついていけたこともあり、自分にはちょうど良いくらいの内容の勉強会でした。
個人的には質疑応答での議論が活発だったのも良かったかと思います。
課題としては、最近データ解析初めてみましたという方には難しそうなので、今後BUGS/Stanのユーザの裾野が広がるような内容も取り扱っていきたいと思います。
次回は初心者セッションもある予定ですし、近いうちに私の発表資料を日本語訳してチュートリアルを作成して公開したいと思います。
関係者の皆様、会場をご提供してくださった株式会社ドリコム様本当にありがとうございました。
ドラマの視聴率を予測する（２）〜stanを用いて、キャストの視聴率を算出する。
こんにちは、進撃のアナリストです。
前回の記事から早4ヶ月、今回もドラマの視聴率データを用いて、キャストの視聴率を算出したいと思います。
今回は前回とは異なりstanを用います。前回はキャスト視聴率は単なる按分平均を用いて算出していました。
今回はMCMCによるベイズ推定により、キャスト自体が持つ視聴率を直接モデルに組み込んで算出します。
参考 ATNDリンク BUGS/stan勉強会 #1
ちなみにこれが私のstanデビューです。
先に概念的な話として、今回のモデルや結果について解説してから、後半で技術的な話に触れます。
概念的な話
モデルのイメージは以下の通りです。
ドラマの初回視聴率＝B1×年度にとる効果＋B2×時間帯による効果＋B3×前回ドラマ最終話（ある場合は前シリーズ）最終回視聴率＋CP木村拓哉×出演ランク＋CP中居正広×出演ランク＋・・・（俳優890人分続く）・・・
（出演ランクは出演順を表す。主演なら5、2番目は4、3番目は3、2番目は2、1番目は1となります。）
要は単なる重回帰モデルです。これまでに使った時間帯等の効果によるモデルに加え、各キャストが出ているか出ていないかのフラグに対してキャストの効果となる係数、CP○○○○を付加したものです。
前回に記載した通り、重回帰やその他のロバストであるとされている回帰モデル手法等を試したのですが全く効果ありませんでした。
CPはキャストパワーの略です。ある特定のキャストが出演している際に、キャストパワー分の視聴率が付加されるイメージです。
今回はこのキャストパワーを俳優別に知ることを目的としています。
キャストパワーは一部のキャストのみ非常に高い値を持つ、ほとんどのキャストは0に近いと考えました（対数正規分布に従う）。
その結果、各キャストが持つキャストパワーBEST30は以下のようになりました。
結果
| 順位 | 俳優 | CP(中央値) | CP標準偏差 |
| 1 | 木村拓哉 | 2.16 | 0.23 |
| 2 | 黒木瞳 | 1.47 | 0.23 |
| 3 | 米倉涼子 | 1.44 | 0.19 |
| 4 | 仲間由紀恵 | 1.36 | 0.19 |
| 5 | 堂本剛 | 1.32 | 0.22 |
| 6 | 小雪 | 1.29 | 0.2 |
| 7 | 鈴木京香 | 1.29 | 0.21 |
| 8 | 豊川悦司 | 1.25 | 0.22 |
| 9 | 大塚寧々 | 1.25 | 0.24 |
| 10 | 竹野内豊 | 1.17 | 0.16 |
| 11 | 高橋克典 | 1.17 | 0.17 |
| 12 | 藤原紀香 | 1.16 | 0.2 |
| 13 | 中居正広 | 1.12 | 0.2 |
| 14 | 夏川結衣 | 1.11 | 0.2 |
| 15 | 坂口憲二 | 1.11 | 0.19 |
| 16 | 志田未来 | 1.1 | 0.19 |
| 17 | 香取慎吾 | 1.08 | 0.19 |
| 18 | 佐藤浩市 | 1.07 | 0.2 |
| 19 | 福山雅治 | 1.07 | 0.21 |
| 20 | ともさかりえ | 1.01 | 0.19 |
| 21 | 木村佳乃 | 1.01 | 0.2 |
| 22 | 森山未來 | 1 | 0.2 |
| 23 | 手越祐也 | 0.98 | 0.2 |
| 24 | 大沢たかお | 0.98 | 0.2 |
| 25 | 山田孝之 | 0.98 | 0.17 |
| 26 | 西島秀俊 | 0.97 | 0.19 |
| 27 | 黒木メイサ | 0.96 | 0.18 |
| 28 | 酒井若菜 | 0.93 | 0.19 |
| 29 | 天野ひろゆき | 0.9 | 0.22 |
| 30 | 田中美佐子 | 0.9 | 0.2 |
前回に比べると少しましになったという印象がありますが、ここに掲載していないものも含めて１を切るあたりから妥当性が怪しくなってきます。
とりあえずドラマに出演している回数が少ない俳優さんは排除したほうが良いかなというのが私の感想です。
しかしながら、その他の係数が想定していた数値と異なっていました。
以下はB1、B2、B3の中央値です。
> median(extract(fitmerged)$beta1)
[1] -1.425487
> median(extract(fitmerged)$beta2)
[1] -1.741203
> median(extract(fitmerged)$beta3)
[1] 2.380221
全て、正の値で１付近で収束することを期待していたのですが、時間帯と年度の効果が思いっきり負の値に振り切れています。
それを補うように前回ドラマの視聴率の係数が２を越えています。。。
時間帯や年度の効果をキャスト効果と同様に１を基準に考えたいのでこの辺りは完全に失敗です。
技術的な話
モデルを表すstanファイルは以下の通りです。
data {
int N;
int C;
vector[N] rating;
vector[N] trend;
vector[N] timezone;
vector[N] ratinglastdrama;
matrix[N, C] castmat;
}
parameters {
real beta1;
real beta2;
real beta3;
vector<lower=0>[C] mu;
}
transformed parameters {
vector[C] castbeta;
castbeta <- log(mu);
}
model {
vector[N] q;
beta1 ~ normal(1,1);
beta2 ~ normal(1,1);
beta3 ~ normal(1,1);
for(j in 1:C){
mu[j] ~ normal(-2.5,0.5);
}
for (i in 1:N) {
q[i] <- beta1 * timezone[i] + beta2 * trend[i] + beta3 * ratinglastdrama[i] + dot_product(castmat[i],castbeta);
rating[i] ~ normal(q[i], 2);
}
}
最初はBUGSを使っていたのですが、遅すぎて話にならずstanに変更しました。
実行時間の目安は以下の通りです。
| 言語 | 試行回数 | チェイン数 | 時間 |
| BUGS | 1500 | 3 | 36h |
| stan | 1500 | 4 | 30m |
圧倒的に早くなり、また収束のスピードもstanの方が早く感じました。
環境は以下の通りです。
用いた環境
課題
今回のはまだまだモデルとして未成熟で、stanを動かすテストです。次回までに以下のような課題を少しずつクリアーしたいきたいと考えています。
- 変数の取り入れ方・・・時間帯や年度の効果に関して今は代表値を用いて一つの変数にしていますが、フラグを利用した変数に変換して各それぞれの時間帯や年度のフラグに基づく効果を純粋に抽出してもよいかと考えています。色々工夫が必要なんだろうなと思います。
- 状態空間モデル化・・・俳優が旬であるかどうかをモデルに組み込みんだり、月9などに主演した俳優が高い評価を得てさらに次のドラマに抜擢されるような仕組みを導入したいと思っています。
stanを動かすまでにかなり苦労したのですが、これであれば現実的な計算時間でモデルをいじることができるので今からこれを使ってモデルの妥当性を上げていく作業を行いたいと思います。
