ネスティッドロジットモデルとIIA特性

離散選択モデルに関するこのブログ内の記事はこちら
[タグ : DiscreteChoiceModel]

以前、ロジットモデルとログサム変数についてのエントリーを書きました(その1その2)。このエントリーではネスティッドロジットモデル (入れ子型ロジットモデル, Nested logit model, NLM)について書いてみようと思います。NLM について考えるためには「そもそもなぜロジットモデルじゃいけないの?」ということを理解しなければいけませんが、これは巷でも言われているとおり、IIA 特性というロジットモデルの問題点について理解する必要があります。

選択肢の類似性

例題を使って考えてみましょう。

ある小さな町のオフィス街には牛丼屋しかありません。ランチの選択肢としては「牛丼」「弁当持参」の二つの選択肢しかなく、半分の人は牛丼を、残りの半分の人は弁当を持参していました。そんなある日、この小さな町のオフィス街に新しい牛丼屋が誕生しました。働いている人たちは「別の店作ってくれよ」と思いましたが、こればかりはどうしようもありません。さて、新しい牛丼屋ができた後のランチの選択の割合はどのように変化したでしょうか?

さて、常識的に考えると、この問題の答えは弁当持参の人があいかわらず半数近くで、「牛丼」というカテゴリーでレッドオーシャン(激しい競争)になる、と考えるのが妥当でしょう。典型的な回答は

牛丼屋A:25%, 牛丼屋B:25%, 弁当持参:50%

となるでしょう。新しい牛丼屋の味が気に入った弁当持参の人たちがいるとするならば、多少は牛丼に流れるので

牛丼屋A:26%, 牛丼屋B:26%, 弁当持参:48%

という回答もアリでしょう。選択肢が三つになったのだから

牛丼屋A:33.3%, 牛丼屋B:33.3%, 弁当持参:33.3%

という回答も考えられますが、普通、この回答はどこかおかしい、と考えるでしょう。これがおかしい、と思うのは上の状況では私たちは「選択肢の類似性」を暗黙的に仮定していたからなのです。ようするに「今まで弁当持参していた人は牛丼屋がもうひとつできたぐらいで行動を変えたりしない」ということですね。

“ネスティッドロジットモデルとIIA特性”の続きを読む

【C++】 基本型・クラス間の型変換を定義するには?

c++ では static_cast<double>(int) とかやると型を変換できますが、それを自分で定義したクラスに対して行う方法について書いてみます。
コンストラクタを使ってやる方法もありますが、この方法だと int, double などの基本型(組み込み型)には使えません。たとえばあるクラス X を int 型にキャストしたい場合などですね。

X::operator int() を定義する

こういうときは次のような方法が有効です。

struct X {
  operator int(){ /* 変換のためのコード */}
};

このような operator T を定義しておくと、

X x;
int y=x;
int z=static_cast<int>(x);

というような書き方ができるようになります。

サンプル

上で書いてあることとほとんど同じですが、コンパイルできるコードを一通り書いておきます。
さまざまな型から Hoge クラスにキャストしたり、また、Hoge クラスから別の型へのキャストをするようなサンプルになっています。

#include &lt;iostream&gt;
#define castmsg(from, to) std::cout << "Cast from "#from" to "#to << std::endl;
struct Fuga {
  int fuga;
};
struct Hoge {
  int hoge;
  operator int(){   // Hoge -> int へのキャスト
    castmsg(Hoge,int);
    return hoge;
  }
  operator Fuga(){  // Hoge -> Fuga へのキャスト
    castmsg(Hoge,Fuga);
    Fuga f;
    f.fuga = hoge;
    return f;
  }
  Hoge(int h){      //  int -> Hoge へのキャスト
    castmsg(int,Hoge);
    hoge=h;
  }
  Hoge() : hoge(0) { }  // デフォルトコンストラクタ
};

int main(){
  Hoge hoge;
  hoge.hoge = 1;
  std::cout << "=== Type cast test (by static_cast) ===" << std::endl;
  /* Hoge -> int */
  int h = static_cast<int>(hoge);
  /* Hoge -> Fuga */
  Fuga fuga = static_cast<Fuga>(hoge);
  /* int -> Hoge */
  Hoge hogera = static_cast<Hoge>(4);
  std::cout << "Hoge->int  : " << h << std::endl;
  std::cout << "Hoge->Fuga : " << fuga.fuga << std::endl;
  std::cout << "int->Hoge  : " << hogera.hoge << std::endl;
  std::cout << "=== Type cast test (implicit) ===" << std::endl;
  /* Hoge -> Fuga */
  fuga = hoge;
  /* int -> Hoge */
  hogera = 6;
  std::cout << "hogera: " << hogera << std::endl;    // ここでも暗黙の Hoge->int が使われる
  //std::cout << fuga << std::endl;    // これはエラー!(Fuga->int は存在しない)
  std::cout << "fuga.fuga: " << fuga.fuga << std::endl;
}

