elsur.jpn.org >

2019年11月15日 (金)

 用事があって、日帰りで長野の軽井沢というところに行った。有名な避暑地だが、これまでご縁がなく、記憶の限りでは初めての訪問である。
 思ったより早く帰れることになったのだが、夕方の軽井沢駅は池袋駅並みにごった返しており、切符の都合で1時間半ほど時間を潰さなければならなくなった。駅の北側に少し歩いて、運よく喫茶店に空席をみつけた。
 本でも読もうと思ったんだけど、ウェイターのお兄さんがコーヒーと一緒になにかをテーブルに置き、「時間つぶしにでもどうぞ...」という。それは木製の小さなパズルであった。ボードの上にイラストが描かれた小さなコマがいくつかあり、使い込まれてかすかに茶色い光沢を帯びている。コマを上下左右にずらし、目指す模様をつくるパズルだ。子供の頃にやったことがある。
 「ははは」と愛想笑いしたけれど、そんなに暇そうに見えたかなあ、とちょっと気になった。もっともこれは私の気の回しすぎで、他のテーブルのお客さんにも渡していたようである。隣のテーブルの母娘連れも、親戚の悪口などを云いながらひとしきりこのパズルで遊んでいた。

 鞄から本を取り出したが、こんな子供だましのパズルも懐かしいものだ、などと思い、頬杖をつき、戯れに指でコマを適当に滑らせてみた。完成したら本でも読もうか、と。
 ... それから1時間半。読書どころか、ぶっつづけでこのパズルに取り組む羽目になった。
 解けない!解けないよ!
 このパズル、めちゃくちゃ難しいじゃん!!

 時間切れでお店を出るときにはイライラが極限に達しており、レジの横に並べられたこのパズルをひっつかんで買ってしまった。1300円、思わぬ出費である。店主らしき方が笑っていた。
 帰りの新幹線でもこのパズルがぐるぐると頭をめぐって離れなかったのだが(「あの4つのリンゴを隣り合わせにして...」云々)、いやまてよ、と我に返った。レジ前の陳列には「日本一難しいパズル」と謳われていた。このパズル、たぶんほんとに難しい奴なのだ。いっぽう、私はあんまし頭の良いほうじゃない。パズルなんて大の苦手だ。なにも私が頑張って解くことはないんじゃないか。
 むしろ、パズルを解くプログラムを書いたほうが早いんじゃないか?

問題
 買ってきたパズルはこちら。

DSC_0037.JPG

 コマをボードの外に出さずに移動させ、左上のこびとを右下の角に3つの柵でとじこめる、というのがゴールである。
 パズルももちろん著作物であり、勝手にブログに載せてはいけないと思ったのだが、検索してみたところ、喫茶店の方が画像を公開しておられることに気が付いた。これって、喫茶店の店主の方のオリジナル作品なのかしらん? すごいなあ。
 お店の名前は「喫茶館丹念亭」。とてもよい雰囲気のお店で、コーヒーもシナモンケーキも美味しかった。いつかまた軽井沢に行く機会があったら立ち寄りたいと思います。

 うちに帰ってビールを飲みながら再び試してみて、このパズルにはまった理由がなんとなくわかってきたような気がした。
 このパズル、なんとなくコマをすべらせていると、なんとなく盤面が変化する。思いもよらぬ形で大きなコマが動いたときなど、ちょっとした快感がある。しかしよくみると、全く展望が開けていない。同じような局面をぐるぐると回っているだけで、そこから全く脱出できないのである。まるで迷路にはまり込んでしまったような気分だ。
 よくわからんが、こういうパズルを考案するのってきっとすごい才能と努力の結晶なのであろうと思う。ほんとに尊敬します。

 プログラムで解くにあたっての方針は次の通り。

  • Rで書く。この課題にRが向いているからではなくて、単に私がRに慣れているからである。
  • このパズルに限らず、このタイプのパズルすべてに適用できるプログラムにする。Wikipediaによれば、こういうパズルのことをスライディング・ブロック・パズル、ないしスライド・パズルと呼ぶらしい。知らなかった。
  • 初期状態(つまり、最初のコマの配置)から出発し、そこからコマを動かして到達し得る状態を、力づくですべて探索する。なぜなら、どういう評価関数を書いたらいいのか、私にはよくわからんからである。
  • 得られたすべての状態からなるネットワーク・グラフを生成する。パズルの解を求めるコードは書かない。だって、グラフが手に入っていれば、ネットワーク分析のソフトを使えば最短経路を探索できるじゃん?
  • 朝までに書き終える。

 我ながら頭の悪さがにじみ出る方針である。特に最後の奴がそうだ。結局、翌朝までには書けませんでした。

 上で状態のネットワーク・グラフと書いたが、これはこういう意図である。
 このパズルについていえば、最初にできるのは右上の柵を上に動かすことだけである。2手目にできることは2つある。右端列の下から2番目のリンゴを上に動かす、ないし、右から2つめのリンゴ2個のピースを上に動かしていくことである。3手目にできることは...
 ... と考えていくと、どんどん枝分かれしていくが、単純に増えていくだけではない。元の状態に戻ったり、別の手順でたどり着いた状態に一致したりすることもありうる。

exec4.png

 上図は5手までの状態を表現するネットワーク・グラフである。紫色のノード(円)は初期状態を表している。24Uとは「行2, 列4にあるコマを上に動かす」こと、つまり右上の柵を上に動かすことを意味している。
 このように、5手までで14個の状態に到達できるわけだ。到達できる状態の数はどんどん増えていく。これをずーっと続けていけばさ!いつかはゴールに辿りつくよね!
 後述するが、こういう過度に楽観的な態度、問題の難しさを最初に客観的に推し量ることを怠る姿勢が、人生の悲劇を招くのである。

これがパズルの解答だ
 初期状態から到達し得る状態をすべて調べてみた。ただし、目標条件を満たす状態に辿り着いたら、そこからさらに探索するのはなし、ということにした。つまり、何も考えずにピースを動かしてはいるものの、ゴールに辿り着いたことに気がつかないほどアホではないプレイヤーが、いったいどのような盤面に辿り着きうるか、網羅的に調べたわけである。
 軽井沢で買ったパズルの場合、手元のパソコンで約8時間で探索が終了した。8時間。映画を4本見られる時間、本を数冊読める時間、精力絶倫な王様なら子供をたくさん作れる時間だ。

 すべての状態のネットワーク・グラフを描いてみた。それではご覧ください。軽井沢駅北口の喫茶店の店内でパズルのコマを動かしているみなさん! あなたはいま、この図の中のどこかにいますよ!

exec4_2_1.png

 ネットワーク・グラフどころか、なにか模様のようなものしかみえないが、これはカメラをすごく引いているからで、目を凝らせば見えるであろう微小な点が個々の状態を表している。状態の数は977,412個。うち1個は初期状態、14個はゴールの条件を満たす条件である。

 これではあまりに見にくいので、初期状態からゴールまでの経路(正確に言うと、14個のゴールそれぞれへの最短経路)の途中に位置する状態、計3,028個だけを抜き出して描いてみると...

exec4_2_2.png

青いノードが初期状態、赤いノードがゴール状態である。青いノードをお探しですか? 図右側のちょっと下寄りにあります。拡大してみましょう。

exec4_2_3_rev.png

 要は、この青いノードから辿って行けば、誰でもゴールにたどり着けるわけだ。
 ああよかった、これで自分でパズルを解く必要がなくなった。夜もぐっすり眠れるというものだ。

 ついでにいうと... 野暮なのでパズルの解答は書きませんけど、このパズルにおけるゴールへの最短経路は217手である。
 買ってきたパズルには親切にも模範解答を図示した紙が添付されているのだが、その手数は254手(181手と記載されているが、ここではひとつのコマの1マスの移動を1手と数える)。おかげさまで、模範解答よりもすこし手数の少ない解答を見つけたことになる。ほんの少しだけ達成感があるが、いやいやそれよか、パズルを考案する人のほうがずっとすごいですよね。

だからどうした
 だからどうしたんだこの暇人め、といわれると困ってしまうのだが、今回は自分なりにいろいろ学ぶところ大きかった。
 帰りの新幹線のなかでぼんやり思ったのは、どうせ盤面の状態は限られているだろう、ということであった。盤面の大きさは4x5マス、コマの数は10個。喫茶店内でコマをあれこれ動かしているときも、なんだか同じ状態を行きつ戻りつしているような気がしてならなかった。だったら私にも全探索できちゃうんじゃない?と思ったのである。
 振り返るとこれは無謀な賭けであった。組み合わせの数の大きさというものは人間の直感を遥かに超える。今回のパズルはたまたま初期状態から到達できる状態の全探索が可能であったし、あとで試算したところでは、そもそもこのパズルで可能な状態の数は1,281,550個しかない。しかし、問題が少し異なるだけで、状態の数は爆発的に増大し、まともな時間では到底全探索できないものとなっていただろう。
 なんについてもいえることなんだろうけど、最初に問題の大きさを推し量ることって、すごく大事ですね。どうもそういう姿勢が私には欠けているようだ。反省。

なぜ難しいのか
 それにしても不思議なのは、このパズルがなぜこんなに難しいのかという点である。実際に遊んだことがない方に、解のネタバレを避けつつ説明するのはすごく難しいんですが、とにかくほんとに難しいんですってば。イライラするんですってば。
 難しい理由のひとつは、もちろん、ゴールまでの手順の多さであろう。しかしそれだけが要因ではないという気がしてならない。なんというか、いくらコマを動かしても、同じような局面をぐるぐると回っているように感じられるのである。

 こうした難しさを、ネットワーク・グラフの形状によって定量的に示すことはできないだろうか。
 同じような局面からなかなか脱出できないということは、ネットワーク・グラフのあちこちにノードの密な固まりがあって、その固まりから外に出て行くエッジが少ない、ということだろう。とすれば、ネットワーク・グラフ全体を通じてノードが局所的に固まっている程度を測ればよいのではないか。このパズルでは、その値がとても高くなっているのではないか?

 思い悩んだ末、次のようなことを試してみた。
 いま、ネットワーク・グラフのノードをいくつかのグループにわけたとする。ネットワーク分析の慣例に従い、これをコミュニティと呼ぶことにする。
 ノードの間のエッジは、同一コミュニティ内でつながるエッジと、コミュニティをまたがるエッジとにわかれる。前者のエッジの比率が大きいとき、コミュニティはネットワークのなかの密度の高い固まりをうまく捉えていることになる。
 そこで、コミュニティ内のエッジの比率が、そのグラフのエッジを一旦すべて消し、かわりに同じ数のエッジをランダムに張ったときに得られる同じ比率と比べてどれだけ大きいかを調べる。これをモジュラリティという。うーん、下手な説明だが、要するに、モジュラリティが大きいことは、ノードがコミュニティの中で高密度に固まっていることを示しているわけである。
 そこで、ネットワーク・グラフを徐々にたくさんのコミュニティに分割していき、そのたびにモジュラリティを測ってみよう。
 どうやって分割するか。いろんな方法があるんだけど、ここでは手っ取り早く、Clausetたちの手法を使うことにする。これ、計算がすごく速いんです。

 比較のために、すごくシンプルなスライドパズルについても同じことを試してみる。盤面は3x3, コマは8つ、サイズはすべて1x1であるパズルを考える。つまり、下図を3x3に縮めたようなパズルです。

