The One with ...

授業・研究用のメモなどを公開しています

微分方程式と状態空間モデル

イントロ

微分方程式をベースにした統計モデルに関して興味深い論文を読んだので紹介します.

Naik, Prasad A., Murali K. Mantrala and Alan G. Sawyer, 1998, Planning Media Schedules in the Presence of Dynamic Advertising Quality, Marketing Science, Vol. 17, No. 3, pp. 214-235

https://www.jstor.org/stable/193228

広告のタイミングと売上の関係を時系列データで分析するという内容です. この中で特に私が注目したのは以下の部分です.

In Table 2, we briefly describe the extant models we have chosen to compare with the proposed model, namely, the models of Vidale and Wolfe [ 1957 ] , Nerlove and Arrow [1962], Brandaid [ Little 1975 ], Tracker [ Blattberg and Golanty 1978 ], and Litmus [ Blackburn and Clancy 1982 ]. It is interesting to note that these five well-known models can be cast into the same state-space form that was employed to estimate the proposed model.

筆者らは,微分方程式をベースにした理論モデルが状態空間モデルで表現できる(埋め込むことができる)と述べています.

はじめはよく分からなかったのですが,彼らのアイデアを理解するにつれ,次第に「なるほど確かにこういう方法は有用かもしれない」と思い至りました.

Naik らのアイデア

論文内で参照している five well-known models はいずれも微分方程式モデルです.Naik らのアイデアの骨子は,微分方程式を陽表的に解くのではなく,離散化した微分方程式を状態空間モデルの状態方程式として使う,というものです.

たとえばVidal-Wolfeモデルは $$ \frac{ds}{dt}+\left( \frac{rA}{M} + \beta \right)s= rA $$ という微分方程式なので,これを陽表的に解くと,特殊解は $$ s=f(t) $$ という形で表すことができます(厳密には条件分岐しますがちょっと省略して書いてます).したがってこの特殊解を決定論的パラメータとしてもつ観測誤差関数を定義すれば,カルマンフィルタでパラメータを推定したり,ベイズモデル化して予測分布を計算できます. 一方でNaikらの方法は離散化によって微分方程式を近似するため,わざわざ特殊解を特定する必要がありません.それゆえ計算がむずかしい微分方程式モデルに対して適用可能です. integrate_ode( )を使って微分方程式stanに計算させることも可能ですが,簡単な条件下では「微分方程式+状態空間モデル」のほうが手軽です. *1

パラメータリカバリ

Naikらのアイデアをパラメータリカバリで確認してみましょう. まず分析用のデータは次の微分方程式 $$ \frac{dy}{dt} = r y \frac{m - y}{m},  t=0 のとき y=y_0 $$ から得た特殊解 $$ y= \frac{\exp\{ rt \}m y_0 } { m- y_0+ \exp\{ rt \} y_0} $$ に誤差を足した値であると定義します.

Rコードで書くとこんな感じです.

r <- 0.1  #増加率パラメータ 推定対象
m <- 100  # 上限
y0 <- 5  #初期値
n <- 100  #データ数.
sigma <- 1
for (t in 1:n) {
Y[t] <- rnorm(1, (exp(r*(t-1))*m*y0) / (m - y0 + exp(r*(t-1))*y0),sigma)
}

このデータを,つぎの三つの統計モデルで予測します.

  1. 決定論的関数+観測誤差モデル
  2. ローカルレベルの状態空間モデル
  3. 離散化した微分方程式状態方程式とする状態空間モデル

それぞれのstanコードは以下の通りです.*2

まず決定論的関数+観測誤差モデル

data {
  int n;//期間
  real Y[n];//観測値
  real y0data;//初期値
real m;
}

parameters {
  real<lower=0,upper=2> s_v;  // 観測誤差の標準偏差
  real<lower=0,upper=1>r;//増加パラメータ
}

model {
target += uniform_lpdf(s_v|0,2);
target += uniform_lpdf(r|0,1);
target += uniform_lpdf(mu|0,100);
for(i in 1:n) {
target += normal_lpdf(Y[i]|(exp(r*(i-1))*m*y0data) / (m - y0data + exp(r*(i-1))*y0data), s_v);
  }
}

