2024-04-25

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

（はてなblogにVSCodeから直接コピペはできない．ショートカットキー（win+Shift+s）でスクリーンショットをとればコピペできる．画像として保存する手間が省ける）

おまけ

\バックキーの後にTeXコマンドを入力すると，juliaソースにUnicode入力できる． TeXに慣れているとギリシア文字などの入力が手軽で，論文と計算との接続がよりスムーズになるだろう．

2023-10-05

VSCode上でのTeX設定メモ

三日後に忘れるのでメモを残しておく

前準備

VSCodeインストール
LaTeX-workshop extensionインストール
TeXliveインストール
Github copilotインストール

設定

こちらのサイトを参考にして.latexmkファイルを作成し，ルートフォルダに置く． Winの場合，C:\Users\USERNAME直下でよいらしい．

www2.yukawa.kyoto-u.ac.jp

VSCodeでLaTeX-workshopのsettings.jsonを修正する．やはり上のサイトを参考にして，修正後は以下のように設定する．

        {
            "name": "latexmk",
            "command": "latexmk",
            "args": [
                "-synctex=1",
                "-interaction=nonstopmode",
                "-file-line-error",
                "-halt-on-error",
                "-pdfdvi",
                "-outdir=%OUTDIR%",
                "%DOC%"
            ],
            "env": {}
        },

変更点はエラー時の設定と，pdfdviを使うところ．

上記の設定で実行したところ，よく分からないが，TeXStudioからコンパイルするよりも格段に早くなった． TeXStudioでのビルド設定がよくないのかもしれない．（余計なDVIチェーンとか除けば早くなる？）

中間ファイルを消す

「Build LaTeX project」を押すとビルドされる。「Build -」の中にある「Clean up auxiliary files」を押すと中間ファイルを削除することができる．

ショートカットキーの登録

TeXStudioで使っているSCを使いたいので，VSCodeで登録する

キーバインド設定のjsonファイルを次のように編集する．

// Place your key bindings in this file to override the defaults
[
    {
        "key": "ctrl+4",
        "command": "type",
        "args": { "text": "$$"},
        "when": "editorTextFocus"
    },{
        "key": "ctrl+shift+m",
        "command": "type",
        "args": { "text": "\\[  \\]"},
        "when": "editorTextFocus" 
      },{
          "key": "ctrl+e",
          "command": "type",
          "args": { "text": "\\begin{align*}\n& \\\\\n& \n\\end{align*}"},
          "when": "editorTextFocus"
    }
]

ちょっと面倒だけど，一回設定すれば後はjsonファイルを使い回しできるから，よしとしよう．

2023-09-18

oTree: ラウンド毎のパラメータのランダム化

ラウンド毎にパラメータをランダム化する方法

最初はよく分からないので

def set_parameters(group: Group):#グループオブジェクトを引数にした関数，引数groupのタイプ（クラス）がGroupであることを指定
    rvec = C.Rvec.copy()
    random.shuffle(rvec)
    group.R = rvec[group.round_number-1]
    group.PROFIT = group.R*100-100
    nvec2=C.nvec.copy()
    random.shuffle(nvec2)
    group.n = nvec2[group.round_number-1]

のようにグループオブジェクトの関数としてWaitPageで実行させていた．

これでも一応シャッフルするが，復元抽出になってしまい，被験者が同じパラメータに暴露されることになり効率が悪い．

解決法

groupやsubsessionで定義しても怒られるので，結局，

#class の前に定義
v1=[3,5,6,4]
v2=[1,2,4,9]
random.shuffle(v1)
random.shuffle(v2)


class C(BaseConstants):
    NAME_IN_URL = 'single_rd'
    PLAYERS_PER_GROUP = None # 被験者は一人
    NUM_ROUNDS = 4
#    NUM_SUCCESS = random.randint(1,5)
    Rvec = v1#セッション毎に変わるグループ共通の利益率
    nvec = v2#random.shuffle(seed2) # N=10とおく

で解決した．被験者毎にシャッフルしたパラメータをラウンド数で呼び出して重複なく使う． oTreeは独自のルールが多すぎて，体系的に全然理解できない．

まあ動けばいいか．

2023-01-29

補機バッテリーの残量確認方法

冬場になると車の補機バッテリーの警告灯が毎朝点滅する．

自粛要請時代に一度バッテリーが上がってしまったので，残量を確かめるために市販の車用バッテリーチェッカーを購入した．

以下その方法のメモ

前準備

うちの車は補機バッテリーがトランク側に設置してあるため，エンジンルームを開けてプラグを接続することが出来ない．

トランクにはバッテリ端子をのぞける小窓があるが，奥側（プラス）側の端子が届かない（というか暗すぎて見えん）．説明書を見ると「ヒューズ交換時にはトランク内張を外せる」とあるので試してみる．内張は差し込み式のプラスチックピンで固定してあり，これを外すには小さめのマイナスドライバーでピンの内側の円の側部をテコの要領で上に押し上げる必要がある（説明書の図が雑すぎて構造が分からず，このピンの解除にかなり手こずった．一個外せば要領が分かるので後は簡単）

内張は右側だけを外せばあとは上に押し上げるとバッテリ上部に手が届く．やれやれ．

手順

まずプラス側の端子カバーを上に押し上げ，赤色のプラグを接続する（といってもカバー上部がフレームに干渉するため，斜めに突っ込んでギリギリ保持できる程度にしか接続できない）．
マイナス端子に黒いプラグをつなぐ．するとテスターの電源が入る．エンジンがかかっている状態だと14Vなのでバッテリーテストは実行できず，充電テストしか実行できない．
バッテリー残量を知りたい場合はエンジンを切り，テストを実行する．ACC規格がわからないのでJISを選択して，搭載バッテリーの型番を選ぶ（うちの車はディーラーの記録から46B24らしい）．
テストの実行結果を確認，健全度と充電容量が表示されるので，写真などで結果を控えておく．