15-Puzzle.jpg
By Micha L. Rieser - Own work, Public Domain, Link

このパズルについて、初期状態から辿り着くことができるすべての状態(181,440個)のネットワーク・グラフを作った。で、こちらについてもコミュニティ数を増やしながらモジュラリティを調べた。

 結果は下図の通り。

exec6.png

 丹念亭のパズルではコミュニティ数816、3x3のパズルではコミュニティ数76のとき、モジュラリティが最大になった。
 ここで注目したいのはモジュラリティの違いである。丹念亭パズルのモジュラリティは最大で0.993と、法外に高い。つまり、ネットワーク・グラフの描画からは読み取れないが、あたかも少量の水に溶かした小麦粉のダマのように、ノードが小さくて密度の高い固まりを形成しているのではないかと思う。そのため、やみくもにコマを動かしても、その固まりの外側になかなか出られない... ということではないかしらん?

疑問点
 ... などと書いてはみたものの、正直なところ、あまり自信がない。

  • 3x3パズルは丹念亭のパズルよりも状態数が小さいし(9!=362880通り)、ネットワーク・グラフの作成にあたって目標状態を設定せずに全探索しているので、適切な比較になっているかどうかわからない。
  • Clausetたちの手法は凝集的なコミュニティ検出手法なので、所与のコミュニティ数についてモジュラリティ最大のコミュニティを検出していないかもしれない。
  • なにより、ここでモジュラリティに注目したのは良いやり方だったのだろうか。ネットワークがクラスタ化している程度を調べる際にはクラスタ係数と呼ばれる指標を使うことが多いと思うんだけど、スライド・パズルの状態のネットワーク・グラフはクリークを持たないので(3手で元に戻ることはできないから)、クラスタ係数は0となる。こういうネットワークについて、それがクラスタ化されている程度を調べるためにはどうしたらいいんだろう?

 ううむ... このへんで時間切れということにしておくけど、いろいろ勉強が足りないことがわかった。

rSlidePzlパッケージ
 軽井沢駅前の喫茶店でパズルを渡された副産物として、任意のスライド・パズルについて初期状態から到達しうる状態を全探索しネットワーク・グラフを生成するRパッケージ rSlidePzl を作りました。
 ああ、私には聞こえる。任意のスライド・パズルについて初期状態から到達しうる状態を全探索しネットワーク・グラフを生成したいと思っている、世界三千万人の善人男女の歓喜の叫びが。(嘘です聞こえません)

 パッケージはgithubからインストールできる。
 このパッケージを使えば、お手元のスライド・パズル(9x9マス以内)について、初期状態から到達し得る状態をすべて探索し(幅優先探索)、みつかった状態からなるネットワーク・グラフをつくることができる。
 せっかくなので、初期状態からゴールへの最短手順を出力する関数もつくっておいた。また、得られたネットワーク・グラフを(それが小さめであれば)簡単に描画できる関数もつくった。

 探索すべきすべての状態を探索し尽くす場合、とんでもなく時間がかかることを覚悟する必要がある。これを避けるため、探索する手の深さ(何手目まで探索するか)や、探索する状態の総数について上限を指定し、探索を途中で打ち切ることもできる。
 このパッケージを使って得られるネットワーク・グラフは igraph パッケージの igraph クラスのオブジェクトである。igraph パッケージが提供する多様な関数によって、その性質を心ゆくまで調べることができる。

 コード例を示す。

### devtoolパッケージをインストール
### install.packages("devtools")

### rSlikePzlパッケージをインストール
# devtools::install_github("shigono/rSlidePzl")

### rSlikePzlパッケージをロード
library(rSlidePzl)

library(tidyverse)
library(igraph)

### セッティング
oSetting <- makeSetting(
## 盤面のサイズ (行サイズ, 列サイズ)
boardsize = c(4,5),
## コマのサイズ (行サイズ, 列サイズ)
piecesize = list(
A = c(2, 2), # コマタイプA.こども パズルに添付されている解答ではNo.1
B = c(2, 1), # コマタイプB.りんごとサク No.5
C = c(1, 2), # コマタイプC.横長のサク No.7
D = c(1, 1), # コマタイプD.りんご, No.2,3,9,10
E = c(2, 1), # コマタイプE.縦長のサク, No.6
F = c(2, 1), # コマタイプF.縦長のりんご No.8
G = c(1, 2) # コマタイプG.横長の林檎 No.4
)
)

### こどものピースの位置(1,1), (1,2), (2,1), (2,2)について考える
anNumStates <- sapply(
list(c(1,1), c(1,2), c(2,1), c(2,2)),
function(anLoc){
### 可能な状態の数を算出する
getNumStates(
makeState(
list(
makePiece(type = "A", loc = anLoc),
makePiece(type = "B", loc = c(NA,NA)),
makePiece(type = "C", loc = c(NA,NA)),
makePiece(type = "E", loc = c(NA,NA)),
makePiece(type = "F", loc = c(NA,NA)),
makePiece(type = "G", loc = c(NA,NA))
)
),
oSetting
)
}
)
# - こどものピースが位置(1,1)(左上端)であるときの状態数は、
# 左下端(3,1), 右上端(1,4), 右下端(3,4)にあるときの状態数と同じはずなので、4倍する
# - こどものピースが位置(1,2) にあるときの状態数は、
# (1,3), (3,2), (3,3) にあるときの状態数と同じはずなので、4倍する
# - こどものピースが位置(2,1) にあるときの状態数は、
# (2,4) にあるときの状態数と同じはずなので、2倍する
# - こどものピースが位置(2,2) にあるときの状態数は、
# (2,3) にあるときの状態数と同じはずなので、2倍する
# 各状態においてりんごのピースを置く位置は6 x 5 / 2 通りあるから...
print(sum(anNumStates * c(4,4,2,2)) * 6 * 5 / 2) # 状態数は1281550

### 開始時点の盤面の作成
oStart <- makeState(
list(
makePiece(type = "A", loc = c(1,1)), # パズルに添付されている解答ではNo.1
makePiece(type = "B", loc = c(1,3)), # No.5
makePiece(type = "C", loc = c(2,4)), # No.4
makePiece(type = "D", loc = c(3,1)), # No.2
makePiece(type = "D", loc = c(3,2)), # No.3
makePiece(type = "E", loc = c(3,3)), # No.6
makePiece(type = "F", loc = c(3,4)), # No.8
makePiece(type = "D", loc = c(3,5)), # No.9
makePiece(type = "G", loc = c(4,1)), # No.4
makePiece(type = "D", loc = c(4,5)) # No.10
)
)
stopifnot(isValidState(oStart, oSetting))

### 目標条件の作成
oGoalCondition <- makeState(
list(
makePiece(type = "A", loc = c(3,4)), # No.1
makePiece(type = "B", loc = c(1,3)), # No.5
makePiece(type = "E", loc = c(3,3)), # No.6
makePiece(type = "C", loc = c(2,4)) # No.7
)
)
stopifnot(isValidState(oGoalCondition, oSetting))

### 5手だけ検索
oGraph <- makeGraph(oSetting, oStart, oGoalCondition, verbose = 1, max_depth = 5)

### 描画
set.seed(123)
plotGraph(oGraph, method = "GGally")

### すべての手順を探索。手元のパソコンで8時間くらい
oGraph <- makeGraph(oSetting, oStart, oGoalCondition, verbose = 1)

### 目標条件を満たした状態を表示
print(V(oGraph)[V(oGraph)$status == 4])

### 目標条件を満たした状態への最短経路を表示
print(getShortestPath(oGraph)$transition)

### コミュニティ検出。fast_greedy法でmodularityを最適化する
oGraph_undirected <- as.undirected(oGraph)
oCommunity <- cluster_fast_greedy(oGraph_undirected, modularity = TRUE)

### 描画。Rでは大変なのでGephi (https://gephi.org/) で描く。

### グラフに描画のための属性を追加する

### state_color属性: {1(初期状態),2(目標状態),3(目標への途中),4(ほか)}
### Gephiでこの属性をcolorに指定すると、ノードを色分けできる
anColor <- rep(4, length(V(oGraph))) # 4: 下記以外の状態
anColor[getReachability(oGraph) == 1] <- 3 # 3: 目標状態に至るパス上にある状態
anColor[V(oGraph)[V(oGraph)$status == 4]] <- 2 # 2: 目標状態
anColor[1] <- 1 # 1: 初期状態
V(oGraph)$state_color <- anColor

### state_size属性: {1(目標への経路上にない),2(ある)}
### Gephiでこの属性をsizeに指定すると、目標への経路上のノードを大きくできる
V(oGraph)$state_size <- c(2,2,2,1)[anColor]
V(oGraph)$state_comm <- membership(oCommunity)

### state_x属性, state_y属性: x座標とy座標
### GephiのGeo layoutでこれらの属性を緯度経度に指定すると、レイアウトを再現できる。
### コミュニティ内のエッジに大きな重みを与え、Fruchterman-Reingoldアルゴリズムで
### 最適化している。Gephiでレイアウトを作成してもいいけど、こっちのほうが速い
oLayout <- layout_with_fr(
oGraph,
weights = ifelse(crossing(oCommunity, oGraph), 1, 10),
niter = 1000
)
### 座標が大きな値だとメルカトル図法では歪んでしまうので、-1から+1に基準化する
oLayout <- norm_coords(oLayout)
V(oGraph)$state_x <- oLayout[,1]
V(oGraph)$state_y <- oLayout[,2]

### graphml形式で出力。このファイルをGephiで読み込む
write.graph(oGraph, file = paste0("./full_states.graphml"), format = "graphml")

### 3x3パズルのセッティング
oSetting <- makeSetting(
boardsize = c(3,3),
piecesize = list(
A = c(1, 1), B = c(1, 1), C = c(1, 1), D = c(1, 1),
E = c(1, 1), F = c(1, 1), G = c(1, 1), H = c(1, 1)
)
)
oStart <- makeState(
list(
makePiece(type = "A", loc = c(1,1)),
makePiece(type = "B", loc = c(1,2)),
makePiece(type = "C", loc = c(1,3)),
makePiece(type = "D", loc = c(2,1)),
makePiece(type = "E", loc = c(2,2)),
makePiece(type = "F", loc = c(2,3)),
makePiece(type = "G", loc = c(3,1)),
makePiece(type = "H", loc = c(3,2))
)
)
stopifnot(isValidState(oStart, oSetting))

### 3x3パズルのネットワーク・グラフ作成
oGraph_Nine <- makeGraph(oSetting, oStart)

### 上と同じやり方でコミュニティ検出
oCommunity_Nine <- cluster_fast_greedy(as.undirected(oGraph_Nine), modularity = TRUE)