generated quantities{
    real y_pred[n];
  for(i in 1:n){
    y_pred[i] = normal_rng((exp(r*(i-1))*m*y0data) / (m - y0data + exp(r*(i-1))*y0data), s_v);}
}

次にローカルレベルの状態空間モデル

data {
  int n;//期間
  real Y[n];//観測値
}

parameters {
  real<lower=0,upper=2> s_w;  // 過程誤差の標準偏差
  real<lower=0,upper=2> s_v;  // 観測誤差の標準偏差
  real<lower=0,upper=100> mu[n];// 状態の推定値(水準成分)
}

model {
target += uniform_lpdf(s_w|0,2);
target += uniform_lpdf(s_v|0,2);
for(i in 2:n) {
target += normal_lpdf(mu[i]|mu[i-1], s_w);//いわゆるローカルレベルモデル
  }
for(i in 1:n) {
target += normal_lpdf(Y[i]|mu[i], s_v);
  }
}

generated quantities{
    real y_pred[n];
  for(i in 1:n){
    y_pred[i] = normal_rng(mu[i], s_v);}
}

最後に微分方程式状態方程式とする状態空間モデル

data {
  int n;//期間
  real Y[n];//観測値
  real m;//上限
}

parameters {
  real<lower=0,upper=2> s_w;  // 過程誤差の標準偏差
  real<lower=0,upper=2> s_v;  // 観測誤差の標準偏差
  real mu[n];// 状態の推定値(水準成分)
  real<lower=0,upper=1>r;//増加パラメータ
}

model {
target += uniform_lpdf(s_w|0,2);
target += uniform_lpdf(s_v|0,2);
target += uniform_lpdf(r|0,1);
for(i in 2:n) {
target += normal_lpdf(mu[i]|mu[i-1]+r*mu[i-1]*((m-mu[i-1])/m), s_w);
//状態方程式を微分方程式(ロジスティック)で代替する
  }
for(i in 1:n) {
target += normal_lpdf(Y[i]|mu[i], s_v);
  }
}

generated quantities{
    real y_pred[n];
  for(i in 1:n){
    y_pred[i] = normal_rng(mu[i], s_v);}
}

結果を比較します.

微分方程式の特殊解+誤差モデル

f:id:hamada7418:20210325125902p:plain
微分方程式特殊解+誤差

この統計モデルはデータを生成した真の分布と一致しているので,まあ当然の結果と言えるでしょう.

ローカルレベル状態空間モデル

f:id:hamada7418:20210324152538p:plain
ローカルレベル状態空間モデル

これは1期前の状態+システム誤差で現在の状態が決まるモデルです.

離散化した微分方程式状態方程式とする状態空間モデル

f:id:hamada7418:20210324152618p:plain
離散化した微分方程式状態方程式とする状態空間モデル

今回実験してみたかったモデルです.いい感じにフィットしているようです.

各統計モデルの自由エネルギーを計算したデータへのフィットを確認します.

  1. 微分方程式の特殊解+誤差モデル 143.1769
  2. ローカルレベル状態空間モデル 228.1158
  3. 微分方程式状態方程式とする状態空間モデル 148.0088

特殊解を使った正解のモデルと,微分方程式+状態空間モデルで近似したモデルの自由エネルギーはほぼ等しい値でした.

時間以外の説明変数があるデータのほうが,Naikらのアイデアは真価を発揮するのではないかと予想します.

*1:ただし連立微分方程式は integrate_ode( ) でないと無理かもしれません

*2:状態空間モデルのtarget記法はもしかしたら間違っているかもしれません.ここが間違っていれば自由エネルギーの計算結果が少し変わります.予測分布の計算には支障がないはずです.

ケースファンの回転数調整

PCがうるさい

なんだかPCがうるさい.HWMonitorで確認したらケースファンの回転数が900RPMでした.

マザーボードの設定を変えたら静かになったのですが,多分3日後にはやり方を忘れてしまうので記録を残しておきます.

4ピンではなく3ピン接続だった