/* 出力 */
/*
=== Type cast test (by static_cast) ===
Cast from Hoge to int
Cast from Hoge to Fuga
Cast from int to Hoge
Hoge->int  : 1
Hoge->Fuga : 1
int->Hoge  : 4
=== Type cast test (implicit) ===
Cast from Hoge to Fuga
Cast from int to Hoge
hogera: Cast from Hoge to int
6
fuga.fuga: 1
*/

暗黙の型変換を抑制する

上の例で言うと、ユーザー定義クラスの Hoge に int が = 演算子を使って代入できるようになっていますが、場合によってはこの挙動は危険だ、という場合もあるでしょう。そんなときはコンストラクタ Hoge(int) の前に explict 修飾子をつけましょう。
これをつけることで

Hoge x=1;

のようなコードをコンパイル時エラーにすることができるようになります。

[R] RSQLite の使い方 (2)

前回は R で SQLite を利用するための基本中の基本を示しました。
今回は高速化の方法など少しだけ踏み込んだ内容について記します。

トランザクションの利用

トランザクションを使って大量のクエリを高速に処理します。トランザクションって何?って人はググってみてください。前回も書きましたが、dbWriteTable などはトランザクションを利用しないので、場合によってはこちらを使ったほうが速い、のかもしれません(未確認)。

library(RSQLite)
driver=dbDriver("SQLite")
dbname="test.db"
con=dbConnect(driver,dbname)

# 準備
test.tbl=as.data.frame(matrix(rnorm(100),nc=2))
colnames(test.tbl)=c("random1","random2")
rs=dbWriteTable(con,"random",test.tbl,row.names=F)

# トランザクション開始
dbBeginTransaction(con)
rs=dbSendQuery("DELETE FROM random WHERE random1>0 AND random2>0;")

# DELETE される行が 10 行より多ければロールバックする(DELETEしない)
if(dbGetInfo(rs)$rowsAffected > 10){
  cat("Rollback!")
  dbRollback(con)
} else { # そうでなければコミットする(DELETE が反映される)
  cat("Commit!")
  dbCommit(con)
}

# コミットされたかどうかチェック
#(この場合コミットされていないのでデータベースの内容はそのまま)
dbGetQuery(con,"select count(*) from random;")
#=>   count(*)
#=> 1       50

C/C++-API でいうところの sqlite3_prepare & sqlite3_bind_XX の利用

これも高速化のテクニックのひとつです。あらかじめパラメータつきのクエリを準備しておいてパラメータに値を次々に代入して SQL 文を実行するための方法です。

library(RSQLite)
driver=dbDriver("SQLite")
dbname="test.db"
con=dbConnect(driver,dbname)

# 準備(上で作った random テーブルがすでにあれば必要なし)
test.tbl=as.data.frame(matrix(rnorm(100),nc=2))
colnames(test.tbl)=c("random1","random2")
rs=dbWriteTable(con,"random",test.tbl,row.names=F)

# Prepare された SQL ステートメントを使ってデータ挿入
dbBeginTransaction(con)
keys=data.frame(r1=runif(10), r2=runif(10)) # 適当なデータ

# セミコロンに続いてデータフレームのカラム名で指定
pquery="INSERT INTO random (random1, random2) VALUES (:r1,:r2);"
rs=dbSendPreparedQuery(con, pquery, keys)
dbCommit(con)

拡張関数の利用

RSQLite.extfuns パッケージを使うと sin, exp, floor のような数学関数,trim のような文字列関数,stdev, mode, median のような統計関数が使えるようになります。
これらの関数を使うことでより複雑なクエリを生成することが可能になり、結果として処理速度の向上が見込めるかもしれません。

下の例では有名なアヤメデータ(iris)に対して、

  • 種(Species)ごとにグループ分けして がく片の長さ(Sepal.Length)の標準偏差の小さい順にデータフレームを作成

という処理をしています(注:plyr パッケージの ddply 関数を使えばほぼ同じことができます)。

library(RSQLite)
library(RSQLite.extfuns)   # 拡張関数用のパッケージ
driver=dbDriver("SQLite")
dbname="test.db"
con=dbConnect(driver,dbname)
init_extensions(con)       # 拡張関数が使えるようにするためのおまじない

# 準備(iris データを使う)
data(iris)
rs=dbWriteTable(con,"iris",iris,row.names=F)

# Species ごとにグループ分けして Sepal.Length の標準偏差の小さい順にデータフレームを作成
dbGetQuery(con,"select Species, stdev(Sepal_Length) as slstdev from iris group by Species order by slstdev;")

[R] RSQLite の使い方 (1)

RSQLite(+RSQLite.extfuns)パッケージを使うことで R から SQLite3 へ非常に簡単にアクセス可能になる。SQLite 本体もどうやらパッケージに含まれているようで、事前にインストールする必要がないため、パッケージをインストールするだけで非常に簡単に使い始めることができます(パッケージ自体も 1MB もないくらいのサイズです。如何に SQLite が小さいかがよくわかります)。

