遅延価値関数と階層ベイズを用いた
男心をくすぐる女の戦略.R


遅延価値観数と階層ベイズを用いた男心をくすぐる女の戦略.R from MrUnadon
author: Mr.Unadon (見習い飯炊き兵)
動作環境:Mac OS Sierra 10.12.1; R version3.3.1; rstan 2.10.1



はじめに


本記事は第60回Tokyo.R のLTセッションへのエントリーです。

発表スライドはページ上部を御覧ください。

スライドタイトルの「遅延価値観数」は誤字で、正しくは「遅延価値関数」です。すみません。


内容への補足: コード類の詳細


ここでは、LTでは紹介しきれなかったStanコードを紹介します。

データは2つの選択肢からワクワクする方を選ぶ。選択のデータです。

それを補足する情報として、モデルさんの魅力評定得点、デートまでの遅延日数の条件データがあります。

これらのデータを用いて、「時間で変化するデートの魅力」を推定します。

まずは、データを確認しておきましょう。

######################################################
#  Tokyo.R 2017/04/22
#  LT presenter: Mr.Unadon
#  『男心をくすぐる女の戦略.R』
######################################################
invisible({rm(list=ls());gc();gc()})

library(tidyverse)#データ操作
library(rstan)#HMCベイズ推定
library(formattable)#HTMLのtable

#Stanのチェーンを複数コアで並列化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#実験データ一括読み込み
Dat<-do.call(rbind,tmp<-list.files()%>%lapply(read.csv))%>%
    mutate(trial=rep(c(1:nrow(tmp[[1]])),length(tmp)),ID=rep(1:length(tmp),each=nrow(tmp[[1]])))%>%
    select(ID,Choice1Delay,NowEval,DelayEval,DelayTime)

formattable(Dat[1:100,])
ロングデータ形式のデータになっています。
Table1. 実験データ(抜粋した参加者の順序は匿名性を守るためランダム化してます)

Stanに渡すデータ


続いて、Stanに渡すデータのリストとパラメータを宣言します。 最も関心のあるパラメータは、Vt(将来のデートの魅力値)です。

#Stanデータ
datastan<-list(meanA=mean(Dat$NowEval),#Idol Attractive mean
               DelayTime=Dat$DelayTime, #デートまでの日数
               DelayEval=Dat$DelayEval, #遅延選択肢のモデルさんの魅力値
               NowEval=Dat$NowEval,#即時選択肢のモデルさんの魅力値
               N=max(Dat$ID), #参加者数
               r=length(Dat$ID), #行数
               c=length(Dat), #列数
               ID=Dat$ID, #参加者ID
               choose=Dat$Choice1Delay)#選択データ。遅延選択肢が1。

#Stanパラメータ
pars=c("beta", #softmax coefficients
       "gamma", #discount rate
       "mu", #hierarchical mean
       "sigma", #hierarchical SD
       "theta", #prob choose later
       "Vt") #values of later dating


Stanコード


さて、Stanコードを書いていきましょう。

スライドとパラメータの名前が対応していなくて申し訳ないですが、

選択データ(choose)をベルヌイ分布で予測し、その選択確率パラメータthetaをソフトマックス行動選択で定義しています。

ソフトマックスに入る二つの選択肢の価値の片方に将来のデートの魅力(Vt)が入っており、

Vtをわくわく方程式で定義しています。

model<-"
data {
    int N; //number of subjects
    int r; //nrow of data set
    int c; //ncol of data set
    int choose[r]; //choose delay =1
    int ID[r]; //integer ID
    real NowEval[r]; //Attractive sooner idol
    real DelayEval[r]; //Attractive later idol
    int DelayTime[r]; //wating time(day)
    real meanA;
}
parameters {
    real <lower=0, upper=1> beta [N]; //softmax coefficients
    real <lower=0, upper=1> gamma[N]; //discount rate
    vector <lower=0,upper=1>[2]mu;
    vector <lower=0>[2] sigma;
}

transformed parameters{
    real <lower=0,upper=1>theta[r];
    real Vt[r];
    for (i in 1:r){
        Vt[i] = DelayEval[i]+(DelayEval[i]*(gamma[ID[i]]^DelayTime[i])-meanA);
        theta[i] = 1/(1+exp(beta[ID[i]]*(NowEval[i]-Vt[i])));
    }
}

model{
    for (i in 1:r){
        choose[i] ~bernoulli(theta[i]);
    }
    for (j in 1:N){
        beta[j]  ~ normal(mu[1],sigma[1]);
        gamma[j] ~ normal(mu[2],sigma[2]);
    }
    mu ~ normal(0,100);
    sigma ~ cauchy(0,2.5);
}
"

#stan fit
fit<-stan(model_code = model,data=datastan,pars=pars,
              iter = 5000,warmup=2000,chains = 3,thin = 3)


ここまでで、 rstanを用いた推定が完了しました。


Stanの推定結果を表で確認


推定結果を確認していきます。まずは、収束とパラメータ推定値を表でまとめるところから参ります。

#
# 推定結果の整理
#

#Diagnose でGerman-Rubin testの確認
library(shinystan)
launch_shinystan(fit)

#formattableで最小実効サンプルサイズ等の確認
df<-data.frame(summary(fit)$summary)
#RhatのたまにあるNaNは数値認識されないのでゼロで置き換え
df$Rhat[is.nan(df$Rhat)]<-0
#小数点第3ケタまでにする
for(i in 1 : length(df)){df[,i]<-digits(df[,i],3)}
#対数尤度の項lp___を取り除く
df<-df[-length(df$mean),]