ケースを開けて,ケースファンのコネクタを外してみると4ピンではなく3ピンでした. 3ピンの場合は確か電圧で回転数を変更できるはずなので,BIOSもといマザボUEFIを起動します. これはMSI社製のマザボの設定画面です.起動時にdelキーを押すとこの画面が出てきます(メーカーによってこのキーがF2だったりします).

Fan InfoからPUMPを選択.ここはマザボによって異なるはずなので,ケースファンに対応しているものを選択してください.私はCPUクーラーの事情でケースファンをPUMP用コネクタに接続しているので,それを選びます(3ピンの受けに接続してたら,そもそもこの騒音問題は起こらなかったと推測).

f:id:hamada7418:20210320123936j:plain

マザボの設定

ファン制御の設定をよく見ると,PUMP用のコネクタが4ピンなので,PWMになっています.PWMで制御できない3ピンコネクタを挿しているせいでしょうか,回転数が1000RPMを超えています.どうやらここが騒音の元凶のようです.

f:id:hamada7418:20210320121637j:plain

PWMからDCに変更します.340RPMに下がり,かなり静かになりました.主観的幸福感が0.5ポイントほど上昇しました.

f:id:hamada7418:20210320122042j:plain

メーカーによって多少違いはありますが,設定は似ているはずなので,ファンがうるさいなと思ったら試してみてください(ただし回転数を下げると温度は上がるので,その辺は適当にナニをアレしてください).

CPUファン

このマザボに使用しているCPUファンは横方向にファンを2個つけて風を前から後ろに流すタイプなのですが,ケース後方側のファンの取り付け方向が間違っていました.これをついでに直したおかげで,静かになったうえに温度も下がりました(オイ).

オブジェクトやファイル名をループ内で割り当てる(R)

ファイル名を直書きするのがめんどうくさい

50個のデータセットを使い複数のモデルで推定して自由エネルギーを計算する必要があったのでコードを書いてみました.

多分忘れるだろうから,ここに貼り付けておきます.

fevalues.m0 <- c()
for (i in 1:50) {
  stanfit <-stan(file='model0.stan',iter=10500,warmup=500,data=f2(i),chains=4)
  saveRDS(stanfit,c(paste('fit/fe', i,'_0.rds', sep = '')))
  fevalues.m0 <- append(fevalues.m0,-logml(bridge_sampler(stanfit,method="warp3"))
  )
}  

ここで関数f2( )はstan用に各データセットをlist型に整形する関数です. stan( )の引数dataに関数の出力を渡しても問題ないようです.

途中でRがコケてもstanfitが残るようにsaveRDSで逐次保存します.このときpaste( )forループで回してデータセットに対応したファイル名を自動作成します(10個くらいなら直打ちのほうが早いでしょう.50×3個のファイル名はさすがに書く気がしない). saveRDS内でpaste( )を使っても問題ないようです.

stanfitは上書きされて1つしか残らないため,メモリ不足でRがクラッシュすることを防ぎます(多分)

計算した統計量(ここではブリッジサンプリングから計算した自由エネルギー)は空のベクトルfevalues.m0appendしていきます.

モデルを変える場合は,stan( )のモデル部分だけ変え,上記のfor文をコピペして使い回します. 1回のループにまとめることもできますが,かえって可読性が悪くなると判断しました. このくらいで分割した方が,数ヶ月後の自分には理解しやすいでしょう.

呼び出したstanfitからの計算

readRDSで保存したstanfitを使ってbridge_sampler( )を適用すると,「modelがないぞ」とrstanに怒られます.

コンパイルしたモデルは保存する設定にしとるやろがい!」と悪態をつきたくなる気持ちをぐっとこらえ,stanfitに対応したモデルをchains=0で再コンパイルします(これは時間はかかりません).

mod<-stan(file='model0.stan',iter=10500,warmup=500,data=f2(i),chains=0)
-logml(bridge_sampler(stanfit,mod,method="warp3"))

これなら計算してくれます.

他に解決策があるのかもしれませんが,今回はこれでOKとします.

数理モデル本の続編を書きました

はじめに

『その問題,やっぱり数理モデルが解決します』をベレ出版さんから刊行しました.以下,その内容をざっと紹介します.

