# 数理社会学・なぜタバコがやめられないか

01章：なぜタバコがやめられないか
モデルの説明：
べッカーの依存症モデル

```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)
```

```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)
```