### モジュラリティ比較
oGraph_undirected <- as.undirected(oGraph)
anModularity <- sapply(
1:1000,
function(nNum){
modularity(oGraph_undirected, cut_at(oCommunity, no = nNum))
}
)
oGraph_undirected <- as.undirected(oGraph_Nine)
anModularity_Nine <- sapply(
1:1000,
function(nNum){
modularity(oGraph_undirected, cut_at(oCommunity_Nine, no = nNum))
}
)
dfPlot <- data.frame(
nNumComm = 1:1000,
nMod_1 = anModularity,
nMod_2 = anModularity_Nine
) %>%
gather(sVar, gValue, c(nMod_1, nMod_2)) %>%
separate(sVar, c("sVar1", "nVar2")) %>%
group_by(nVar2) %>%
mutate(bMax = if_else(gValue == max(gValue), 1, 0)) %>%
ungroup() %>%
mutate(fGraph = factor(nVar2, levels = c("1", "2"), labels = c("丹念亭パズル", "3x3パズル")))
g <- ggplot(data = dfPlot, aes(x = nNumComm, y = gValue, color = fGraph))
g <- g + geom_line()
g <- g + geom_point(data = dfPlot %>% filter(bMax == 1), size = 2, alpha = 0.5)
g <- g + labs(x = "コミュニティ数", y = "モジュラリティ")
g <- g + scale_color_discrete(name = NULL)
g <- g + theme_bw()
print(g)

### 以上

 githubを使うのも、Rパッケージを公開するのもはじめてなので、いろいろ不備があるかもしれませんが、そのへんはひとつ、ご愛敬ということで...

雑記:データ解析 - 軽井沢駅前の喫茶店の店員さんに渡されたパズルを解く (ためのRパッケージを作った)

 ふつうのデータ分析のなかで個体のクラスタリングをやるように、ネットワーク・データの分析の際にもノードのクラスタリングをやりたくなることがある。この分野ではクラスタじゃなくてコミュニティと呼ぶことが多いと思うけど。
 個体のクラスタリングと同様、コミュニティ検出の方法もいろいろありすぎて、私のような素人は困惑してしまうのだが、中にはやたら時間とメモリを食う奴もあれば、逆に法外に早くできちゃう奴もある。時間がかかる分にはあきらめがつくが、あまりにパッと答えが出ちゃうのは、ありがたいけどなんだか気持ち悪い...という気がする。好き勝手云ってすいません。

Clauset, A., Newman, M.E.J., Moore, C. (2004) Findng community structure in very large networks. arXiv:cond-mat/0408187v2. 
 というわけで、igraph の cluster_fast_greedy() のドキュメントで引用されていた資料を読んでみた。
 Rのigraphパッケージには, コミュニティ検出の手法としてedge betweenness, fast greedy, label prop, leading eigen, louvain, optimal, spinglass, walktrapの8つが載っているが、このうちfast greedyというのは100万ノードくらいあるグラフであっても数分でコミュニティを返してくるので、なんだか気味が悪いのである。お前真面目にやっとんのかといいたくなる。

 ノードの隣接行列を$A$とする。ノード$v$と$x$がつながっているかどうかを$A_{vw}$とする。$i=v$のときに1、そうでないときに0になる関数を$\delta(i,j)$とする。グラフのエッジ数を$m=(1/2) \sum_{vw} A_{vw}$とする。
 ノード$v$が属するコミュニティを$c_v$とする。エッジのうち、両端が同じコミュニティに落ちている割合は、
 $\displaystyle \frac{\sum_{vw} A_{vw} \delta(c_v, c_w)}{\sum_{vw} A_{vw}} = \frac{1}{2m} \sum_{vw}A_{vw} \delta(c_v, v_w)$
 となる。この量はネットワークをうまく分割できているときに大きくなるけれど、これ自体はコミュニティ構造の指標にはなれない。すべてのノードが同じコミュニティに属していれば1になるから。
 そこで、この量からランダムグラフにおけるこの量の期待値を引く。ノードの次数を$k_v = \sum_w A_{vw}$として、ランダムグラフにおいて$v$と$w$の間にノードがある確率は$k_v k_w/2m$だから、
 $\displaystyle Q = \frac{1}{2m} \sum_{vw} \left( A_{vw} - \frac{k_v k_v}{2m} \right) \delta(c_v, c_w)$
 これがみなさまご存知の modularity。経験的には、Qが0.3以上のときネットワークには顕著なコミュニティ構造がある。[←へー]

 ネットワークのコミュニティへの分割がうまくいっているときmodularityが高いなら、modularityが高くなるような分け方を探せばいいことになる。しかし、すべての分割を通じて大域的に最大のmodularityを探すのは大変なので、なんらか近似的な最適化手法を考えたい。
 というわけでですね、かつて我々はNewman(2004)でこういう手法を提案いたしました。まだコミュニティに属していないひとりきりのノードたちからはじめて、Qがもっとも大きくなるような併合を探す。これを$n-1$回繰り返すとついに全員が単一のコミュニティに入ることになる。こうして、ノードを葉としたデンドログラムができる。[←要はあれですかね、F比を評価関数にして凝集型階層クラスタリングをやるようなもんですかね]

 さて、この提案の折には、隣接行列を整数の行列として持っておいて毎回行と列をマージしていけばいいや、と思っていた。しかし巨大な疎行列の場合、ほとんどの要素が0なので、この作業は無駄である。
 そこで本日はもっと効率的な方法を御提案します。[←ええええ? そういう論文なの?]

 ... というわけで、いちおう最後までめくったけど、論文の主旨は計算を速くするための工夫の話だったので、メモは省略。よくわからんが、ネットワークの隣接行列を更新するんじゃなくて、ネットワークのコミュニティの隣接行列を更新すればよくね? というような話らしい。
 要するに... やっていることは階層的クラスタリングなんだけど、ネットワーク・データの場合は元の距離行列が疎なので、計算をめっちゃ速くするテクニックがあります、ということらしい。ふうん。

論文:データ解析(2018-) - 読了:Clauset, Newman, Moore (2004) すごく大きなネットワークにおけるコミュニティ検出

2019年11月13日 (水)

 ここんところこういう覚え書きばかりで、自分でも嫌になっちゃうんだけど...

 無向グラフにおいて、ノードiとj, jとkの間にエッジがあるときに、iとkの間にもエッジがあることを推移的であるという。有向グラフでは、iからj, jからkへのエッジがあるときに、iからkへのエッジもあることを推移的であるという。
 あるグラフにおいて、推移的であるかもしれない2つのノード(上でいうiとk)が実際に推移的である割合を推移性 transitivity という。えーと、Newman, Watts & Strogatz(2002)が提案したのだそうだ。
 Rのigraphパッケージでいうと、transivity(g, type ="global") で無向グラフの推移性を求めることができる。なお、transitivity()は有向グラフも無向グラフとして扱うようだ。エッジの向きを考慮したいんなら自分で求めろ、ってことでしょうか。

 あるノードからみて、それと隣接するノードからなるサブグラフを考え、その密度(つまり、ありうるエッジの数に占める実際のエッジの数の割合)を求める。これをそのノードのクラスタ係数 clustering coefficient と呼ぶ。igraphパッケージではtransivity(g, type = "local")で求めることができる。
 あるグラフについてノードのクラスタ係数を平均した値を、グラフのクラスタ係数という。igraphパッケージではtransivity(g, type = "average")で求めることができる。
 このクラスタ係数というのはWatts & Strogatz (1998)が提案したのだそうだ。どっちもワッツさんたちである。なんだかなあ。

 推移性とグラフのクラスタ係数はよく似ているんだけど、ちょっと異なる。
 いったいどう異なるのか。なにがなんだかわかんなくなってイライラしてきたので、仕事を中断してメモを取った。

 準備。
 無向単純グラフ$G = (V,E)$について考える。つまり、エッジに向きはないし、自己ループもないし、2ノード間には辺があるかないかのどっちかで、エッジに重みなどというややこしいものはない。
 ノード$v$に隣接するノード数を次数$d(v)$とする。

 $G$の3つのノードからなる完全サブグラフを三角形と呼ぶ。ノード$v$を含む三角形の数を$\delta(v)$とする。グラフ全体について考える場合、ひとつ三角形が増えるだけで$\delta(v)$の合計は3増えるから、グラフ全体の指標としては
 $\delta(G) = \frac{1}{3} \sum_{v \in V} \delta(v)$
がふさわしい。
 ノード$v$を中間に持つ長さ2のパスを$v$のトリプルと呼ぶ。その数を$\tau(v)$とすると、$n$個から$k$個を取り出す組み合わせを$combn(n, k)$と書くならば$\tau(v) = combn(d(v), 2)$だ。合計すると
 $\tau(G) = \sum_{v \in V} \tau(v)$

 本題に戻ろう。
 推移性とはなにか。それはトリプルの総数に占める三角形の総数だから、
 $\displaystyle T(G) = \frac{3\delta(G)}{\tau(G)}$
である。

 グラフのクラスタ係数とはなにか。
 次数$d(v) \geq 2$のノードについて、ノードのクラスタ係数$c(v)$とは、トリプルの数に占める三角形の数、つまり$c(v) = \delta(v)/\tau(v)$である。グラフのクラスタ係数というのはこれを平均した値だから、次数2以上のノードの集合を$V'$として
 $\displaystyle C(G) = \frac{1}{|V'|} \sum_{v \in V'} c(v)$
