The One with ...

思いついたことや作業メモなどを公開しています

Pythonによるジニ係数の計算コード

内容

ジニ係数Pythonを使って計算します.よく知られた定義ではなく,簡略化された定義式を使うと計算が速くなる,という例を紹介します.

ジニ係数の直感的意味

ジニ係数は,所得や資産などの非負の資源分配の不平等度を測定する指標です. たとえば,社会の全員に全く同じ量の富が分配されている場合にはジニ係数Gは0の値をとります.

y=[1,1,1,1,1,1]
G(y)=0

もし,たった1人が, その社会に存在する全ての富を独占している状態だと,ジニ係数は1の値をとります.

y=[1,0,0,0,0,0]
G(y)=0

つまり,ジニ係数は,完全に平等な状態であれば0,不平等が大きくなるほど1に近づく指標です.

ジニ係数の定義

よく知られたジニ係数の定義は 所得ベクトル 
y=(y_1, y_2, \ldots, y_n)
に対して,

\displaystyle{
G=\frac{1}{2 n ^2 \mu}\sum \sum  | y_i - y_j | \qquad  ただし \mu=\frac{y_1+y_2+\cdots+y_n}{n}
}

です(なお \sum \sum_{i=1}^{n}の略です.はてぶろ上では2重和がうまく表示できないので省略しました).この定義を使った場合,最大数は 1-1/nです. 最大値をぴったり1にしたい場合は

\displaystyle{
G=\frac{1}{2 n(n-1) \mu}\sum \sum  | y_i - y_j |
}

を使います.ジニ係数の定義の直感的な意味は,「全てのペアの所得の差の絶対値を合計して,基準化した量」です. これを計算するPythonコードを考えてみましょう.定義通りにコードに翻訳すると,例えば次のような関数になるでしょう.

import statistics
import math

def gini(y):
    m = statistics.mean(y)
    n = len(y)
    a = 2*m*(n*(n-1))
    ysum = 0
    for i in range(n):
            for j in range(n):
                ysum = ysum +(math.fabs(y[i]-y[j]))
    
    return(ysum/a)

ためしに計算してみましょう.

y=[1,0,0,0,0,0]
gini2(y)

Out[ ]: 1.0

y=[1,1,1,1,1,1]
gini2(y)

Out[ ]: 0.0

完全不平等の場合に1,完全平等な場合に0をたしかに出力しています.

y=[1,2,3,4,5,6]
gini2(y)

Out[ ]: 0.3333333333333333

一様分布の場合は,約G=1/3でした.

ところでジニ係数の定義を見て分かるとおり,所得ベクトル要素の差の絶対値を2重に足しています. |x_i-x_j|と|x_j-x_i|は常に同じ値なので 計算としてはずいぶん無駄な作業をしています. 上記の工夫のないコードでも,データが少なければ特に問題は無いのですが,データが増えてくると若干計算速度が気になります.

たとえばrange( )関数を使って,データを生成しながら計算してると,長さ2000くらいの短いリストでも出力までに時間がけっこうかかります.

import time #時間計測用
y=range(2000)    
start = time.time()
[gini(y),time.time() - start]

Out[ ]: [0.33366683341670833, 2.4708046913146973]

なんということだ!2.47秒もかかっている. これではシミュレーションでジニ係数を反復計算する際に面倒です.

ジニ係数の別の定義式

そこで,計算を早めるためにジニ係数の別の定義式を使います. 次の定義式は上述した式と同値であることが証明されています.ただしこの定義を用いる場合は,所得ベクトルを昇順に並び替える必要があります.

  \displaystyle G= \frac{ \sum i y_i}{ \sum y_i}-  \frac{n+1}{n}

ただし  y_1 \le y_2 \le \cdots \le y_n

en.wikipedia.org

この式を使えば総和が2重になっている部分が簡略化されるので,ぐっと計算量が減ります.最大値を1にしたい場合は調整用に n/(n-1)を掛けます.

  \displaystyle G=( \frac{ \sum i y_i}{ \sum y_i}-  \frac{n+1}{n} )  \frac{n}{n-1}

上記の定義にもとづくコードは例えば以下のようなものになるでしょう.

#sortして計算の効率化
def gini2(y):
    y.sort()
    n = len(y)
    nume = 0
    for i in range(n):
        nume = nume + (i+1)*y[i]
        
    deno = n*sum(y)
    return ((2*nume)/deno - (n+1)/(n))*(n/(n-1))

速くなったかどうか,確かめて見ましょう.

y=list(range(2000))
start = time.time()
[gini3(y),time.time() - start]

Out[9]: [0.3336668334167085, 0.0]

速くなりました!

以上です.

ggplot2: 複数の関数をplotして比較する

関数のパラメータを変えながら,そのplotを比較したい,という場面に日々遭遇します. 例えば平均の異なる正規分布の密度関数を比較したい場合,Mathematica(というウルフラム社の数式処理ソフト)なら,

Plot[{
  PDF[NormalDistribution[0, 1], x],
  PDF[NormalDistribution[1, 1], x],
  PDF[NormalDistribution[2, 1], x]}, {x, -6, 6}]

でOKです.

