Talent versus Luck modelのパラメータリカバリ(Julia)

この記事の内容

運と才能のシミュレーション(Pluchino et al 2018)を数理モデルとして一般化した論文

Hamada, Hiroshi, 2023, Luck of Outcome in the Talent Versus Luck Model, Advances in Complex Systems Vol. 26, No. 04n05.

をもとにjuliaでパラメータリカバリを再現するためのメモ．

Juliaの導入

こちらの記事を参考にJuliaをイントールする

qiita.com

さらにVSCodeのJulia用エクステンションをインストールする．

www.julia-vscode.org

VSCode経由でJuliaを書けば，コードの書き方が分からないとき　Github Copilot　に教えてもらうことができる．

コードを書く

RStanとStan用に書いたコードを流用してjulia用のコードを作成する．基本的な文法は同じなので，分からない箇所だけ Github Copilotに教えてもらうとよい．

翻訳した結果がこちら

#パッケージの読み込み．
using Turing
using Distributions
using StatsPlots
using DataFrames
using Plots


# talent,capital,Sn=log(capital)*log(c0*b)^(-1) の生成，c_0=1,b=1.2,m=人数，it=繰り返し回数
# 一般性を損なうことなくc0=1を仮定する
function TvLsim(m, it, mu, b)
    t = rand(Normal(mu, 0.1), m)
    c = ones(m)#　初期資本c0=1とおく
    Sn = ones(m)#　agentごとに総成功回数を記録する．    
    for j in 1:it
        for i in 1:m
            if rand() < t[i]
                c[i] *= b
                Sn[i]=log(c[i])*log(b)^(-1) #Snの理論値はHamada2023参照.log(b)はc0=1を仮定したのでlog(c0*b)の略
            end
        end
    end
    return DataFrame(talent = t, capital = c, Sn=Sn)
end

# 人工データを生成．以下で使うのはSnだけだが，後のためにdataframeを作っている
m,it,mu,b=1000, 100, 0.6, 1.2
dataf=TvLsim(m, it, 0.6, 1.2)
Sn=dataf[:,:Sn];

dataf

#才能の平均muの事後分布を計算するmcmc用統計モデル
@model function TvL_model_mu(Y)
    # 事前分布の設定．sは総成功回数の標準偏差．μ∈(0,1)は才能分布の平均位置
    s ~ Uniform(0,100)
    μ ~ Uniform(0,1)
    
    # サンプルサイズ
    N = length(Y)
    it=100
    for i in 1:N
        Y[i] ~Normal(it*μ, s)
    end
end

#sの理論値を確認 Hamada2023参照
it, μ, σ = 100, 0.6, 0.1
sqrt((it^2)*(σ^2)+it*(μ-μ^2))

# モデルのコンパイル
model = TvL_model_mu(Sn)

# MCMCサンプリング
chain = sample(model, NUTS(), 3000)

# サンプル結果の要約
summary(chain)
println(chain)

最初はjupyter notebookで書いてみたが，.jlファイルをVSCodeで編集したほうが圧倒的に書きやすかった．ノートブックを公開しないならこちらで十分だろう．

結果

目的となる真のパラメータ値はμ=0.6とs=11.1355である．（理論的にはs=sqrt((it^2)*(σ^2)+it*(μ-μ^2))だが，この明示的な関数の全てのパラメタは資本量のデータだけでは識別できない．今後プロキシ変数を検討する必要がある）

Stanで計算した結果と概ね同じ出力を得た．コンパイルの時間が短いので，その分ストレスがない．サンプリングの時間は大差ないかもしれないが，モデルを修正しながら試す段階ではjuliaは手軽で良いと思う．

モデルに渡す外生変数を気軽に変更できるのもよい． Stanの場合，渡すデータをR側のlistで指定してからStanコードでもそれを引数として明示しないといけないので若干面倒である．細かなことだが，そういう手間がjuliaでは省けるのが結構便利だと感じた．

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ess_per_sec 
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64       Float64

           s   11.0044    0.2467    0.0042   3427.4607   2005.7838    0.9998      692.1366
           μ    0.6036    0.0036    0.0001   3369.6896   2060.8297    1.0000      680.4704