である。
 なお、次数が1までのノードについてもなんらかクラスタ係数を決め(0だか1だか)、単純平均してしまうという方法もある。

 ここで加重クラスタ係数というのを考えてみる。各ノードにウェイトとして正の実数$w(v)$が振られているとしよう。グラフの加重クラスタ係数を
 $\displaystyle C_w(G) = \frac{1}{\sum_{v \in V'} w(v)} \sum_{v \in V'} w(v) c(v)$
とする。
 仮に、トリプルの数をウェイトにしたら、つまり$w(v) = \tau(v)$としたらどうなるか。$c(v) = \delta(v)/\tau(v)$より、総和記号の右の$\tau(v)$は消えて、
 $\displaystyle C_\tau(G) = \frac{\sum_{v \in V'} \delta(v)}{\sum_{v \in V'} \tau(v)} = T(G)$
となる。
 つまり推移性とは、各ノードが持つトリプルの数を重みにした加重クラスタ係数のことだ。全ノードの次数が同じとき、ないし全ノードのクラスタ係数が一致しているとき、グラフのクラスタ係数と推移性は等しくなる。

 以上、次の論文のイントロ部分よりメモ。
 Shank, T., Wagner, D. (2005) Approximating clustering coefficent and transitivity. Journal of Graph Algorithms and Applications. 9(2), 265-275.

 なるほどね...
 グラフの特性を記述する際、推移性を使うべきなのだろうか、グラフのクラスタ係数を使うべきなのだろうか。たとえばある社会が「友達の友達もまた友達だ」的性質を持っているかどうかを調べるならば、ノードをベースに考えて後者を使い、たくさんの離散的状態を遷移していくなにかについて、そこでの遷移が「3期後には元に戻っちゃう」的性質を持っているかどうかを調べるならば、エッジをベースに考えて前者を使う、ということかなあ...

雑記:データ解析 - ちょっとした覚え書き:グラフの推移性とクラスタ係数はどうちがうか

2019年11月 8日 (金)

 最近なんだかこんなことばっかりなんだけど...
 Rで環境 (environment) と呼ばれるオブジェクトを操作したいとき、関数名がよく分からなくて困ることが多い。rlangパッケージというのを使うと、もっと使いやすい関数がわかりやすい名前で提供されていてありがたいんだけど、従来の関数との対応がわからなくて、それはそれで混乱する。
 このたびついに業を煮やして、神Hadleyが与えたもうたAdvanced R 2nd edition, Chapter 7. Environments を参照しつつメモを取った。
 すいません、私の私による私のためのメモです。

環境そのものの操作

  • 環境をつくる: rlang::env(); new.env()
  • 環境 e を表示する: rlang::env_print(e), print(e). ただしprint()では中身はわからない
  • 環境 e の親環境: rlang::env_parent(e), parent.env(e)
  • 環境 e のすべての先祖の環境: rlang::env_parents(e)
  • 空の環境(すべての環境の始祖): rlang::enpty_env()
  • グローバル環境(いわゆるワークスペース): rlang::globalenv(), .GlobalEnv
  • サーチパス(グローバル環境の先祖の環境たち)の名前を返す: search()
  • サーチパスの環境たちを返す: rlang::search_envs()
  • baseパッケージのパッケージ環境 (そのパッケージが外部に提供する環境): rlang::base_env()
  • 現在の実行環境(関数の実行中はその関数の実行環境、そうでない場合はグローバル環境): rlang::current_env(), environment()
  • 現在の関数を呼び出している関数の実行環境: rlang::caller_env(), parent.frame()
  • ある関数が実行されるときに生まれる実行環境の親は、その関数の「関数環境」である(呼び出し元の関数の実行環境ではない。ここ、ときどき勘違いしちゃいますね)。関数環境とは、その関数がつくられたときにそいつがbindしたというかキャプチャした環境のこと。対話的につくった関数の関数環境はグローバル環境、関数1のなかでつくった関数2の関数環境は関数1の実行環境となる。
  • 関数 fの関数環境: rlang::fn_env(f), environment(f)
  • パッケージのなかの関数の関数環境はそのパッケージの名前空間環境。その親はそのパッケージのimports環境。その親はbaseパッケージの名前空間環境。その親はグローバル環境。

環境の中身の操作

  • 環境 e のなかの名前をみる: rlang::env_names(e), names(e), ls(e, all.names = T)
  • 環境 e が要素 x を持っているかを調べる: rlang::env_has(e, "x"), exists("x", e)
  • 環境 e の要素 x にアクセスする: e$x, e[["x"]]
  • 環境 e の要素 x の値を得る(存在しないときにはエラーを発生させる): rlang::env_get(e, "x"), get("x", e, inherits = FALSE)
  • 環境 e の要素 x に値 1 をbindする: rlang::env_poke(e, "x", 1), rlang::env_bind(e, x = 1); assign("x", 1, e).
  • 環境 e の要素 x に表現myfun()をbindするが、そのときは評価せず、最初にアクセスしたときだけ評価する: rlang::env_bind_lazy(e, x = myfun()), delayedAssign("x", myfun(), assign.env = e)
  • 環境 e の要素 x に関数myfunをbindする。アクセスするたびに評価する: rlang::env_bind_active(e, x = myfun), makeActiveBinding("x", myfun, e)
  • 環境 e の要素 x を消す: rlang::env_unbind(e, "x"), rm("x", envir = e)。なお、e$x <- NULL では消せない。
  • 永続付値 <<- は、先祖の環境たちのどこかにある変数を変更する。どこにもみつからなかったらグローバル環境に変数をつくる。現在の環境にある変数を変更したり、変数をつくったりすることはない

雑記:データ解析 - ちょっとした覚書:Rの環境を操作する関数

2019年10月29日 (火)

 話はちがうが(なにからだ)、Stanを使っていると、ときどきデータ型について混乱してしまうことがある。数日前も同僚が書いたコードをみていて、しばしわけがわからなくなってしまった。私が悪いんです、すいません。
 松浦「StanとRでベイズ統計モデリング」の9章を熟読すればよい問題なのだが(自宅にも職場にも置いてあるんだし)、毎度毎度めくっているのもどうかという気がしてきたので、Users Guide (2.19)の15-16章を眺めてメモを取った。

15. Matrices, Vectors, and Arrays

  • Stanのデータ型は、int, real, vector, row_vector, matrixの5つ。
  • どの型も配列にできる。たとえばreal b[M,N]は2次元の配列。
  • サイズを動的に変えることはできない。
  • 2次元配列と行列を比べると、行列のほうがメモリがちょっと小さい。行列はメモリ上で固まってストアされているが配列はそうでない。
  • 行列は列優先順でストアされているので、matrix[M,N] aならさきに1:Nでループしたほうが速い。配列は行優先順ではいっているので、real b[M,N] ならさきに1:Mでループしたほうが速い。そして行列のほうが速い。
  • matrix[M,N] aに対して a[m] を呼ぶと行ベクトルが返るけど、これは遅い。row_vector[N] b[M]に対して b[m]を呼んだほうが速い。
  • 行列の行列演算は速い。
  • コンテナとしてみたとき、vector, row_vector, スカラーの一次元配列は同じ(ただし整数が入れられるのは配列だけだけど)。

16. Multiple Indexing and Range Indexing

  • たとえば int c[3] に(5,9,7)をいれといて、ind indx[4]に(3,3,1,2)をいれとくと、c[idxs]は(7,7,5,9)を返す。こういうのがmultiple indexing。
  • その特殊ケースとしてrange indexing (slicing)がある。c[1:2]とかね。c[2:size(c)]ってのもあり。c[:2]はc[1:2]の略記になる。c[]とc[:]は同じ。

16章にはほかにもいろいろ書いてあるんだけど、なんだかめんどくさくなってきたので、まあいいや、別の機会にしよう。

雑記:データ解析 - ちょっとした覚書:Stanのデータ型について

2019年10月23日 (水)

 たとえば、いまサイズ4の標本があり、ある変数の値が小さい順に100, 200, 300, 300であったとする。標本中央値は250ですわね。では40パーセンタイルは?
 これは意外にややこしい話で、以前気が付いたんだけど、ソフトによって出力が異なるのである。
 こんな細かいことに気を揉んでも一銭の得にもならん、くわばらくわばら...と思っていたが、このたび勤務先のお仕事の関係で関連するお問い合わせを受け、まあこれも機会であろうと思い、ちょっと調べてみた。

Hyndman, R.J. & Fan, Y. (1996) Sample Quantiles in Statistical Packages. American Statistician, 50, 361-365.
 この話になるとよく引用される論文はこれだという印象だけがあったのだが、実際に手に入れてみたら、4pの短いコラムのような記事であった。

 いわく。
 分布関数を$F(x)$として、分位数とは
 $Q(p) = F^{-1}(p) = inf \{x: F(x) \geq p \}, \ \ 0 < p < 1$
だが、これを標本においてどう定義するか。

 以下、$i$番目の定義を$\hat{Q}_i(p)$とする。標本を$\{X_1, \ldots, X_n\}$とし、その順序統計量を$\{X_{(1)}, \ldots, X_{(n)}\}$とする[小さい順に並べて$3$番目の値を$X_{(3)}$と書きますってことでしょうね]。
 $u$を超えない最大の整数を$L(u)$, $u$を下回らない最小の整数を$U(u)$と書く[←原文ではあまり見たことない記号を使っているので、このメモではこのように略記する。以下、$u$に負値が来ることはなさそうだから、要は$u$の端数を切り捨てたら$L(u)$, 切り上げたら$U(u)$ね]

 標本分位数の定義はソフトによってさまざまである。一般化して書くと、

 $\hat{Q}_i(p) = (1-\gamma) X_{(j)} + \gamma X_{(j+1)}$
ただし $\frac{j-m}{n} \leq p \leq \frac{j-m+1}{n}$、$m$は実数、$0 \leq \gamma \leq 1$。
 $\gamma$は $j = L(pn + m)$と$g = pn + m - j$の関数である。

 [著者らは頭が良すぎるせいか抽象化しすぎているので、具体例を使って、私にもわかるように書き換えよう。
 標本$\{100, 200, 300, 300\}$(サイズ$n=4$と書く)の40パーセンタイル($p=0.4$と書く)とは、$j$番目の値と$j+1$番目の値の間の区間のどこかにある点である。その点からその区間の左端までの長さが区間の幅の$(1-\gamma)$倍であり、その点からその区間の右端までの長さが区間の幅の$\gamma$倍であるとしよう。ソフトによるちがいは、結局は$j$と$\gamma$の決め方の問題にすぎない。
 $j$の決め方のひとつは、$np = 4 \times 0.4 = 1.6$ の端数を切り捨てた値とする、というものである。このデータ例だと$j=1$、つまり100と200の間ということになる。ソフトによっては端数を切り捨てる前に$m$を足すものもある。たとえば$m=0.5$なら$j=2$、つまり200と300の間ということになる。
 $\gamma$の決め方にはいろいろあるんだけど、どの決め方にせよ、それは$j$と、切り捨てた端数$g$のふたつによって決まる。]

 さて、標本分位数$\hat{Q}_i(p)$に期待される特性として次の6つを挙げることができる。

  • P1: 連続量であること。
  • P2: $Freq(X_k \leq \hat{Q}_i(p)) \geq L(pn)$であること。[Hyndmanさんのブログ記事の誤植リストに従って修正した。上の具体例でいうと、40パーセンタイルの左側にある値の個数が、$0.4 \times 4 = 1.6$の端数を切り捨てた値、すなわち1以上であってほしい。そりゃそうだ]
  • P3: $Freq(X_k \leq \hat{Q}_i(p)) = Freq(X_k \geq \hat{Q}_i(1-p))$ [上の具体例でいうと、40パーセンタイルより左の値の個数は, 60パーセンタイルより右の値の個数と同じであってほしい。うん、そりゃそうあってほしいわね]
  • P4: $\hat{Q}^{-1}_i (X)$が一意に定義されるとき、$k=1,\ldots,n$について $\hat{Q}^{-1}_i (X_{(k)}) + \hat{Q}^{-1}_i (X_{(n-k+1)}) = 1$ [えーと、たとえば左からk番目の値が20パーセンタイルだったら、右からk番目の値は80パーセンタイルであってほしいということだ。P3と同じく、データの扱いが左右対称であってほしいということであろう]
  • P5: $\hat{Q}^{-1}_i (X)$が一意に定義されるとき、$\hat{Q}^{-1}_i (X_{(1)}) \geq 0$かつ$\hat{Q}^{-1}_i (X_{(n)}) \leq 0$ [最小値が0パーセンタイルというのは変だし最大値が100パーセンタイルというのも変だということであろう]
  • P6: $\hat{Q}_i(0.5)$が通常定義される標本中央値(奇数個だったら真ん中の値、偶数個だったら真ん中にある二つの値の中点)であること。

 大変ながらくお待たせしました。それでは!選手の入場です!

 定義1.

$m=0$とする。$g > 0$のとき$\gamma = 1$, $g = 0$ のとき$\gamma = 0$とする。

 この定義は以下の性質を持つ。
 $Freq(X_k \leq \hat{Q}_1(p)) = U(pn)$
 $Freq(X_k \geq \hat{Q}_1(1-p)) = L(pn+1)$
 P2を満たす。また$n$が奇数の時P6を満たす。他の特性は満たさない。
 [上の具体例でいうと、$j$は$0.4 \times 4 = 1.6$の端数を切り捨てた値, すなわち$1$。切り捨てが起きているので$\gamma = 1$、よって答えは$j+1$番目の値200となる。それ以下の値の個数は2、これは$1.6$の端数を切り上げた値に等しい。なお、60パーセンタイルは300、それ以上の値の個数は2となり、これは$1.6+1$の端数を切り捨てた値に等しい]

 定義2.

$m=0$とする。$g > 0$のとき$\gamma = 1$, $g = 0$ のとき$\gamma = 0.5$とする。

 この定義は以下の性質を持つ。
 $Freq(X_k \leq \hat{Q}_2(p)) = Freq(X_k \geq \hat{Q}_2(1-p)) = U(pn)$
 P2, P3のみを満たす。
 [定義1とほぼ同じなのだが、$pn$が端数を持たなかったときの扱いが異なる。たとえば上の具体例での25パーセンタイルは、定義1では100となるが、定義2では100と200の中点すなわち150となる。なんだか奇妙な定義だが、たとえば20パーセンタイル以下の値の個数と80パーセンタイル以上の値の個数が常に一致し$U(0.2n)$となる、つまりP3の特性を持つというメリットがあるわけだ]

 定義3.

$m=-0.5$とする。$g > 0$のとき$\gamma=1$とする。$g=0$の場合、jが偶数ならやはり$\gamma=0$, 奇数の場合は$\gamma=1$とする。

 この定義はP2のみを満たす。
 [上の具体例でいうと、40パーセンタイルを求めるとして、$j$は$np - m = 1.6-0.5= 1.1$の端数を切り捨てた値1となる。切り捨てが起きたので2番目の値200を採用する。切り捨てが起きなかったときがなんだかややこしいが、まあその値だかその次の値だかを使うわけだ]

 以上の定義では、分位数は$j$番目の値か$j+1$番目の値かその中点であり、つまり非連続的である(P1を満たさない)。
 ここからは区間ごとに線形関数を考えて内挿するタイプの定義で、すべてP1を満たす。
 以下、$\hat{Q}_i (p_k) = X_{(k)}$と書く。
 [つまりですね、横軸に$X_{(k)}$、縦軸に$\hat{Q}_i (p_k) = X_{(k)}$となる$p_k$をとり、点$k=1, \ldots, n$が布置する散布図を書いて、点と点を直線でつなぐ。で、たとえば水平線$Y=0.4$との交点を求めてこれを40パーセンタイルとしようという算段である。問題は$p_k$をどう定義するかだ]

 定義4.

$p_k = k/n$として線形に内挿する。

 この定義はP1, P2を満たす。
 [上の具体例でいうと、40パーセンタイルを求めるなら、点$(100,1/4)$と点$(200,1/4)$を結ぶ線分と、水平戦$Y=0.4$との交点をもとめるわけね。40パーセンタイルは160となる。素直な発想だが、このやり方だと50パーセンタイルは200になり、通常の標本中央値と合致しなくなる、つまりP6は満たされない]

 定義5.

$p_k = (k-0.5)/n$として線形に内挿する。

 P1からP6のすべてを満たす。
 [上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1.5/4)$と点$(300,2.5/4)$を結ぶ線分と、水平線$Y=0.4$との交点をもとめるわけだ。40パーセンタイルは 210。50パーセンタイルは250となり、通常の標本中央値と合致する]

 定義4-5には分布の分位数を推定するという視点が欠けている。いっぽう定義6-8は、$F(X_{(k)})$のなんらかの特性値を$p_k$にするというもの。いま$F(X_k)$が一様分布とすれば、$F(X_{(k)})$はベータ分布$Beta(k, n-k+1)$になるので、そのなんらかの特性値を使う。[←この理屈がいまいちぴんとこないんだけど... $F(X_k)$は実際には一様分布じゃないっすよね...]

 定義6.