f:id:hamada7418:20190418122206p:plain

凡例が欲しければ,Plot関数のオプションに

PlotLegends -> "Expressions"

と指定します.

f:id:hamada7418:20190418122513p:plain

え?これだけ?というくらい簡単です.

同じことをRでやってみましょう*1

R上での関数のplot

もっとも簡単なコードとして

curve(dnorm(x),-4,4)

を試してみます.

f:id:hamada7418:20190418125100p:plain

ファイルとして出力するには

png("hoge/pics.png")
curve(dnorm(x),-4,4)
dev.off()

です.hogeディレクトリにpics.pngが保存されます *2

複数の関数を1つの図に重ねてplotしましょう.par(new=TRUE)でplotを上に重ねて書きます.

curve(dnorm(x),-4,4)
par(new=TRUE)
curve(dnorm(x,1,1),-4,4)

f:id:hamada7418:20190418131007p:plain

グラフを較べる目的ならこれで十分でしょう. 縦軸ラベルの重なりが気になる場合は片方のラベルを消します.

curve(dnorm(x),-4,4)
par(new=TRUE)
curve(dnorm(x,1,1),-4,4,ylab="")

f:id:hamada7418:20190418142109p:plain

線種を変えて,凡例もつけましょう.

curve(dnorm(x),-4,4,lty=1)
par(new=TRUE)
curve(dnorm(x,1,1),-4,4,lty=2,ylab="")
legend("topright",legend=c("mean=0","mean=1"),lty=c(1,2))

f:id:hamada7418:20190418142448p:plain

見やすくなりました.

ggplot2をつかう

ggplot2を使って,グラフの見た目を調整します.

library("ggplot2")
x <- seq(-4,5,0.1)
m0 <- dnorm(x,0,1)
m1 <- dnorm(x,1,1)
data <- data.frame(x=x,y1=m0,y2=m1)
ggplot(data,aes(x))+
    geom_line(aes(y=y1))+ 
    geom_line(aes(y=y2),linetype="dashed")

f:id:hamada7418:20190418165519p:plain

基本的な方針は以下のとおりです.

  1. 独立変数のベクトルxをつくる
  2. dnorm(x,0,1)の出力をベクトルm0に代入
  3. dnorm(x,1,1)の出力をベクトルm1
  4. 上記のベクトルをdata.frameにまとめる
  5. geom_line( )をで曲線をplotする

以上です.ラベルが不要であれば +theme(axis.title=element_blank( ) ) で消去します.このときxlab=" "等を使うと余分な空白が入るので,axis.title=element_blank( ) )を使うとよいでしょう. 背景のグレーが不要ならtheme_bw( )を使って,背景を白黒印刷用に白くします. 以上をまとめためのが,次のコードです

ggplot(data,aes(x))+
    geom_line(aes(y=y1))+ 
    geom_line(aes(y=y2),linetype="dashed")+
    theme_bw() +theme(axis.title=element_blank(),text=element_text(size=9)) 

f:id:hamada7418:20190418165530p:plain

以上です.もっと簡単な方法が分かれば,メモを更新します.

*1:Stanをつかって決定論的な関数を統計モデルに変換する際に, R上で関数比較をすると作業が効率化される場合が時々あります.やり方をすぐに忘れるのでここに方法をメモしておきます

*2:作画デバイス上のエラーが原因で,作成した画像ファイルを開けないことがあります. その場合はdev.offやgraphics.offで作画デバイスを消します

物語の機能とはなにか?


登場人物

青葉:数学苦手女子

花京院:数学好き男子


文学部をでたけれど

青葉と花京院は駅前の喫茶店で,コーヒーを飲みながら本を読んでいる.

さきほどから青葉は文庫版の小説を,花京院は英語で書かれた論文のコピーを読み続けていた.

テーブルを挟んで向かい合わせに座ったまま,お互いに一言も喋らない.

駅前の喫茶店は相変わらず静かだ.純粋にコーヒーの味でいえば,シアトル系チェーン店のほうが青葉の好みだった.でも,おちついて本を読んだり会話ができるという意味でやはりこの古い店を選んでしまう.そしてなにより,会社帰りにここに立ち寄れば花京院と話ができる.約束をしているわけではないので,運がよければ,の話ではあるが.

やがて青葉は開いていた頁にしおりをはさんでそっと閉じると,小説をテーブルの上に置いた.

「ねえ,花京院くん」

 花京院は手にしていた論文から目を離さずに「なに?」とだけ答えた.テーブルの上には計算用紙が散乱している.おそらく書いた本人でさえ,数ヶ月後には判読できないほど乱雑な筆跡で,数式が書き込まれている.

「わたしたちって文学部の出身でしょ?」

「そうだね」

青葉と花京院は某大学文学部の出身だ.2人は同期で同じ研究室を1年ほど前に卒業した.

「会社の子と話をしていて,ふと思ったのよ.文学部を卒業したものの,結局《文学とはなにか》っていうテーマでぜんぜん勉強をしなかったなーって」

「僕らは文学部の中でも数理行動科学研究室だからね.たしかにそういう授業は受講しなかった」