精度

テスターの精度はいまいち分からない．最初に計ったとき16%と出てこれはやばいと思い，次に計ったら6%だったので終わったと思ったら，次は58%になり，20kmほど走らせた後計測したら85%だった．ある程度暖まらないとまともに計測でできないのか，プラス端子側のプラグ固定が不安定なせいなのか，バッテリーに対応していないのかは，よく分からない．

今後も測定を繰り返して安定性を確認するしかないかな...

2023-01-29

MSIマザボでの起動ドライブの設定

特定のアプリケーションのエラーでwindwosが強制終了し場合に，再起動すると「起動ドライブ内にOSが見つからない」というエラーがしばしば発生する．以下，その対処方法．

MSIのロゴが出たらf11キーを連打してUFEIを起動する．
起動ドライブの優先順序が表示されるが，この状態ではOSがインストールされたドライブがない．
アドバンストモード，→boot settingsへと進み，画面を下にスクロールすると，ドライブ一覧が表示され，起動順序を指定できる（このときタブをスクロールできることに気づかず，指定場所を何度も探した）．
OSがインストールされたドライブを一番最初に指定する（OSの入っていないドライブが勝手に最初に指定されることがエラーの原因）

起動ディスクはM.2接続のSSDだが，他にもM.2接続しているSSDがあるため，どうやらそちらを起動ディスクに設定してしまうらしい．

2022-09-28

TeXlive更新

久しぶりにTeXLiIveをアップデートした． 3日後に忘れるのでメモを残しておく．

TeXStudioの設定

設定→コマンド→latexを platex.exe -synctex=1 -interaction=nonstopmode %.tex に変更する．この変更で日本語文書クラスのコンパイルができる．

beamerでminjisがないと怒られたら，jsclassesをtexliveshellでアップデートもしくは，そもそも入ってなければインストールしておく．

TeXStudioのショートカットキー設定

とりあえず現環境から

仮空きスペース
\vspace{2cm}
%|
\vspace{2cm}

dots:
$\cdots \cdots$

引用:
\begin{quote}%\ %\\end{quote} 

インライン数式 $...$:
$%|$

alignat:
\begin{alignat*}{3}%\ & = & \quad  &  \\ %\ & = &  &  \\ %\ & = &  &  \\ %\ &= & & %\ \end{alignat*}

複数行数式:
\begin{align*} %\ & = \\ %\  & = %\ \end{align*}

ディスプレイ数式 \[...\]:
\[ %| \]

\begin{verbatim:
\begin{verbatim}%\ %\ \end{verbatim}

などをコピペして登録．こういう設定，エクスポートできそうな気もしますが，頻繁にやることでもないので手作業で設定します．

ビルド＆表示は txs:///dvi-pdf-chain

とりあえずこれで，コンパイルはOK．

2022-09-26

ラウンド毎に異なるパラメータを設定する方法

前回の記事でワンショットの剥奪ゲームのトリートメントを作成しました．同じ被験者にゲームの設定を変えつつ複数回繰り返す方法を考えてみましたが3日後には忘れそうなので，例によってメモを残しておきます．

グループに共通の値を設定する

結論から言うと，次のコードで一応の目的は達せられます．

class C(BaseConstants):
    NAME_IN_URL = 'relative_deprivation'
    PLAYERS_PER_GROUP = 5
    NUM_ROUNDS = 2
    Rvec = [3,5]#セッション毎に変わるグループ共通の利益率
    nvec = [2,4]#セッション毎に変わるグループ共通の成功者数

class Subsession(BaseSubsession):
    pass

class Group(BaseGroup):
    x = models.IntegerField()
    R = models.IntegerField()#セッション毎に変わるグループ共通の利益率
    PROFIT = models.IntegerField()
    n = models.IntegerField()

def set_parameters(group: Group):#グループオブジェクトを引数にした関数，引数groupのタイプ（クラス）がGroupであることを指定
    group.R = C.Rvec[group.round_number-1]
    group.PROFIT = group.R*100-100
    group.n = C.nvec[group.round_number-1]

class MyWaitPage(WaitPage):
    after_all_players_arrive = set_parameters

このような面倒なコードになる理由は

ズバリ

Conventionally, group-level treatments are assigned in creating_session:

for g in subsession.get_groups():
     g.treatment = random.choice([True, False])

However, this doesn't work when using group_by_arrival_time, because groups are not determined until players arrive at the wait page. (All players are in the same group initially.) Instead, you need to assign treatments in after_all_players_arrive.

だからです．

この実験では行為者の選択前に所与の条件（グループ固定の属性である利益率と成功人数）を、表示する必要があります．条件表示の前にcreating_sessionを使ってグループオブジェクトの属性に値を渡すとエラーがでるというわけです．そこで条件表示の前にafter_all_players_arrive = set_parametersを使ってグループ属性を定義します．（sessionレベルの変数の定義しようと試行錯誤しましたが，うまくいきませんでした）．

結局，

定数クラスで水準のベクトルを準備する
set_parameters関数を定義し，group属性のround_numberをつかって，ラウンドに対応したパラメータを取得
最初の入力フォーム表示前にwaitpageを追加してafter_all_players_arrive = set_parametersでパラメータをグループの属性に代入

という手順でラウンド毎のパラメータ設定を実行しました．やれやれ．

多少面倒ですが，otreeは「手軽にオブジェクトの属性やメソッドを定義できるPythonの特徴」を生かした，便利なパッケージと言えるでしょう．

The One with ...

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