$p_k = k / (n+1)$として線形に内挿する。

 これは$F(X_{(k)})$の期待値である。P3以外の5つの特性を満たす。
 [上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,2/5)$と点$(300,3/5)$を結ぶ線分と$Y=0.4$との交点を求めるわけだ。40パーセンタイルは200, 45パーセンタイルだったら225]

 定義7.

$p_k = (k-1) / (n-1)$として線形に内挿する。

 これは$F(X_{(k)})$の最頻値である。P1,P2,P4,P6を満たす。
 [上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1/3)$と点$(300,2/3)$を結ぶ線分と$Y=0.4$との交点を求めるわけだ。答えは220]

 定義8.

$p_k = (k-1/3)/(n+1/3)$として線形に内挿する。

 これは$F(X_{(k)})$の中央値の近似である。P3以外の5つの特性を満たす。
 [上の具体例でいうと、40パーセンタイルを求めるなら、点$(200,1.666/4.333)$と点$( 300,2.666/4.333)$を結ぶ線分と、$Y=0.4$との交点を求めるわけってことか...答えは206.66]

 定義9は視点が違って、$F$が正規分布だと仮定し、$F(E(X_{(k)}))$を求めるというもの。定義6では$F(X_{(k)})$の期待値を考えたが、ここでは$X_{(k)}$の期待値を考える。次の定義で近似できる由。

$p_k = (k-3/8)/(n+1/4)$として線形に内挿する。

 P3以外の5つの特性を満たす。
 [上の具体例でいうと...めんどくさいので省略するが、40パーセンタイルは207.5となる模様]

 特性だけでいえば優れているのは定義5である。いっぽう推定という観点からは、定義8は$F$の形状を問わず$Q(p)$のmedian-unbiasedな推定値であり、これがもっとも優れている。
 思えば、標本分散の定義もいろいろありうるが(直観的にはnで割るのが分かりやすいし、正規分布ならn+1で割ればMSE最小になる)、不偏性を重んじてn-1で割るというのが標準になっている。いっぽう標本分位数の定義は標準化されていない。由々しいことだ。ソフトメーカーのみなさん、定義8に統一しませんか。

 ... 細かいところはよくわかんないんですけど、勉強になりましたです。
 ところで、第1著者Hyndman先生は2016年のブログ記事でこの論文について振り返っている。いわく、この論文はある時期からすごく引用されるようになった(そのきっかけのひとつはWikipediaに論文へのリンクが載ったからだが、実は先生がご自分で載せたのだそうだ。ここ笑っちゃった)。しかし思いの丈は果たされなかった。各ソフトウェアは標本分位数のデフォルトを定義8に統一するのではなく、むしろこれまで実装していなかった定義も実装するようになってしまった... とぼやいておられる。面白いなあ。

 現時点での各種ソフトはどうなっているだろうか。
 Rの場合。そもそもRのstats::quantile()を書いた人のひとりがHyndman先生その人であるからして、この論文の定義の番号に合わせ、たとえば quantile(x, type = 8) という風に書けばよい。ただしデフォルトは先生の意に反し、S言語での定義に合わせた定義7となっている。

 私は使っていないけど、時々お問い合わせを受けるのがSPSSの挙動である。どうなっているのか。
 手元にあるIBM SPSS 24のStatistic AlgorithmsというPDFによると(ここからダウンロードできる)、単変量の記述統計量を求めるEXAMINEコマンドは、分位数として5つの定義を実装している。FAQページにいわせるとデフォルトはHAVERAGEという定義だ。
 PDFによればHAVERAGEの定義は次の通り。

標本サイズを$W$, $p$パーセンタイルを$x$, 重複を取り除いた$i$番目の値を$y_i$, その度数を$c_i$, 累積度数を$cc_i$と表記する。$tc_2 = (W+1)p$とし、$cc_i$で区切られる区間のうち$tc_2$が落ちる区間の下限を$cc_{k_2}$とする。$g^{*}_2 = tc_2 - cc_{k_2}$, $g_2 = g^{*}_2 / c_{k_2+1}$として、

  • $g^{*}_2 \geq 1$のとき $x = y_{k_2+1}$
  • $g^{*}_2 < 1$ かつ $c_{k_2+1} \geq 1$のとき$x = (1-g^{*}_2)y_{k_2} + g^{*}_2 y_{k_2+1}$
  • $g^{*}_2 < 1$ かつ $c_{k_2+1} < 1$のとき$x = (1-g_2)y_{k_2} + g_2 y_{k_2+1}$

うむむ、ややこしい。
 話がややこしくなっているのは、SPSSではローデータの個々の行に対して整数でない頻度ウェイトを与えることができるからである。
 話を簡単にするために、頻度ウェイトは1以上の整数としよう(頻度ウェイトとは標本におけるその行の頻度であるから、本来は1以上の整数になるはずである)。このとき$y_i$の度数$c_i$は1以上の整数になるので、3本目の式は無視できる。
 1本目、$tc_2$が区間下限$cc_{k_2}$より1以上離れていたら区間の上限を採用するというのは、たとえば$\{100, 200, 300, 300\}$の65パーセンタイルを求めるとき、$tc_2 = 0.65 * (4+1) = 3.25$が落ちる区間は累積度数2と累積度数4の間だけど、下から3人目以降はあきらかに300なんだから、$tc_2$が3以上だったら四の五の言わずに300だ、ということだと思う。SPSSでいう$y_i$はHyndmanさんのいう$X_{(i)}$とちがって、データの重複を取り除いたうえで小さい順に並べたときの$i$番目の値だから、こういう但し書きが必要になるわけだ。
 話のポイントは2本目の式である。最初に挙げた具体例$\{100, 200, 300, 300\}$で考えると、その40パーセンタイルは、$tc_2 = 0.4 * (4+1) = 2$, これが落ちる区間の下限は2だから$g^{*}_2 = 0$で、2番目の式より200となる。45パーセンタイルなら$tc_2 = 0.45 * (4+1) = 2.25$, $g^{*}_2 = 0.25$となり、0.75*200 + 0.25 * 300 = 225。どちらの例でもHyndmanさんのいう定義6に一致する。
 それもそのはずで、Hyndmanさんの説明は「点$(200,2/5)$と点$(300,3/5)$を結ぶ線分と水平線$Y=0.45$との交点を求める」であり、SPSSの説明は「点$(200,2)$と点$(300,3)$を結ぶ線分と水平線$Y=0.45 \times 5=2.25$との交点を求める」である。データ点の縦座標を5で割るか、欲しい$p$に5を掛けるかの違いで、結局同じ事をいっている。
 このように、頻度ウェイトが正の整数であれば、SPSSの出力は定義6だと思う。

 Excelはどうなっているんだろう... 上の具体例で求めてみたところ、セル関数PERCENTILE.INCは40パーセンタイルとして220を返し、PERCENTILE.EXCは200を返した。調べる気力が尽きたが、これ、なにやってんですかねえ。