「文学って結局なんなの?わたし,よく分かんないんだよね.4年間も文学部にいたのに」

「ずいぶんざっくりした質問だなあ.それは僕の専門じゃないよ」

「でも考えたことはあるでしょ?だって花京院くんは,わざわざ理学部から文学部に転部してきたんだから」

文学と論文をつなぐもの

花京院は大学2年生のとき,理学部応用数学科から文学部数理行動科学研究室に転部したという,少し変わった経歴の持ち主だった.

「うーん,文学とは何かを答えることはできないけど,文学や物語の持つ機能なら少し考えたことがあるよ」

「ほおー」

「たとえば,僕がいま読み直しているこの論文」花京院は手に持った英語論文をテーブルの上に置いた.

Gale, David and Lloyd Shapley, 1962, College Admissions and the Stability of Marriage, The American Mathematical Monthly, 69: 9-15.

「これはゲールとシャプレーがはじめてDAアルゴリズムを発表した歴史的な論文だ.わずか6頁で,それまでに誰も考えてこなかった問題を明確に提示して,それを解決するための概念と命題を証明した模範的な論文だよ」

「これ,まえにマッチングの説明をしてくれたときにつかった論文だよね」青葉は花京院が置いた論文を手に取った.論文にはところどころ証明をフォローするために花京院が書き込んだメモが残っている.

「こういうのはね,文学の持つ一部の機能を極限まで研ぎ澄ました形態の作品だといえる」

「ちょっとなに言っているかわからない」

「すべての研究には共通の目的がある.誰も解いたことのない新しい問題を提示して,それを解くという目的だ.科学論文というのはね,研究内容の情報を圧縮して,必要なことだけを,効率良く文章化したものなんだ」

「ふむふむ.まあそりゃそうだね」

「情報を伝えるための文章という意味で,論文も文学の一形態だと僕は考えている.だから一般的な文学作品や小説にあって,論文にないものを考えると,逆に文学の本質が分かると僕は推測した」

「う-ん,あいかわらず花京院くんはヘンなこと考えるなあ」青葉は腕を組んで,目をつむった.

「きみは専門的な理論や研究内容や手法を,物語の形式で紹介した本を知ってる?」花京院が質問した.

「たとえば?」

「たとえば・・・・・・,そうだね,ぱっと思いつくところだと

あたりかな」

「『もしドラ』だけ知ってるよ.本屋でちょっと立ち読みしたことがある.表紙が可愛いんだよね」

「『もしドラ』はピーター・F・ドラッカーの経営哲学を,高校野球を舞台にした物語で紹介する入門書だ.原典であるドラッカー自身の『マネジメント』は日本で翻訳版が100万部売れた.ものすごいヒット作だよ.でも『もしドラ』はそれを上回り,なんと270万部も売れた」

「へえー.すごいねえ.そんなに売れたんだ.『もしドラ』は確かにおもしろいし,読みやすいもんね.買ってないけど」

「この現象は,文学や物語の機能を劇的に示している」

「どういうこと?」青葉は首をひねった.

ドラッカーの教えを伝えるという意味で,情報量は原典の『マネジメント』のほうが明らかに多い.にもかかわらず多くの人々は『もしドラ』の方を手にとった.なぜか?」

「そりゃあ『もしドラ』の方が読みやすいからでしょ」

「その通りだけど,問題はなぜ読みやすいのか?という点だ」花京院は質問を言い換えた.

花京院にそう言われて,青葉は考えた.

(読みやすい理由? そんなの《読みやすい簡単な文章》で書かれているから以外にないじゃん)

「えーっと,ストーリーがおもしろいから?」

「それも大きな理由の一つ.ほかには?」

花京院に促されてすこし考えてみたものの,ストーリー以外の答えを見つけることはできなかった.しかたなく,コーヒーを飲みながら青葉は花京院の説明を待った.

「可愛い女の子がでてくるからだ」

青葉は思わず口に含んだコーヒーを吹き出しそうになった. 「え,ちょっと.それ真面目に言ってるの?」

「真面目だよ.もう少し抽象的に言えば,《キャラクタ》の存在ともいえる」

「キャラクタ・・・・・・」

キャラクタの機能

「『もしドラ』にあって『マネジメント』にないもの,それは物語とキャラクタだ.人々を『もしドラ』や『ソフィーの世界』や『数学ガール』に引き込む力の正体の一つ.それが《キャラクタ》なんだよ」花京院が説明を続ける.

「うーん,でも物語風に説明するんだから,キャラクタが出てくるのは当たり前な気がするんだけど」

「その通り.したがって問題の本質は,キャラクタとはなにか?ということだ」花京院はさらに問い直した.

「なるほど.じゃあキャラクタってなんなの?」青葉は身をのりだした.

「キャラクタの正体とは・・・・・・・.ずばり,人だ」

「はあー」青葉はがっくりと肩をおとした.「だからそんなのあたりまえじん」

「いや,きみはキャラクタという現象の,不思議さをまったく理解していない」そう言って花京院は意味ありげに笑いながら,鞄の中から一冊の漫画をとりだした.

