社会を＜モデル＞でみる：数理社会学への招待
01章：なぜタバコがやめられないか
モデルの説明：
べッカーの依存症モデル
効用関数
仮定１：今日の消費量x(t)は、昨日の消費量x(t-1)に依存する
仮定２：時間が経てば魅力が減る(0<=C<=1) 仮定３：満足度はu*(x(t),C*x(t-1))で表し、依存状態はもっとも高い効用uを得るとする 仮定４：はじめはあまり依存せず、途中にはまりだし、最後には慢性的にはまる つまり、昨日の消費が増えるにしたがって最大の満足を得るには今日もたくさん消費する必要がある。しかし、消費が増えると魅力が減りそれほど増えなくなる。 ＃実装例（ベッカーのモデルを忠実には再現できていない） 今日の消費量をy、昨日の消費量をxとした場合、それらの関係を以下の式で表す。 y=K/(1+exp(r*K*(shift-x))) rをgainと呼び、今日と昨日の消費量の変化率である（rが大きいと昨日と同じ満足度を得るにはより多くの消費が必要） Kは最大消費量（あまり意味なない） shiftは今日の消費量に対する抑制・ペナルティー効果

smoker<-function(K=20,r=0.02,shift=9,nDays=10) {
  #------  input args -------------------
  #K: max # of consumption 
  #r:  satisfaction decay  
  #shift: shift in satisfaction function
  #initial amount of consumption = 10
  #-----------------------------------------
  x.seq=seq(0,K,0.01);
  y.seq=K/(1+exp(r*K*(shift-x.seq)))
  plot(x.seq,y.seq,type="l",lwd=2,xlab=c("consumption at T"),
   ylab=c("consumption at T+1"),xlim=c(-1,21),ylim=c(-1,21))
  lines(x.seq,x.seq,type="l",col="blue",lty="dashed",lwd=2)
  x=10;
  y=K/(1+exp(r*K*(shift-x)))
  lines(c(x,x),c(0,y),type="l",col="red",lwd=2)
  lines(c(x,y),c(y,y),type="l",col="red",lwd=2)
  for (i_days in 2:nDays) {
    x.old=x;x=y;
    y=K/(1+exp(r*K*(shift-x)))
    lines(c(x,x),c(x,y),type="l",col="red",lwd=2)
    lines(c(x,y),c(y,y),type="l",col="red",lwd=2)
  }
}
# example 
par(mfrow=c(2,2))
smoker(K=20,r=0.02,shift=9,nDays=10)
smoker(K=20,r=0.01,shift=9,nDays=10)
smoker(K=20,r=0.02,shift=11,nDays=10)
smoker(K=20,r=0.01,shift=10.5,nDays=10)

実装例　その2 (2015年度）

addiction<-function(s,r,x){
# initialization
k=100
xs=seq(0,k,0.1);
plot(xs,k/(1+exp(k*r*(s-xs))),type='l',lwd=2,xlim=c(-3,k+5),ylim=c(-3,k+5),
 xlab="consumption @ time T",ylab="consumptino @ time T+1")
lines(c(0,k),c(0,k),col='blue',lty=3,lwd=3)

# main
for (i_day in 1:14) {
  u=k/(1+exp(k*r*(s-x)))
  if (i_day ==1){
    lines(c(x,x),c(5,u),col='red',lty=2,lwd=2)
  } else { lines(c(x,x),c(x,u),col='red',lty=2,lwd=2)}
  lines(c(x,u),c(u,u),col='red',lty=2,lwd=2)
  lines(x,u,type='p',cex=2.5,pch=20,col="magenta")
  text(x,0,'|',cex=1.5,pch=18,col="darkgreen",font=2)
  x=u
 }
}

> addiction(20,0.0005,30)