論文:データ解析(2018-) - 読了:Hyndman & Fan (1996) 標本分位数の出力はソフトによってちがう

2019年10月21日 (月)

 ちょっと事情があって調べ物をしていたら、American Business History CenterというWebサイトに、アメリカの消費者パネルデータの歴史についてのすごく面白い記事が投稿されていた。仕事そっちのけでメモをとってしまったのでブログに載せておく。

 著者のFred Phillipsさんはどうやら現在ニューメキシコ大の教授らしい。著書に"Market-Oriented Technology Management"があり、このブログ記事はその一部の改訂版である由。

 いわく。
 メーカーは自社の出荷についてはわかっているが、(他社製品を含めた)消費者の購入実態についてはわからない。そのため、FMCG(日用消費財)メーカーはパネル調査を頼りにする。消費者をパネルとして抱え込んでおいて複数時点で調査する方法だ。
 80年代以降、アメリカの消費者パネル調査はA.C. Nielsen Co.とMarket Research Corporation of America (MRCA)に支配されていた。Nielsenは消費者のメディア習慣のデータを、MRCAは消費者の購買データを売っていた。アメリカの主なFMCGメーカーは、たいてい両方のデータを買っていた。
 これはNielsenとMRCAの栄枯盛衰の物語である。

 1936年、Arthur C. Nielsen, Sr.はMITからラジオ聴取を記録する装置のライセンスを受けた。彼は1942年、Nielsen Radio Indexを開始した。
 戦後のTV普及に伴い、NielsenはTV視聴率を郵送の日記調査で調べるようになった(この日記調査は80年代まで続いた)。さらにNielsenは据え置きデバイスを併用するようになった。
 といっても、据え置きデバイスでは世帯内の誰がTVを観ているのかわらない。視聴者にいちいちボタンを押してもらい、誰が観ているのかを調べようとしたが、これは協力率が低かった。ひとりひとりに発信器つきのメダルをぶら下げてもらうというのも試したがこれも無理だった。据え置きデバイスに熱センサーを付けてみたが、犬や赤ちゃんやストーブやトースターに反応してしまった。

 1941年、ゲーム理論研究者のOskar MorgensternはIndustrial Survey Inc.を創業する。この会社はまもなく社名をMRCAに変えた。
 [なによりこの記載に度肝を抜かれた。フォン・ノイマン & モルゲンシュテルン「ゲームの理論と経済行動」のモルゲンシュテルンだよね? まじか。この人、ナチス占領でアメリカに亡命したウィーンの知識人だから、きっと生活費を稼ぐためにいろんなことをしたんだろうなあ...]
 MRCAは訪問調査を行っていたが、50年代の都市化とともに訪問調査は難しくなり、日記調査に軸足を移した。
 70年代、女性が職場に進出し、離婚が増え、世帯サイズが小さくなり、日記調査への協力率は下がった。しかしMRCAのデータは依然として出荷データのトレンドに対応していた。この頃、David B. Learner(心理学出身、広告代理店BBDOの元幹部)がMRCAを買い、70年代末には単独の所有者となった。
 [おいおいちょっと待って... Clancy, Krieg, & Wolf (2006)の歴史の章には、David Learnerは「60年代にBBDOでDEMONというSTMモデルを開発するも、うまくいかずBBDOを去ってピッツバーグのハイテク企業の社長になる人」として登場するんですけど?]

 さて70年代末。いよいよ小売店に清算用スキャナが普及し始める。
 1979年、野心的スタートアップ Information Resources, Inc. が登場する。この会社はアメリカ全国から人口構成が全国のそれに近い複数の街を選んで"pod markets"とし、そこのすべてのスーパーにスキャナを配り、世帯を抽出して調査を掛けたうえ、IDカードを使って精算時にスキャンするように求めた。さらにケーブル局と提携し、街を半々にわけて別のCMを流せるようにした。
 IRIにもいろいろ問題はあった。そもそもすべての製品にUPCコード[日本で言うJANコードね]が振ってあるわけではない。味覚やブランドへの選好には地域差があるので、人口構成が似ているからと言って購買データを全国にあてはめることはできない。パネルのメンバーはIDカードを家に置き忘れるかもしれない。pod marketsがどこにあるかは公開されていたから、競合はプロモーションを打って広告テストを無効化するかもしれない。スキャナ・パネル・データからは原則として競合の価格はわからない。スキャナがあるのはスーパーだけで、他の小売チャネルはわからない。当時のスキャナ・データは結構エラーが多かった。
 というわけで、メーカーは当初はIRIに飛びついたものの、いったんはMRCAに戻ってきた(典型的なハイプ・カーブである)。結局はIRIが勝つことになるわけだけど。

 80年代は新技術の時代であった。RDD, マイコン, コンピューター産業調査のためのSurvey-on-a-disk [おそらくフロッピーを配ってコンジョイント分析をやるというような話であろう]、自動コールセンター。
 NielsenはIRIに対抗し、世帯パネルに手持ちスキャナを送るという作戦にでる。
 実はMRCAも同じことを考えていた。80年代初頭、私[著者のFredさん]はMRCAの副社長で、家庭でのスキャンを実験したんだけど、なかなかうまくいかなかった。子どもがはしゃいでしまって、シリアルの箱をなんどもスキャンするのである。[←ははは]
 社長のDaveはいつも云っていたのだが、彼はMRCAの経営に満足しており、上場する気はさらさらなかった(上場などしようものならウォール街のアナリストと延々電話する羽目になる)。スキャナを配るには資金が必要だ。Daveはスキャナ・パネルと競合するのをやめて、データ収集ではなく企業へのデリバリーに特化しようとした。社名をMRCA Information Servicesに変え、調査対象カテゴリを拡大した(金融サービスとか)。旗艦商品は即時的かつインタラクティブな市場調査ツール。おかげさまで顧客も戻ってきた。
 私は1988年にMRCAを去った。円満退社であった。

 90年代はさらなる新技術の時代、そして技術市場の国際化の時代であった。Nielsenは欧州と日本に進出する。日本ではそれまで、朝日新聞のような日刊紙の消費者パネルが使われていた。
 fax、そして電子メールによる調査が始まった。TVとWWWの融合が始まった。画像認識によるTV視聴者判定、音声認識による調査が登場した。
 NielsenとIRIは技術を向上させ、いよいよMRCAからシェアを奪い始める。もっとも同時期、IRIが600万ドルかけたシンジケート広告調査は大失敗したし、1996年にはTV局4社が共同でNielsenデータの不正確さを批判する広告を出したのだけれど。

 IRIとNielsenは、店舗スキャナ・データとスキャナ・パネル・データを併売するようになる。後者はなんと実質無料。独禁法に反するのではないかという疑いもあったが、MRCAに法廷闘争の体力はなかった。クライアントもスキャナ・パネルの怪しさを知っていたが、低価格には抗えなかった。こうして日記調査は急速にシェアを失う。
 近視眼、拒否、警告、安堵、自己満足、パニック、そして諦め。MRCAはまさにこの路を辿った。MRCAは資金不足で税の不払いに陥った。幹部社員は有罪を認め、短期収監と罰金刑に服した。[←ええええええ]
 かくしてMRCAは滅びた。

 といっても、日記調査そのものは滅びていない。日記調査はいまでも、コードが付いていないカテゴリ(服飾・靴・宝石、家具など)の購入や、食品の消費を追跡するための最良の手法である。
 Nielsenもまた買収・分割された。1984年、Dun & Bradstreet CompanyがNielsenを買収。1996年、TV視聴率測定のNielsen Media Researchと消費者購買・映画興収のAC Nielsenに分割。1999年、オランダの出版社VNUが前者を買収。2001年、VNUは後者も買収して併合。
 皮肉なことに、IRIもかつてのMRCAと同様、データ収集ではなく分析を志向する企業となった。

 90年代初期までには、食品の約60%が小売店での清算時にスキャナを通ることになった。スキャナ・データによって、メーカーから小売店へと権力がシフトした。全国テレビ局から地域の広告主、そして消費者へと権力がシフトした...[中略]
 調査も変わった。厳密な無作為標本ではなく任意型標本が支配的となった...[中略]

 20世紀にMRCAとNielsenが築いたデータベースは、いまからみれば笑ってしまうくらい小さなものだった。しかし、それらから現代に通じる教訓を得ることができる。

  • データの買い手は、「高品質だが高い」データよりも「疑わしいが安い」データよりも好む。
  • 悲惨な戦略的決定を行うことは容易である。
  • 企業はこれからも、新技術によって馬鹿げた失敗を繰り返すだろう。かつては熱センサーとメダル。いまならドローンやソーシャルメディア。
  • 人々はこれからも、たいした見返りもないのに個人データを企業に手渡すだろう。

 [...以下略]

 。。。ところで、モルゲンシュテルン創業のIndustrial Survey Inc.がMRCAに社名変更するというくだり、Jones & Tadajewski(2011)によれば次のような事情らしい。
 そもそもMRCAという会社を創業したのは最初期のマーケティング研究者Percival Whiteである。彼はPauline Arnoldという人のArnold Research Service (ARS)を買いMRCAと名付けた。MRCAは30-40年代の市場調査会社の最大手のひとつとなる。さて、Whiteは会社を娘に継がせようとしたんだけど、結局は1951年、Samuel Bartonという人が経営するIndustrial Surveys Companyに売る。合併後、MRCAという社名が継承された。... えーと、ってことはモルゲンシュテルンは51年までには会社を手放していたんだろうね。

 10/23追記: 勢い余ってPauline Arnoldの伝記(Jones, 2013) にも目を通してしまった。
 Paulineは1894年生まれ、セント・ルイスの高校の音楽教師だったが、1926年(ってことは御年30才)、ARSを創業。なぜ音楽の先生がそんなことをはじめたのかはよく分からない由。ARSははじめて地域ごとのフィールド・スタッフを全米規模で組織した会社であり、1936年の時点で1000都市に3000人のインタビュアーを擁していた[←戦前にそんな会社があったのね...]。またARSはcoincidental radio survey、すなわち「いまラジオのCMを聴いてました?」と訊く電話調査のパイオニアであった。それまでラジオ広告への接触調査はあとから思い出してもらうやり方で行われていた。
 Paulineは1928年にNYのマーケティング・コンサルタントPercival Whiteに出会い、1931年に結婚(あ、なるほど...)、NYに移る。ARSはPercival White Inc.と合併してMRCAとなる。夫が社長で妻が副社長。
 さて、夫は次第にMRCAの経営に関心をなくし、軍需用の機械製造みたいなことをはじめちゃったもんで(このへんもちょっと面白い話なんだけど割愛)、妻が社長、夫の連れ子のMatilda(1911年生)が副社長になる。
 このMatildaさんという人もただのお嬢ではなくて、気まぐれな父に代わって経営に尽力し、確率抽出による調査設計によって1946年にAMAから賞をもらったりする。で、一時は会社を継ぐつもりだったんだけど父が約束をまもらないもんだからうんざりし、1949年に大学に移る。両親は結局1951年にMRCAを手放した、という次第。
 Matildaさんはのちに社会学の教授となる。Matilda White Riley (2004年没)、老年学のパイオニアとして知られる研究者だそうだ。