「あ,それ私が貸してあげた『ハンター×ハンター』じゃん」 花京院がとりだしたのは,人気少年漫画の単行本だった.彼は頁をめくって,そこに描かれた登場人物の1人を指さした.

「彼は生きているかな?」花京院が質問した.

「そりゃあ,生きているよ」青葉は,やや戸惑いながら答えた.

「ほんとに?どうしてそう思うの?」

「どうしてと言われても,話したり動いたりするんだから,そりゃ生きているにきまっているじゃん」

「確かにこのキャラクタは生きている.喋って,考えて,戦って,怒って,恐れて,泣いて,僕たちと同じように生きている.ところできみは,人工知能やロボットは,僕たちと同じように心をもって生きていると思う?」花京院はパラパラと漫画の頁をめくった.

「うーん,いまの段階では,まだ心をもつとはいえないんじゃないかな.よく知らないけど」

「じゃあドラえもんは?」

ドラえもんには心があるよ.ちょっとそれ,ずるくない?」青葉は質問に答えながら,なんだか騙されたような気持ちになった.

「ロボットに心があるかどうかを疑う人でも,ドラえもんに心があることは疑わない.なぜならドラえもんは命を持ったキャラクタだからだ.でも見方を変えてみると,キャラクタは紙の上に印刷された絵や文章に過ぎない.物理的には紙の上に印刷されたインクでしかないにもかかわらず,僕らはキャラクタを自分たちと同じように心を持った人として認識している.実際には作者が作り出した人工の創造物であるにもかかわらず,僕らにはその想像上の人物が動いて生きているように見える.われわれにとってあまりにもあたりまえの事実だから,それを疑ったりしないけど,これは錯覚に過ぎない」

青葉は少しずつ花京院の真意を理解していった.

「うーん,まあ確かに花京院くんの言うとおりなんだけどさあ.それをいっちゃあ,お終いっていうか.それは言わないお約束なんじゃないの」

「でも僕らは無理してその《お約束》につきあっているわけじゃない.ここが重要だ.作中に人物が登場し,なんらかの台詞を喋ったら,その登場人物がそのとき考えていることを喋ったという風に必ず錯覚する.この錯覚は不可避だ.ではなぜこんな錯覚を起こしてしまうのだろうか?」

「分からないなー.そもそも,そんなこと,考えたことなかったよ」

《物の理解》と《心の理解》

「ポール・ブルームという心理学者が『神は偶然の産物か』という論文のなかで,この問題と同型の問題を考えている.読んだことある?」

「ないない.初めて聞いたよ」

「おもしろい論文なんだよ.彼は宗教には多様な形態が存在するにも関わらず,多くの人々が魂や死後の世界や奇跡や万物主による世界(宇宙)の創造を信じるのはなぜか,という問題を認知科学の観点から考えたんだ.ようするに人が神様を信じるのはなぜか?と考えたんだよ」

「へえー,そんな論文があるんだ.で,どうして信じちゃうの?」

「ブルームの説の主発点はこうだ.まず人間は

  1. 物理的世界を理解できる
  2. 他者の心を理解(推測)できる

と仮定する.物理的世界の理解というのは,重力によって物体が落下したり,摩擦によって物体の運動が停止したりする現象の物理法則を人間が理解できるという意味だよ.こういう物理法則の理解は比較的知能が高い動物や生後六ヶ月の幼児でも可能であることが知られているんだ」

「ふむふむ.たしかにそうだね」

「次に人は,より高次元の能力として,他者の意図を解する能力を獲得する.これによって人は他者の感情などを類推できるようになるんだ.前者を《物理的世界の理解》,後者を《社会的世界の理解》と呼ぶことにしよう」

「なるほど.ようするに人は,モノの動きと心の動きが理解できるってことだね」

「そのとおり.《物理的世界の理解》と《社会的世界の理解》の二重性は,《物体としての肉体》と,《物体とは別の心という存在》の別々の知覚を可能にする.そしてわれわれの社会的世界の理解が暴走した結果,本来は存在しない意図や欲求があたかも存在するかのように知覚してしまうんだ」

「ちょっとなに言っているかわからない」

「例えば,そうだね・・・・・・.カフカの『変身』や,中島敦の『山月記』は呼んだことある?」花京院が聞いた.花京院の口から文学作品のタイトルが出るのは少し意外な気がした.ただし,彼は語らないだけで意外と文学好きなのかもしれないなと青葉は考えている.

カフカは読んだことないけど,山月記は知ってるよ.現国の教科書に載ってた」

「あれは,意識は人のままで肉体だけが虎に変化してしまう話だ.『変身』も主人公の精神はそのままで,姿だけが虫に変わってしまう.僕らは現実にはそんなことは,ありえないことだと知っているけど,物語としては容易に《こころ》と《からだ》の分離を受け入れることができる.他にも肉体間で精神だけが入れ替わったり,精神だけが肉体を離れる物語の例がたくさんある.筒井康隆の『転校生』や新海誠の『君の名は』,古い例だと『源氏物語』がそうだ」

「たしかに,そういうお話しにたいして,誰も『いや,入れ替わらないでしょ』ってツッコまないよね.でも言われてみれば不思議だよね.現実には肉体から精神が離れたりしないってことは知っているのに,フィクションとしてならそれを自然に受け入れるのって」