クエリーを投げる、トランザクション、などの基本的機能は完備されています。さらにデータフレームを一気に INSERT するための dbWriteTable などのおかげで便利に使えます。以下に使い方を示していきます。

一番シンプルな使い方

# (0) 準備:ドライバに SQLite を指定し,データベースをオープン.
#           コネクション・オブジェクト (S4) が帰ってくるので
#           保持しておく.
dbname="test.db"
library(RSQLite)
driver=dbDriver("SQLite")
con=dbConnect(driver,dbname)

# (1) dbGetQuery : すでに存在するデータベースを open してテーブルを取得
# (hoge というテーブルを持っていると仮定)
tbl=dbGetQuery(con,"SELECT * from hoge;")

# (2) dbSendQuery : テーブルが帰ってこない(SELECT 以外の?)クエリーは
#     dbSendQuery を使うべき → エラー処理が可能に
#     (dbGetQuery でも NULL が返るだけでエラーにはならない)
rs=dbSendQuery(con,"CREATE TABLE new_tbl(id INTEGER, hoge TEXT);")
dbinfo=dbGetInfo(rs)
print(dbinfo)

# dbinfo の中身
#=> $statement
#=> [1] "CREATE TABLE new_tbl(id INTEGER, hoge TEXT);"
#=> $isSelect
#=> [1] 0
#=> $rowsAffected
#=> [1] 0
#=> $rowCount
#=> [1] 0
#=> $completed # 正常に終了したか?
#=> [1] 1
#=> $fieldDescription
#=> $fieldDescription[[1]]
#=> NULL

データフレームの取得、挿入

R のデータフレームはリレーショナルデータベースと似た構造になっているので、write.table, read.table のようなノリで読み込み、書き込みを行うことができます。

library(RSQLite)
driver=dbDriver("SQLite")
dbname="test.db"
con=dbConnect(driver,dbname)
test.tbl=as.data.frame(matrix(rnorm(100),nc=2))
colnames(test.tbl)=c("random1","random2")

# データフレームをテーブルとして一括書き込み
rs=dbWriteTable(con,"random",test.tbl,row.names=F)

# 書き込めているかテスト
dbGetQuery(con,"SELECT * FROM random limit 5;")

#=>   row_names    random1    random2
#=> 1         1  0.2354093  1.7083581
#=> 2         2 -0.8101699 -0.6867071
#=> 3         3  0.4537572  1.8399874
#=> 4         4 -2.1129855  1.0211205
#=> 5         5  0.9778609 -0.9177596

# テーブルをデータフレームとして一括取得
dbReadTable(con,"random")

(注意) マニュアルによれば dbReadTable, dbWriteTable はトランザクションを利用しないようです。

続きを書きました。

ランダム多項式

以前、ランダム行列について書いたことがありました(その1その2)。

ランダム行列ときたらランダム多項式だろう、ということで数値実験して遊んでみました。
ここでランダム多項式といっているのは
r_0+r_1x+r_2x^2+\cdots+r_nx^nで係数 \(r_i\) がある確率分布に従うようなもののことを言っています。係数は複素数なものを考えます。

ここでは係数の分布は実部と虚部ともに独立な標準正規分布の多項式と、独立な指数分布の多項式の二種類を考えることにしましょう。

ランダム多項式の根を R の関数 polyroot で求めてプロットすると以下のようになります(ソースコードはこのエントリーの一番下に記しておきます)。
前者が標準正規分布で、後者が指数分布のものとなります。
randompoly_normal.pngrandompoly_expon.png
なんとも不思議なことに両方共の例でほとんどの根が半径1の円の上に一様に分布しています。(どうやら)この現象にはきちんと理論的な解明がなされているようで、”random polynomial” とかでググるといろいろ出てきます。分布に関するかなりゆるい条件でこのような円則が現れるようです。

ところでこのランダム多項式、どのような状況で使うのでしょうか…

すぐ思いついた(物理的な意味のない)例では、ランダムな係数の線形漸化式の特性方程式。この漸化式はこの結果をみればわかるとおり特性方程式が絶対値が1より大きな根を持つと考えられるので確率1で発散しそうですね。

# Coefs: standard normal distribution
plot(polyroot(rnorm(201)+rnorm(201)*1i),xlim=1.5*c(-1,1),ylim=1.5*c(-1,1),xlab="Re",ylab="Im")
title("Roots of random 200th degree polynomial\n(the coefficients are complex-valued and\n drawn from standard normal distribution)")
# Coefs: exponential distribution
plot(polyroot(rexp(201)+rexp(201)*1i),xlim=1.5*c(-1,1),ylim=1.5*c(-1,1),xlab="Re",ylab="Im")
title("Roots of random 200th degree polynomial\n(the coefficients are complex-valued and\n drawn from exponential distribution)")