f:id:hamada7418:20200918225126j:plain

本書は次の4つのパートから構成されます.

  • 1~3章 数式の理解のしかた
  • 4~8章 複占モデルを中心としたゲーム理論
  • 9~12章 統計的推測と回帰モデル
  • 13~14章 広告の微分方程式モデル

はじめから順番通りに読んでもいいし,興味のあるパートから読んでも構いません.

前作『その問題,数理モデルが解決します』と同様に,登場人物達の対話を通して数理モデルを紹介します.世界観は共通ですが,内容は独立しているので,今作から読み始めていただいても問題ありません.数学の難易度は前作よりも少し易しくなるよう心がけました.

対象読者

本書の読者として以下のような方を想定しています.

  • 数学って何の役に立つの?と考えている中高生
  • 数理モデル*1に興味のある大学生
  • 統計学を使うけど,仕組みがよく分からないという研究者
  • データ分析について勉強してみようかなと考えている社会人

特徴

「数理」モデルの理解はなかなか困難です.そこでなるべく分かりやすく,かつ,数式はちゃんと使いつつ数理モデルの使い方や作り方を紹介するために,この本を書きました.おそらく,キャッチーな表紙からは想像できないほど本書にはたくさんの数式が登場します.しかしこれにはちゃんとした理由があります. 1つの計算過程を省略せずに書くためには,必然的にたくさんの式が必要です.逆に言うと,途中の変形ステップを省略せずに式がたくさん書いてある本ほど,実は簡単に読めるのです.

多くの初学者は,数式間で省略されたステップが理解できずに挫折します.ですので本書は,(前作同様に)可能な限り式変形を丁寧に書くよう心がけました*2

各章のタイトル・難易度・トピック

タイトル 難易度 トピック
第1 章 しっくりこない数式の読み方とは? 総和記号、添え字の読み方、平均値
第2 章 確率密度関数ってなに? 定義の読み方、グラフの描き方、パラメータの役割
第3 章 確率と積分の関係とは? 正規分布、測定誤差、確率密度関数のつくり方
第4 章 モデルってなに? ☆☆ 競争市場、需要量、代替財、均衡価格、独占モデル
第5 章 競争で損をしない戦略とは? ☆☆ 複占モデル、ゲーム理論ナッシュ均衡
第6 章 自分でモデルをつくる方法 モデルのつくり方、モデルの一般化、モデルで自由に遊ぶ
第7 章 先手が有利な条件とは? ☆☆ 後ろ向き帰納法、展開形ゲーム、情報集合
第8 章 競争に負けない価格設定とは? ☆☆ 代替性、モデルのバリエーション、検証
第9 章 売り上げを予測するには? ☆☆ 回帰分析、最小2乗法、微分偏微分
第10 章 確率モデルでデータを分析するには? ☆☆☆ 確率変数、OLS推定量正規分布の再生性
第11 章 仮説検定ってどうやるの? ☆☆☆ 定量の確率分布、t分布、解釈の注意点
第12 章 観測できない要因の影響を予想するには? ☆☆ 欠落変数バイアス、期待値、不偏推定量、固定効果
第13 章 アプリの利用者数を予測するには? ☆☆ 微分方程式、置換積分法、変数分離形
第14 章 広告で販売数を増やすには? ☆☆☆ 1階線形微分方程式積分因子

メッセージ

私は大学院に入るまでは,「数学は自分の人生に無用だ」と考えていました. しかし今では数理モデルを勉強したり数理モデルを使うことは,人生に欠かせない楽しみの一つです. 本書をきっかけに,数理モデル(数学)っておもしろいかも,と思っていただけたら望外の喜びです.

以下,各章の概要です.ちょっと長いので気になる人だけどうぞ.

各章の内容

第1 章 しっくりこない数式の読み方とは?

難易度☆ 総和記号、添え字の読み方、平均値

ときどき,$\sum$が苦手だ,という学生さんの声をよく聞きます. $$ \sum_{i=1}^{n}x_i=x_1+x_2 + \cdots +x_n $$ は分かるんだけど, $$ \sum_{i=1}^{n}a=an $$ は,なんかしっくりこない,と.