「心理学者Jesse BeringとDavid Bjorklundらの実験によれば,多くの子供はワニに食べられてしまったネズミの肉体的機能が停止する事は認めながらも,食べられたネズミの精神的機能(感じたり,考えたり,欲求をもつこと)は消失せずに残っていると認識するそうだ.つまり心だとか人格というものが肉体を離れて存在する現象をわれわれは幼い頃からごく自然に受け入れることができるんだ」

「なるほど.で,その話と文学がどう関係するの?」

「さっき言ったように,われわれが世界を理解する方法には大きく分けて2つある.物理的法則の理解と,他者に心があるという理解だ.おそらくこの能力は進化的適合の観点から有利だったんだろう.物理現象を正しく把握できたり,仲間の意思を正確に理解できる個体は,それができない個体よりも子孫を残すのに有利だったはずだ. 《物理的世界の理解》と《心的世界の理解》という世界理解の二重性が,精神と肉体(物体)の分離を可能にする.多くの場合,これらを混同したりはしないけど,時々間違った対象にこれらの理解を適応することがある」

「間違うってどういうこと?」

「単なる物理的な対象に心を帰属させたり,逆に心がある他者をモノのように扱ってしまうんだ.前者によってキャラクタが生きているという錯覚が起こる.心的世界の理解力は時に暴走して,心や意思を持たぬモノに精神を(誤って)帰属させる.つまり文学を文学たらしめる最も重要な要素であるキャラクタは,ぼくらの《他者に心がある》という理解の過剰な適用によって生まれる錯覚なんだ」

「うーん,そっかー.錯覚なのか.ちょっと信じられないなー」青葉はすぐには花京院の説明を受け入れることができなかった.

制御不能なキャラク

「僕の予想では,キャラクタが生きてている,という錯覚は作者自身にも起きる」

「え?作った本人が錯覚するの?そんなわけないでしょ.自分で作ったことを自分が一番よく知ってるんだから」

「これについては証拠がある.多くの作家や漫画家が自分の作ったキャラクタが作中で勝手に動くと述べている.たとえば

ネームでは岩鬼が三振するシーンを描いたのに、ペンを入れたら岩鬼がホームランを打ってしまった(水島真司)

実生活でも,明日の予定がないとしても,人はそれなりに動いていくのと同じで,漫画の中でも,そこに描かれた状況があれば,キャラクターは自然に動くベくして動くものです

描いている僕自身ですら「しまった・・・・・・,これは,仗助勝てないかも」と思うぐらいの状況に陥ってしまいました. (荒木飛呂彦).

架空の登場人物達がお話の中で会話をしたり,事件を解決したりしている.不思議な感覚だ(森博嗣).

作者が創りあげたキャラクターだからといって,それは必ずしも100%作者の意思によって動くということではなく,あたかも自分が独自の意思を持って行動を始めるかのような部分ができてくる.そうした変貌を,「キャラクターが独り歩きする」というのです(藤子・F・不二雄).

このように,キャラクタはしばしばその作者にとってさえ,独立した実在として認識される.ようするに自分で作っておきながら自分で制御できず,自分とはちがう他者として生きているように錯覚してしまうんだ」

「うーん,なんか不思議.自分で書いているのにキャラクタを自分でコントロールできないってどういうことなんだろ?」

「つまり僕らは必要以上に,対象の心を読み取ってしまうんだろうね.たとえ空想上の人物でさえも.人類学者のPascal Boyer は,人がしばしばありもしない目的や意図や設計を見出すことを社会的認知の肥大化(hypertrophy of social cognition)と呼んでいる」

「どうしてそんなことが起こるんだろ?」

「ヒトという種のなかで他者への共感性の高い個体は,間接互恵性によって,他者からの協力行動を得やすい.その結果,適合度が高くなり自然選択によって繁栄したと考えられる.共感性が優れて高い個体は,その知覚のエラーにより,しばしば存在しない意図や目的を人以外の自然物のなかに読み取ったんじゃないかな.ブルームはそれこそが神の原初形態だと考えている.僕らはその能力を受け継いだ結果,キャラクターが登場する物語に自然に引き込まれる」

「それが文学ってこと?」

「僕はね,あらゆる作品は,小説でも映画でも漫画でも研究論文でも,それが次の3要素から構成されると考えている.

  • 物語
  • キャラク
  • 情報

の3つだ.この3要素の配分によって作品の性質が決まる」花京院はテーブルの上に置いた計算用紙に3つのキーワードを並べた.

「小説とか漫画とかはそうかな.でも論文にはキャラクタは出てこないでしょ」青葉は花京院が示した3要素を眺めながら言った.

「論文のキャラクタは確かに,その存在を知覚することが難しい.でも確かに存在するんだ.なにか分かる?」花京院が聞いた.

「論文に登場するキャラクタ?うーん全然分かんない.大阪大学のワニ博士とか?」