C<-c("steelblue","tomato")
#結果を表で描画
p<-formattable::formattable(df[1:100], list(
    mean = formattable::formatter("span",style = x ~ formattable::style(color=ifelse(x <0, "steelblue","tomato"))),
    sd = formattable::color_bar("orange"),
    X2.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    X97.5. = formattable::formatter("span",style = x ~ ifelse(x <0, formattable::style(color = "steelblue", font.weight = "bold"),NA)),
    n_eff = formattable::formatter("span",style = x ~ formattable::style(color = ifelse(rank(-x) <= c(length(df$n_eff)-1), "gray", "tomato")),
                                   x ~ sprintf("%.2f (rank: %g)", x, rank(-x))),
    Rhat = formattable::formatter("span", x ~digits(x,2),
                                  style = x ~ formattable::style(color = ifelse(x>=1.1, "tomato", "green")),
                                  x ~ icontext(ifelse(x<1.1, "ok", "remove"), ifelse(x<1.1, "Yes", "No")))
))

p

ふむふむ。良い感じだと思います。
Table2. stanの推定結果

結果を図で確認


続いて、実験結果と解析結果をそれぞれ可視化していきます。

まずは、後日のデートを選んだ比率から。

#------------------------------------------------------------------------
#データの描画
#------------------------------------------------------------------------

vt<-rowMeans(t(rstan::extract(fit)$Vt))
df.plot<-Dat%>%dplyr::mutate(vt=vt)%>%
    dplyr::group_by(DelayTime)%>%
    dplyr::mutate(probChooseLater=mean(Choice1Delay))%>%
    dplyr::mutate(vtm=mean(vt))%>%
               ungroup()
               
#prob
gProb<-ggplot(df.plot,aes(x=DelayTime,y=probChooseLater))+
    scale_x_continuous(expand=c(0.01,0),breaks = c(1,3,7,13,30),labels = c("1day","3day","7day","13day","30day"))+
    scale_y_continuous(limits = c(0.25,0.55))+
    geom_hline(yintercept = 0.5,linetype=2)+
    geom_point(colour="darkblue",size=2)+
    geom_line(colour="darkblue",size=1.1)+
    theme(plot.subtitle = element_text(vjust = 1), 
    plot.caption = element_text(vjust = 1), 
    panel.grid.major = element_line(colour = "gray94"), 
    panel.grid.minor = element_line(colour = "gray94", 
        linetype = "dashed"), axis.title = element_text(family = "serif", 
        size = 18), axis.text = element_text(family = "serif", 
        size = 16, hjust = 1), axis.text.x = element_text(angle = 45), 
    plot.title = element_text(family = "serif", 
        size = 22), panel.background = element_rect(fill = "white")) +labs(title = "後日のデート選択率", 
    x = "デート日", y = "後日の選択率")

#描画
plot(gProb)

遅延日数別、後日のデートを選ぶ選択率がこちら。

DatingResult1.png

続いて、時間で変化するデートの価値をグラフ化します。 コードが回りくどいですがすみません、おつかれもーどです。

#-----------------------------------------------------------------------
#時間で変化するデートの魅力
#-----------------------------------------------------------------------

df2<-df.plot%>%dplyr::select(DelayTime,vtm)%>%
    dplyr:: group_by(DelayTime)%>%
    dplyr:: mutate(vtmm=mean(vtm))%>%
    dplyr::ungroup()%>%
    dplyr::distinct(DelayTime,.keep_all=T)%>%
    dplyr::arrange(DelayTime)%>%
    dplyr::select(vtmm)
df2<-as.data.frame(df2)
mVal<-mean(Dat$NowEval)
set<-rbind(mVal,df2)
delay<-c(0,1,3,7,13,30)
df3<-data.frame(val=set,delay=delay)

gValue<-ggplot(df3,aes(x=delay,y=vtmm))+theme_bw()+
    scale_y_continuous(limits = c(20,65))+
    scale_x_continuous(expand=c(0.01,0),breaks = c(0,1,3,7,13,30),labels = c("Today","1day","3day","5day","7day","30day"))+
    geom_hline(yintercept = df3$vtmm[1],linetype=2)+
    geom_point(colour="#e4007f",size=2)+
    geom_line(colour="#e4007f",size=1)+
    theme(plot.subtitle = element_text(vjust = 1), 
    plot.caption = element_text(vjust = 1), 
    panel.grid.major = element_line(colour = "gray94"), 
    panel.grid.minor = element_line(colour = "gray88", 
        linetype = "blank"), axis.title = element_text(family = "serif", 
        size = 16), axis.text = element_text(family = "serif", 
        size = 14, hjust = 1), axis.text.x = element_text(angle = 45), 
    plot.title = element_text(family = "serif", 
        size = 22), panel.background = element_rect(fill = "white")) +labs(title = "女性とデートする魅力の時間的変化", 
    x = "デート日", y = "主観的な魅力")

#描画
plot(gValue)

よいしょ、 時間で変化するデートの魅力がこちらです。

DatingResult2.png


最後に


色んな数理モデルを考えて、データ生成メカニズムを考えるのもたのしいですね。

デートの魅力は明日がmaxなので(すくなくとも私の周囲の飢えた男子に限りますが)、

女性の皆様は、男性を明日のデートに誘いましょう。

Enjoy!!

Written on April 30, 2017