第1章ではこのような総和記号の意味からはじめて,《分かったような分からないようなしっくりこない数式》に出会ったときに,どうやってその式を理解するのか,その一般的な方法を解説します.

第2 章 確率密度関数ってなに?

難易度☆ 定義の読み方、グラフの描き方、パラメータの役割

統計学や確率論のテキストを開いたところ,

次の2つの条件を満たす関数$f: R \to R$を確率密度関数という.

  1. 任意の$x \in R$について$f(x)\ge 0$
  2. $\int_{-\infty}^{\infty} f(x)dx=1$

といった定義がいきなり出てきて,「ちょっとなに言ってるかわからない」と思ったことはありませんか?

2章では確率密度関数を例に,抽象的な定義を理解する方法を解説します.

統計や数学を学ぶ上で《自分で具体例をつくる》ことはとても大切です. ここでは式を単純化したり,コンピュータを使ってグラフを描いたり,条件を変えて比較したりさまざまな試行錯誤 を通して抽象的概念を理解する方法を学びます.

第3 章 確率と積分の関係とは?

難易度☆ 正規分布、測定誤差、確率密度関数のつくり方

学生の頃,確率密度関数を計算してたら「4.023」や「5.251」のように1よりも大きな値が出てくるので,「積分したら1になる関数なのに,なんでこんな値が出てくんねん」と頭をひねったことがありました.よく考えれば別に不思議なことではないのですが,一見奇妙に思える理由を説明しつつ,積分の直感的な意味や確率密度関数を解説します.

積分と確率の一致が納得できると,確率論の世界が目の前にぱっと広がります.

第4 章 モデルってなに?

難易度☆☆ 競争市場、需要量、代替財、均衡価格、独占モデル

モデルとはなにか,という基本的な話を競争市場を例に説明しました.途中までは数式なしで,最後にほんの少し数式(独占モデル)が出てきます.数式は現実的な『意味』が与えられると,世界を表現する豊かな言葉に変化します.昔,師匠から「数理モデルの勉強を始めてから,無機質な数式が急に色めきたって見えた.まるで白黒の世界からカラーの世界に変わったような気がした」という話を聞かせてもらったことがあります*3.そんな瞬間を読者にも感じてほしいなと思ってこの章を書きました.

第5 章 競争で損をしない戦略とは?

難易度☆☆ 複占モデル、ゲーム理論ナッシュ均衡

『経済学のためのゲーム理論入門』(ギボンズ)でクールノーの複占モデルを初めて知ったとき,「このモデル,すごい!完璧や!」と感動しました.一般に,モデルのよさは $$ モデルのよさ=\frac{インプリケーションの豊かさ}{モデルの複雑さ} $$ で評価できます.複雑すぎるのもよくないし,かといって単純すぎてインプリケーションが自明なモデルもよくありません.単純さとインプリケーションの豊かさを両立させることは困難ですが,クールノー複占モデルは見事なバランスで,両者を実現しています*4.ここではそんなゲーム理論の名作,クールノー複占を例にシンプルなモデルの作りかたを解説します.

第6 章 自分でモデルをつくる方法

難易度☆ モデルのつくり方、モデルの一般化、モデルで自由に遊ぶ

モデルを拡張したり一般化する方法を実例をもとに解説しています.このあたりの話は既存のテキストだと「できて当たり前」の操作として扱われることが多いのですが,自分でモデルをつくるときには,この計算のプロセスを詳しく示す必要があるなと思ってこの章を書きました. 花京院くん曰く,モデルをつくることは「カツカレー」を作ることに似ているそうです.はたしてその意味は?

第7 章 先手が有利な条件とは?

難易度☆☆ 後ろ向き帰納法、展開形ゲーム、情報集合

モデルの拡張やヴァリエーションを複占モデルの展開形で説明します. 4章からスタートして,5,6,7,8章と独占・複占モデルを主題とする変奏曲が続きます. 新しいキャラクターも登場します.

第8 章 競争に負けない価格設定とは?

難易度☆☆ 代替性、モデルのバリエーション、検証