「あれは大学の広報用キャラクタだよ.論文はね,作者自身がキャラクタなんだよ.論文を書くとき,作者は《研究者》というキャラクタを演じながら書くんだ.論文という作品形態は,作者というキャラクタを極端に背景化させて,描いた物語なんだ.情報量がとても多いから,可読性はどうしても低くなる.ようするに情報伝達手法としてはもっとも効率的だけど,難しいから読みにくいんだ.だから情報の濃度を薄めて,そのぶん物語とキャラクタにウェイトをかければ,多くの人にとって読みやすい作品ができあがる」

「なるほどー.そんなふうに考えたことなかったな.あ,これもひょっとして《モデル》ってやつ?」

「数学を使わないモデルだよ.文学という現象の構造を抽象化したモデルだ」

青葉は花京院の説明を聞いて,彼がどうして理学部から文学部にやってきたのか,少し理解できたような気がした.きっと彼にしてみれば研究と文学は連続的につながっているのだろう.

(終わり)

 参考文献


おかげさまで 『その問題,数理モデルが解決します』をベレ出版さんより刊行いたしました.

「数学は人生に不要だ」と考える青葉の抱える問題を,数学好きの花京院が数理モデルで解決する,という内容の入門書です.

f:id:hamada7418:20181221154801j:plain:w200f:id:hamada7418:20181221154634j:plain:w200

https://www.amazon.co.jp/dp/4860645685/

サンプルとして数式の出てこないショートストーリーを書いてみました.

どんな風に2人が対話しているのか,その雰囲気を感じていただければ幸いです.

以上,『その問題,数理モデルが解決します』刊行記念SSでした.


DAアルゴリズムにかんするメモ

RによるDAアルゴリズムの再現

DAアルゴリズム(GSアルゴリズムとも呼ばれます)のコードをRで書きなおして簡単な実験をやってみました. オリジナルのコードは奥村晴彦先生の『C言語による最新アルゴリズム事典』技術評論社を参照しました.

Cではprefix increment,つまりループの中に ++ v を入れることでコードが簡略化されています. 同じ記法はMathematicaでも可能だったので,Mathematicaへの翻訳は簡単でした. しかしRには仕様上prefix incrementがないらしく,その部分で多少の変更を余儀なくされました.

DAアルゴリズムそのものの説明は,ここでは省略します (近々刊行予定の本で説明していますので、よろしければそちらを参照してください(宣伝)).

『その問題、数理モデルが解決します』 https://www.amazon.co.jp/dp/4860645685

###################################################
# GS Algorithm 概略
###################################################
#関数gsは引数として,男性の選好(proposer)と女性の選好(receiver)を
#入力すると,GS Algorithmに基づく男女のマッチングを出力します

#二つの引数proposerとreceiverは男女である必要はありませんが,
#先行研究にならって,以下男女の喩えを使います.

# はじめにproposer, receiverの選好を,集団別にlist型で入力します.
# 例えば男性が求婚する時,proposerとして
# list(c(2, 3, 1), c(1, 2, 3), c(3, 1, 2)) を入力する
# これは
#男性1の選好は  女2>女3>女1
#男性2の選好は  女1>女2>女3
#男性3の選好は  女3>女1>女2
#を意味します.
#同様にreceiverは,女性が持つ男性に対する選好のリストです.

関数内のローカル変数の意味は以下の通りです.カウンタがたくさんあるため,構造は結構複雑です.

########################################
# gs内のローカル変数の意味
########################################
#prefb:男性が持つ,女性に対する選好. prefg:女性が持つ,男性に対する選好.
#nextgirl:プロポーズする順番を記録.各男性が独立に持つ.振られる度に1増加
#tempboy:仮マッチした男性のリスト,最初はn+1(番兵)が入る.
#gid:男性がプロポーズする相手女性のid, 
#b:for用のカウンタ.男性1~nに対応
#s:現在プロポーズ中の男性id
#t:振られた男id.この値をsに代入する事で,
#「今回振られた男t」が自動的に次のsになる.

結果以下のようなコードになりました *1

##################################
# 関数gsの定義
##################################
gs130<- function(proposer,receiver)
{ n<-length(proposer)
prefb = proposer #男性の好みは引数proposerをそのまま入力
z=receiver  #一旦,女性の選好リストを別のzに代入する
for(j in 1:n){
    for(i in 1:n){  
        z[[j]][i]=which(receiver[[j]]==i)}#for j
}#for i
prefg =z #女性の選好を順位表記にする
print("男性の選好") 
print(paste0(matrix(prefb))) 
print(paste0("女性の選好") )
print(paste0(matrix(receiver)) )
for(i in 1:n){
    prefg[[i]] = append(prefg[[i]], n + 1)
}#女性の選好に番兵n+1を追加
nextgirl = rep(1, n)  #個人別のプロポーズ用カウンタ.1位からスタートして,2位,3位と順番に進む
tempboy = rep(n + 1, n)  # while条件用 番人sentinel  
for(b in 1:n){
    s = b
    while(s != n + 1 ){ #sが1,2,3である間ループ*)
        k <- nextgirl[s]
        gid =prefb[[s]][k]#k=1からスタート
        if(prefg[[gid]][s] < prefg[[gid]][tempboy[gid]])#現在の相手順位<仮の相手順位
        {   t = tempboy[gid]#仮の相手を捨てる
        print(c("振られた男=",t))#初回はt=n+1,入れ替わるとt=以前に仮マッチした男性ID 
        tempboy[gid] = s # 仮の相手を入れ替え
        s = t  # 次のs(proposer)を指定する.
        } #if 条件が真の場合の実行文終わり
        else{#もしsが拒否されたらnextgirl[s]に1を足す
        nextgirl[s] <- nextgirl[s]+1
        }
    } # whileここで終わり
}# for ここで終わり

for(i in 1:n)#結果の表示
{print(paste0("男性",tempboy[i], " と 女性",i))}
}#gs ここで終わり