雑記 - 消費者パネルデータの栄枯盛衰

2019年10月18日 (金)

Ferrari, P.A., Barbiero, A. (2012) Simulating Ordinal Data. Multivariate Behavioral Research, 47, 566-589.

 いま、$m$個の順序変数があって、各変数の周辺分布はわかっている。さらに、任意の2個の変数について、(カテゴリに1,2,...と数字を振って得た)相関係数もわかっている。この周辺分布と相関行列にあわせてデータを生成するにはどうすればいいか、という論文。
 なんでそんなの読んでるの?という感じだが、仕事の都合でそういうのが急遽必要になったのである。ああ流されていく人生。

 著者らはRのGenOrdパッケージの中の人。
 イントロは端折るけど、一般に、欲しい特徴を持つデータ行列をシミュレーションする方法は二つあるそうだ。

  • 得られるデータが欲しい特徴を持つようにあらかじめ調整した変数からサンプリングする。
  • サンプリングしてから事後的に調整し、欲しい特徴を持たせる。

 この論文は前者のアプローチ。

 いま、ターゲットである相関行列$\mathbf{R}^{O*}$をつかって、連続変数のデータ行列 $\mathbf{Z} \sim N(\mathbf{0}, \mathbf{R}^{O*})$を生成したとしよう。で、$\mathbf{Z}$の各列について、ターゲットである周辺分布に一致するように離散化したとしよう。ここまでは、まあ、思いつきますね。
 ところが、このデータ行列から相関行列$\mathbf{R}^O$を得ると、それはターゲットの$\mathbf{R}^{O*}$からずれてしまう。そりゃまあそうだ。

 そこで著者らが提案する手順は以下の通り。

  • $\mathbf{R}^{C(0)} = \mathbf{R}^{O*}$を使って$\mathbf{Z}$を生成し、これを離散化して$\mathbf{X}^{(1)}$を得る。
  • $\mathbf{X}^{(1)}$の相関行列$\mathbf{R}^{O(1)}$を求める。
  • $\mathbf{R}^{C(0)}$の要素$\rho_{ij}^{C(0)}$を次式で更新する。
     $\rho_{ij}^{C(t)} = \rho_{ij}^{C(t-1)} f_{ij}(t)$
     $f_{ij}(t) = \rho_{ij}^{O*} / \rho_{ij}^{O(t)}$
    [ああそうか、出てきた$\mathbf{R}^{O(1)}$が目標より大きい(小さい)箇所について、種である$\mathbf{R}^{C(1)}$の該当箇所をちょっと割り引こう(割り増そう)ってことね]
  • 得られた$\mathbf{R}^{C(1)}$がもはや相関行列になっていなかったら、一番近い正則行列を探す。[RではMatrix::nearPD()というのがそれらしい。知らなかったあああ]
  • $\mathbf{R}^{C(1)}$を使って$\mathbf{Z}$を生成し、これを離散化して$\mathbf{X}^{(2)}$を得る。
  • $\mathbf{X}^{(2)}$の相関行列$\mathbf{R}^{O(2)}$を求める。
  • また$\mathbf{R}^{C(t)}$を更新して...という処理を繰り返す。$\rho_{ij}^{C(t)} - \rho_{ij}^{O*}$の最大値が閾値(0.001とか)以下になるか、最大回数を超えたらやめる。
  • 最後に手に入った$\mathbf{R}^{C(t)}$から$\mathbf{Z}$を生成し、これを離散化して出来上がり。

 後半は応用例。この方法でなにかの信頼区間を出すとか、PCAと非線形PCAの性質のちがいをシミュレーションするとか。眠くて読み飛ばした。

 ...意外に単純な話で拍子抜け。要はほしい相関行列を共分散としてもつMVNからサンプリングして離散化し、その相関行列をみて共分散を修正してまたサンプリングして離散化して、というのを繰り返すわけだ。

論文:データ解析(2018-) - 読了: Ferrari & Barbiero (2012) 指定された周辺分布と相関行列を持つ多変量順序データを生成する方法

2019年10月15日 (火)

Hardie, B.G.S., Fader, P.S., Wisniewski, M. (1998) An empirical comparison of new product trial forecasting models. Journal of Forecasting, 17, 219-229.

 消費財の新製品トライアル売上を説明するモデルをいろいろ集め、実データにあてはめて比べてみました、という論文。
 新製品売上予測モデルのバトルロイヤル企画としてはM Competitionsが有名だけど、これはそのトライアル版。まあ誰でも思いつきそうな研究ではあるが(すいません)、レビューのところが勉強になりそうだと思ってめくった。

 イントロ [...省略。トライアル売上の予測がいかに大事かというような、まあそりゃそうでしょという話。直接の先行研究としてMahajan, Muller, & Sharma (1984 MktgSci)を挙げている。あれは認知率の予測じゃなかったっけか。ともかく、読まなきゃいけないなあ]

 CPG(consumer packaged goods)のトライアル予測モデルのうち、パネル/サーベイデータを使い、マーケティング意思決定変数(価格とか広告とか)を使わない8つのモデルを選んだ。
 それでは選手登場です。拍手でお迎えください!

1. 非トライアル者つき指数モデル
 Fourt & Woodlock(1960)は、累積トライアルの天井を$x$, 浸透速度を$r$, 上市からの離散時間を$i$として、時点$i$のトライアル増分を$r x (1-r)^{i-1}$とした。
 連続時間で定式化しよう。ランダムに選んだパネリストについて、その人がトライアルするまでの時間が指数分布 $f(t) = \theta \exp(-\theta t)$に従うとする(つまりある人の非トライアル下での瞬間トライアル率は定数$\theta$)。で、$\theta$の分布は確率$p$で$g(\theta)=\lambda$, 確率$1-p$で$g(\theta)=0$としよう。このとき、累積トライアル率は
 $P(t) = p (1- \exp(-\lambda t))$
 Fourt-Woodlockいわく、このモデルでは累積トライアルの曲線が比較的すぐにフラットになってしまうけど、現実のデータはそうでもない。購入頻度に異質性があるからだろう。

2. 非トライアル者・ストレッチ因子つき指数モデル
 というわけで、Fourt-Woodlock、Eskin(1973)、Kalwani & Silk (1980 JMR) らはこう拡張した:
 $P(t) = p (1- \exp(-\lambda t)) + \delta t$
[あまりに単純な拡張でびっくり。なんすか$\delta$って... 実質的な解釈は難しいよね]

3. 指数ガンマモデル
 いっぽうAnscombe (1961 JASA)は異質性を直接モデル化することを考えた。ガンマ混合分布を考える。
 $g(\theta | r, \alpha) = \frac{\alpha^r \theta^{r-1} exp(-\alpha \theta)}{\Gamma(r)}$
 ここから
 $P(t) = 1 - (\frac{\alpha}{\alpha + t})^r$
これはNBDモデルの待ち時間版だと考えてもよい。

4. 非トライアル者つき指数ガンマモデル
上のモデルで、$\theta$の分布がすごく左に寄っているというかたちで非トライアル者の存在を表現することもできるんだけど、そうすると$\theta$の平均や分散を解釈しにくくなるので、明示的に表現して
 $P(t) = p \left( 1 - (\frac{\alpha}{\alpha + t})^r \right)$
なお、これはもともと Kalwani & Silk (1980) がリピート購買に当てはめたモデルであった。

5. 非トライアル者つきワイブル・ガンマモデル
 Massyという人が60年代に提案した STEAM(Stochastic evolutionaly adoption model)というのがある。トライアルと初回リピートとその後のリピートについてそれぞれモデル化するという複雑な提案である。
 トライアルについていうと、瞬間トライアル率を(定数じゃなくて)時間の関数と見ている。つまり購買間隔はワイブル分布になる。異質性はガンマ分布で表現するので、ワイブル分布のガンマ混合ということになる。ほかにいくつかの仮定を付け加えると下式になる。$c=1$なら非トライアル者つき指数ガンマモデルとなる。
 $P(t) = p \left( 1- \left[\frac{\alpha c}{(t+1)^c + \alpha c - 1}\right]^r \right)$

6. 対数正規-対数正規モデル
 Lawrenceという人は、購買間隔が(指数分布じゃなくて)対数正規分布に従うと仮定した。でもって、その平均が消費者間で対数正規分布すると仮定すると、$\Lambda(t|\mu, \sigma^2)$を$N(\mu, \sigma^2)$の対数として
 $P(t) = \frac{1}{\exp(\mu+\sigma^2/2)}(1-\Lambda(t|\mu, \sigma^2)) + \Lambda(t|\mu+\sigma^2, \sigma^2)$
となる。[ほへぇぇ... 恥ずかしながら初めて聞いたぞ。Lawrence(1979 Euro.Res.; 1982 Euro.Res.; 1985 Mktg Intelligence & Plannning)というのが挙げられている]