ひきつづきモデルの拡張やヴァリエーションを複占モデルで説明します. ここでは戦略集合が『生産量』から『価格』に変わります.やはりコアになっているのは独占モデルで,ここが最後のバリエーションです. データによるモデルの検証についても少し考えながら,次章の統計学の導入へと続きます.

第9 章 売り上げを予測するには?

難易度☆☆ 回帰分析、最小2乗法、微分偏微分

もっとも基本的なデータ分析法のひとつである回帰分析とその原理にあたるOLSを解説します. 微分がよく分からない,忘れてしまったという人はこの章で確認してください.

第10 章 確率モデルでデータを分析するには?

難易度☆☆☆ 確率変数、OLS推定量正規分布の再生性

ここでは確率変数を明示的に導入して,誤差項の分布から推定量の分布を導出します. 回帰分析を確率モデルとして構成することで,係数の検定(次章)が可能となります

第11 章 仮説検定ってどうやるの?

難易度☆☆☆ 推定量の確率分布、t分布、解釈の注意点

回帰係数の検定を解説します. OLSは分かるけれど,係数の検定はよくわからない,という読者に向けて書きました. 確率変数同士を合成して新たに確率変数をつくる,という流れがイメージできると推定量の分布がなぜ正規分布やt分布にしたがうのか分かります.

第12 章 観測できない要因の影響を予想するには?

難易度☆☆ 欠落変数バイアス、期待値、不偏推定量、固定効果

欠落変数バイアスについて解説します.単にそういうバイアスがあるんだよって話ではなく,陽表的な関数として導出する過程を示しています. 本書では触れませんでしたが因果推論や条件付き期待値回帰を勉強する場合にバイアスの理解が役に立ちます.

第13 章 アプリの利用者数を予測するには?

難易度☆☆ 微分方程式、置換積分法、変数分離形

微分方程式モデルの基礎を説明します.積分の計算方法を知らない,あるいは習ったけど忘れてしまった,という方はこの章で復習してください. 置換積分の計算は,確率変数の計算にも役立つので,ぜひ紙に式を書いて納得してください. 積分の計算は面倒ですが,この計算によって確率論の理解がぐっと深まります.

第14 章 広告で販売数を増やすには?

難易度☆☆☆ 1階線形微分方程式積分因子

最後は広告の効果を分析する微分方程式モデルの紹介です.Vidal & Wolfe(1957)の古典的な論文をベースに,新キャラクターのヒスイが発展的なモデルを提案します. 異なる広告タイプを比較するために積分を利用するあたりが面白いモデルだと思います.13、14章では積分が大活躍します.


最後まで読んでいただき,どうもありがとうございました.

*1:数理モデルは,数学を使って現象を表現・説明・予測する方法の総称であり,観察した現象から本質を見抜く洞察やデータを分析するための指針(理論)を与えてくれます.現象が成立するメカニズムを数理モデルで特定できれば,不確実な未来を予測できます.コロナウィルスの感染症対策でSIRモデルおよびその拡張的数理モデルは,しかじかの仮定の下で感染者数はこうなると明確な予想を導きました.

*2:逆にいえば,既に数学を使った議論に慣れている中級者以上の読者には,少しもの足りないかも知れません.もう少し難易度が高くてもOKかつ,ベイズ統計モデルに関心のある方は,『社会科学者のためのベイズ統計モデリング』をご覧ください.共著者の清水さんによる紹介記事はこちらです

*3:私は師匠のつくった数理モデルの論文を読んだとき,最初は意味が分からなかったのですが,計算を繰り返すうちに小さな世界で人が動いている様子を想像できました.いえ,酔っ払ってヘンなモノを見たわけではありません

*4:ゲーム論のモデルにはたくさんの名作がありますが,クールノー複占は個人的にはベスト3に入るくらい好きなモデルです.ちなみに,この章を書くためにクールノー複占をゲーム理論で最初に表現した人が誰かを調べてみました(クールノーはゲーム論が生まれるずっと前に複占モデルを定式化していました).どうやらギボンズよりもフリードマンの方が早かったようです

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でした.