さて,ランダムな選好を関数

lg <- function(n){replicate(n,list(sample(1:n)))}

で発生させると,関数gsは次のように実行できます.

gs130(lg(10),lg(10))

これで10対10の,マッチング結果が出てきます.n=3の場合の計算結果を確認してみましょう.

gs130(lg(3),lg(3))
[1] "男性の選好"
[1] "c(3, 1, 2)" "c(3, 1, 2)" "c(1, 3, 2)"
[1] "女性の選好"
[1] "c(1, 3, 2)" "c(1, 3, 2)" "c(3, 1, 2)"
[1] "振られた男" "4"         
[1] "振られた男" "4"         
[1] "振られた男" "2"         
[1] "振られた男" "4"         
[1] "男性3 と 女性1"
[1] "男性2 と 女性2"
[1] "男性1 と 女性3"

"振られた男" "4" は番兵(センチネル)です.whileループを抜け出す条件として使います.

このアルゴリズムは相手集団に対する完全な選好が必要なのですが,現実的には全員に順位をつけるのは難しいでしょう. 一方でアルゴリズムの原理から,選好によっては10人中5人くらいまで順位をつけておけば事足りる場合もあることは想像できます.

そこで10対10マッチングを実行した場合に,選好リストの下位何番までの情報が必要になるのかを調べてみました.

f:id:hamada7418:20181107180016p:plain

この図は,最も望ましくない相手とペアになったマッチングがどのくらいの確率で生じるのかを表しています. 500回ランダムな選好を発生させてGSアルゴリズムでマッチングを作り,それを100回計測した分布です. 例えば一番左は,proposerから見て一番順位の低い人とのペアが実現したマッチングの割合です.全体の3%程度このようなマッチングが存在します. 逆に言えば,97%は最後まで選好を表明しなくても安定マッチングができるという意味です.

図の横軸-2は1~2番目に低い順位とのマッチング,-3は1~3番目に順位の低い人とのマッチングが生じた割合を示しています.

DAアルゴリズムはとても単純なアルゴリズムなのに,結果は強力です(提案側最適な安定マッチングを実現する).

コードを書くよい練習になるので,みなさんも一度再現に挑戦してみてください.

*1:ちょっと細かい話ですが,関数の内部でreceiverの選好をidではなく順位に変換しています. 例えば,ある女性の選好がc(2,3,1) 男2>男3>男1 のとき,この選好を (男1は)3番目, (男2は)1番目, (男3は)2番目  という順位のベクトルに変換します. つまり c(2, 3, 1)] というreceiverの選好の入力をいったん  c(3, 1, 2) と変換します. これは receiver側で相手を留保するかどうかを順位を参照して決定するためです

Mathematica+TikZ

タイトル通り,MathematicaでつくったグラフをTikZでお化粧してTeXに取り込む方法のメモです.

まず素材となるグラフをMathematicaでちゃっとつくります.

Plot[PDF[LogNormalDistribution[0, 1], x], {x, 0, 5}, 
 PlotStyle -> {Black}, Frame -> True,
 FrameTicks -> {{{0, 0.2, 0.4, 0.6}, None}, {{0, 1, 2, 3, 4, 5}, 
    None}}]

このグラフにlogn_pdfという名前をつけて保存しておきます.

次にTeXのFigure環境の中にTiKZでグリッドを描きます.\includegraphicsで先ほどつくったlogn_pdfを 読み込むのですが,tikzpicture環境の中で読み込むのがポイントです.

\begin{figure}[H]
    \centering
\begin{tikzpicture}
\draw[step=1,black!0] (0,0) grid (12,5);
\node[anchor=south west] (image) at (1.5,-0.5)
    {\includegraphics[width=0.7\linewidth]{logn_pdf}};
\node [right] at (4.432,3.5) 
{$\displaystyle \frac{1}{\sqrt{2 \pi \sigma^2} x}\exp \left\lbrace -\frac{(\log x-\mu)^2}{2\sigma^2}\right\rbrace$};
\node [right] at (0.7,4.8){$f(x)$};
\node [right] at (10,0){$x$}; 
\end{tikzpicture}
{\small 対数正規分布の確率密度関数$f(x)$のグラフ. $\mu=0,\sigma=1$}
\end{figure}

できあがりー.

f:id:hamada7418:20181018230630p:plain

グラフの中と周りにTeXの数式がきちんと挿入されました. R MrakdownJupyter notebookでも同様のことができます. ぜひお試しを.

TeXの柱とノンブルの設定