7.「二重指数」モデル
 Greeneという人の提案。
 $P(t) = \frac{p}{\beta - \alpha} (\beta (1-\exp(-\alpha t) - \alpha (1-exp(-\beta t))$
 これは購買行動についての仮定からではなくて、単なる曲線近似。Greeneさんは累積トライアル曲線が経験的にS字型だと考えたのでこういう提案をした。

8. Bassモデル
 いよいよ本命、Bass(1969)の登場である。ひかえおろう。
 $P(t) = p \left( \frac{1-\exp(-(\alpha+\beta)t}{1+(\beta/\alpha) \exp(-(\alpha+\beta)t)} \right)$
$\beta =0$なら非トライアル者つき指数モデルになる。
 超有名なモデルではあるが、CPGへの適用は実はきわめてレア。

 選手紹介は以上だが、Bassモデルを除いて、普及モデルなら考えるであろう消費者間の影響を考えていない点にご注目。CPGだからね。トイレットペーパーについてクチコミしないからね。
 ここで浮かぶ疑問は、WoMがないんなら、なぜに累積トライアル曲線がS字型になることがあるのか、という点である(モデル5,6,7,8はS字型を表現しうる)。ひとつの説明は、実際の市場はテストマーケットと異なり配荷がコントロールされていない、S字型になるのは配荷のせいで行動のせいじゃないんじゃない、というもの。
 ついでに予選落ちした諸君をご紹介しよう。Massy et al.(1970 書籍)はさまざまなモデルを提案している(ロジスティックモデルとか線形学習モデルとか)。既存製品の購買行動のモデルを新製品に適用するという手もある[Aaker(1971 MgmtSci.), Herniter(1971 MgmtSci.)というのが挙げられている]。

 IRIのスキャンデータを使う。このデータでは新製品の配荷は常に100%である[どういう仕組みなのかねえ...]。4カテゴリ、19新製品。
 それぞれの新製品について、52週間のトライアル購入者がわかる。前半の26週(ないし13週)でモデルを推定し、残りをあてにいく。
 以下、週$t$の実浸透率を$Y(t)$, その推定値を$P(t)$, 前週からの増分を$y(t)$と$p(t)$, 世帯数を$N$, 推定に使った期間の週数を$T$と書く。
 推定は3種類。

  • MLE。尤度は
     $L = (1-P(t))^{N(1-Y(T))} \prod_{t=1}^T p(t)^{Ny(t)}$
    となりますね。モデル2は推定できない。
  • 増分の差の平方和 $\sum_t^T (y(t)-p(t))^2$を最小化するNLS。
  • 累積の差の平方和$\sum_t^T (Y(t) - P(t))^2$を最小化するNLS。

 指標として、テスト期間の累積トライアルのMAPE, ないし52週目の累積トライアルのAPEをみる。

 結果。
 モデル1,2,3,4の成績がよい。つまり、単純にconcaveなモデル(S字型を表現できないモデル)の勝ち。
 26週で学習した場合、優勝はモデル2。ただし13週の場合、成績ががくっと落ちる。13週での優勝は、なんとモデル1。要するに、ストレッチ因子$\delta$は学習期間がある程度長いときに効いてくるわけだ。
 モデル3,4は優勝は逃したが、26週でも13週でもそこそこの成績を示す。
 下位集団のなかで検討しているのはモデル5で、平均すると5位だが、実は26週でも13週でも4位につけている。その次がBassモデルと「二重指数」モデル。堂々の最下位は対数正規-対数正規モデルであった。
 推定方法をみると、MLEと累積NLSが良い。学習期間が長いとMLEが優位となる。もっともモデルとの交互作用もあって、model 3やBassモデルはMLEがよい。
 52週目のAPEを指標としても、だいたい似たような結果。

 浸透限界$p$を想定する6モデルについて$\hat{p}$や$P(52)$を$Y(52)$と比べると、いずれも過小評価が起きやすい。いっぽう残りの2モデルでは過大評価が起きやすい。
 このように浸透限界の推定は難しい。Van den Bulte & Lilien (1997 MktgSci.)は浸透限界はなにか別の方法で推定したほうがいいと述べているが[←へー]、それは耐久財の話であって、CPGでは難しいね。

 結論と考察。

  • 単純なモデルがよろしい。
  • 異質性を明示的にモデル化するのがよろしい。
  • 非トライアル者を考慮する意義は小さい。モデルに入れると過小評価、いれないと過大評価につながる。
  • 推定方法はMLEがおすすめ。学習期間が短い場合は累積NLSもあり。
  • 学習期間の長さは大事。
  • 学習期間におけるモデル適合にはあまり意味がない。

 最近の実務家や研究者は、新製品トライアルモデルはどのモデルでもロバストだしモデル間で大差ないと思っているかもしれないが、そうではない。モデル選択は大事です。
 云々。

 ... なかなか面白かった。
 この研究の最大の知見は「シンプルなモデルがよい」だと思うのだが、これは要するに、「この研究で使ったデータでは、トライアル売上はどの新製品でも上に凸な曲線だった」ということを意味しているのではないかと思う。
 とすると、この知見はどこまで一般化できるだろうか。この論文で使ったデータはジュース, クッキー、スナック、ドレッシングのトライアル売上で、いずれもクチコミは効きそうにないが、消費財にだってクチコミが効くカテゴリはあるのではなかろうか(化粧品とか?)。また、この論文の時代と比べて小売の商品管理はもっと発達しているだろう。コンビニのように、商品の回転率をシビアに観察して棚をひんぱんに入れ替えている業態を考えると、売上曲線の初期の立ち上がりがのちの配荷に影響し、S字型のトライアル売上曲線が生じたりしないかしらん。

 どうでもいいけど、論文の最後に「本論文ではなんら複雑なソフトは使ってません。全部Excelのソルバーでやりました。まあ時間はかかったけどな」と書いていて、笑ってしまった。そこを自慢するのね。

論文:マーケティング - 読了:Hardie, Fader & Wisniewski (1998) 消費財トライアル売上予測バトルロイヤル

2019年10月10日 (木)

 仕事の都合で読んでんだけど、それにしても、なぜ俺こんなん読んでるんだろう... こういう話は人間にとって心底どうでもよいことだと思ったから人文学部に入ったのにさ... わけがわからないよ...

 Reibstein, D.J., Farris, P.W. (1995) Market share and distribution: A generalization, a speculation, and some implication. Marketing Science, 14(3), G190-202.
 この論文、いきなり「一般化」という節から始まって面食らうんだけど、この雑誌の"Empirical Generalizations in Marketing"という特集号(?)に掲載されたようなので、おそらくは冒頭で当該トピックに関する経験的知見の一般化をまとめるという縛りかなにかがあったのだろう。

 一般化。

  • ブランドシェアと流通配荷のクロスセクショナルな関係は下に凸(convex)なパターンを示す。つまり、シェアの高いブランドのほうが、配荷率1ポイントあたりのシェア増大が大きい。
  • シェアの低いブランドは、浸透率が低くリピート率が低い(double jeopardy)。
  • シェアの低いブランドは広告弾力性が低い。
  • シェアの高いブランドはプロモーションの交差弾力性が高い(シェアの低いブランドがプロモーションをやってもシェアの高いブランドのシェアを奪うのは難しい)。
  • シェアの高いブランドは顧客満足が低い。
  • 国レベルで広告するブランドは価格弾力性が低い。

 ここからは、シェアと配荷のconvexな関係をどう説明するかという理論的な話なんだけど、その前に、配荷の指標について。物理配荷率, ACV, PCVについて説明。ACVはブランド扱店の日用品売上を日用品売上総額で割ったもの、PCVはブランド扱店のカテゴリ売上をカテゴリ売上総額で割ったもの(実務家はこの2つをあまり区別していない)。

 話を戻して、理論的説明。
 配荷率とシェアの間の因果関係について。

  • 配荷率→シェア。(1)配荷率が低いと、買いたい人が買えないから、シェアは低くなる(特にロイヤルティが低い時)。これは線形な関係。(2)配荷率が高いと、他のブランドを買いたかったのに買えなかった人がそれを買うから、シェアは高くなる。これはconvexな関係。(3)小売業者はブランドロイヤルユーザを集客したいから、配荷率が高いブランドについて店舗内プロモーションをしたり割引したりするので、シェアが高くなる。
  • シェア→配荷率。(1)小売業者は棚割をシェアに合わせるので、シェアが高いブランドを選ぶので配荷率が高くなる。(2)棚が限られている小売店では首位ブランドしか置かないので、シェアが高いブランドの配荷率が高くなる。
  • 共通の原因。広告やプロモーションなどは、配荷率とシェアの両方を高める。

 関係の形状について。
 最初に次の2点に注意しよう。(1)配荷がゼロならシェアはゼロだし、シェアがゼロなら(棚落ちするから)配荷はゼロになる。(2)シェアはPCVを超えない。
 さて、形状として以下がありうる。

  • 線形で切片ゼロ。シェア=配荷率x扱店におけるシェア、と考えるとそうなる。
  • 上に凸で切片ゼロ。つまり、扱い店舗が増えることによる売上増大が逓減する。メーカーが自社ブランドのポテンシャルが高そうな店を選んで配荷している場合、店舗に対するサポートが大事な場合、探索ロイヤルティの高いブランドの場合。
  • S字型で切片ゼロ。入手可能性がブランドの魅力につながり、かつ一定のロイヤル顧客がいる場合など。
  • 下に凸で切片ゼロ。市場の構造のせいで、クロスセクションでみてそうなるんだという説明と(次項で述べる)、いや時系列的にそうなるんだという説明がある(こっちはいまのところ証拠がない)。

 クロスセクショナルにみて下に凸になる理由。Nuttall(1965)はこう説明している。
 横軸にPCV、縦軸にシェアをとると、各ブランドにおけるシェアの増大は配荷率増大に伴い逓減する(上に凸)。ただしその曲線の高さはブランドによって違う。
 いっぽう小売店の「商品扱いルール」について考えると、シェア増大に伴いPCVは増えるけど、その増大も逓減する(上に凸)。縦横を入れ替えると下に凸になる。
 あるブランドからみたPCV-シェア関係は上に凸だけど、PCVは商品扱いルールで決まるから、結局そのブランドのPCV-シェア曲線と小売からみたシェア-PCV曲線の交点が均衡解になる。各ブランドの均衡解は小売からみたシェア-PCV曲線の上に乗るから、クロスセクションでみるとPCVとシェアの関係は下に凸になる。
 [ちょ、ちょっとまって... 小売の「商品扱いルール」を表す曲線が、シェアを横にとったときに上に凸になるのはなぜ? コンビニみたいに棚が小さい場合、この曲線は上に凸じゃなくてS字型になるのではなかろうか?]

 実証研究... [いくつか紹介されているが、観察研究だった。せめて道具変数とか使った分析があるかと思ったんだけど、それもなさそう]

 インプリケーション。
 実務的には... 配荷とシェアは相互に影響しており、効果は非線形で推定は困難。顧客ロイヤルティが鍵になる(ロイヤルティが高いブランドは、小売にとっては扱うべきで、メーカーにとっては配荷を増やしてもシェアへの効果が小さい)。
 学術的には... 流通チャネルの構造、棚管理意思決定、長期的効果についての研究が必要。そもそもマーケティング・ミックスの中でも流通の研究は乏しいほうだ。
 この論文は消費財の話だったけど他の種類の財についても研究すべき。
 云々。

 ... 仕事の都合で配荷について考える必要に迫られ、なにがなんだか分かんなくなって混乱していたんだけど、Nuttallという人による説明を読んでずいぶんすっきりした(しっかし、1965年の研究なのね...)。
 ことマーケティング・サイエンスに関していえば、仕事の関連で知りたい話は、こういうちょっと古めの論文を探したほうが効率良く探せるような気がする。もっともこれは私の仕事の性質のせいで、たとえばデジタルマーケティング系の華やかなお仕事をなさっていたら、それはもう全然ちがうのでありましょう。

論文:マーケティング - 読了:Reibstein & Farris (1995) 市場シェアと配荷の関係についてこれまでにわかっていること・いないこと

<< 読了:Wilbur & Farris (2014) 配荷率増大が売上に与える効果はシェアが大きい製品において大きい、SKUレベルでみてもやっぱりそうだった
 
validate this page / CSS