jsbookクラスをつかった場合の柱の設定方法を忘れないためにメモを残しておきます.

最終的な出力目標は

偶数頁のヘッダ左側に章タイトル,奇数頁のヘッダ右側にセクションタイトル

偶数頁のフッタ左側に頁番号,奇数頁のフッタ右側に頁番号

であると仮定します.これをデザインAと呼ぶことにしましょう.

\chapterコマンドはデフォルト設定で開始1頁目に\thispagestyle{plain}を実行するので,章のはじめだけ設定したスタイルが表示されないことが しばしばあります.まずはこれを制御します.

直感的に分かりやすいやり方は,\chapterコマンドのすぐ後に自分で\thispagestyle{empty}\thispagestyle{fancy} を使って,\thispagestyle{plain}を上書きしてやることです.

しかしこの方法だと全ての\chapterに同じ命令を書き込む必要があるのであまりスマートではありません. \newcommandでコマンドを定義すれば多少はラクになりますが,章が多いとはやはり煩雑です.

一番簡単な方法は pLaTeX2e 新ドキュメントクラス の指示通りにプリアンブルに

\makeatletter
\renewcommand{\chapter}{%
  \if@openright\cleardoublepage\else\clearpage\fi
  \global\@topnum\z@
  \secdef\@chapter\@schapter}
\makeatother

と書く方法です.この\renewcommandによって,章開始1頁目は自分で設定した柱のスタイルに一致します. どうしてこれでうまくいくのか,実はよくわかっていません.おそらく, jsbook.cls\chapterの定義にある\plainifnotemptyの削除により, \thispagestyle{plain}を実行しないのではないかと推測します.\plainifnotemptyというコマンドは全体のpagestyleemptyでないならplainにせよ,という命令です. ですからこの命令が有効だと,empty以外のpagestyleを使うと,無視されるという現象が生じるのでしょう.

次に, デザインAを実現するには\thispagestyle{fancy}の中身を設定する必要があるので, これを調整します.例としてはこんな感じです.

\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhead{} % clear all header fields
\renewcommand{\chaptermark}[1]{\markboth{第\ \thechapter\ 章\,\, #1}{}}
\fancyhead[RO]{{\footnotesize \rightmark}}
\fancyhead[LE]{{\footnotesize \leftmark}}
\renewcommand{\headrulewidth}{0pt} %ヘッダの罫線を消す 
\fancyfoot{} % clear all footer fields
\fancyfoot[LE,RO]{{\footnotesize \thepage}}%頁番号の表示を偶数頁と奇数頁に分けて指定.
\cfoot{}

\fancyheadfancyfootのオプションは

R: 右側 L: 左側 O: 奇数頁 E: 偶数頁

を組み合わせて指定します.例えば \fancyfoot[LE,RO]{{\footnotesize \thepage}} は,「LEが偶数頁のフッターの左側」「ROが奇数頁のフッターの右側」に頁数を示す,という指定です.

もし\frontmatter\chapterを使っている場合は,

\renewcommand{\chaptermark}[1]{\markboth{第\ \thechapter\ 章\,\, #1}{}}

の部分を

\renewcommand{\chaptermark}[1]{#1}{}}

などに改めて定義してください.こうすることで柱を「第0章 イントロダクション」ではなく「イントロダクション」という見出しに変更できます. 一方で\mainmatterは「第1章 先行研究のレビュー」のように,{\markboth{第\ \thechapter\ 章\,\, #1}を使います.

要するにポイントは

  1. \renewcommand{\chapter}で章開始頁を修正

  2. \frontmatter用に\pagestyleを定義

  3. \mainmatter用に\pagestyleを定義

です.\frontmatter\mainmatterの柱が同じなら定義は1種類でOKです.

tikzのメモ

描き方を忘れないためのメモです.

tikzで関数のグラフを描きたいとき,

\draw [very thick] (始点座標) to [out=5, in=-90] (終点座標);

のようにオプション[out=5, in=-90]で始まりの角度と終わりの角度を指定すると簡単に描けます. 効用関数系のグラフだったら,大体これでいけるでしょう.

f:id:hamada7418:20180830200233p:plain

正規分布のような曲線を描きたい場合は

\draw plot[smooth] (\x, {exp(-0.5* \x * \x)/sqrt(2 * pi)});

こんな感じで関数を記述すればtikzの描画機能だけでグラフを作成してくれます. なおオプションの[smooth]を忘れるとカクカクした線になってしまいます.

f:id:hamada7418:20180830200253p:plain

座標を書く場合は\foreachコマンドで繰り返しを記述すると楽できます.

\draw[thick, ->] (-5,0) -- (5,0) node [below] {$x$};%x軸
\draw[thick, ->] (0,-1) -- (0,5) node [left] {$y$};%y軸
\node at (0,0) [anchor=north east] {O};
\foreach \x in {-4,-3,-2,-1,1,2,3,4}
\draw (\x cm,.2) -- (\x cm,0) node[anchor=north] {$\x$};

最近はMathematicaで図を作成→pdf出力→include→tikzでお化粧 という方法を多用していました.tikzだけで完結するなら,そのほうが楽かもしれません.