In [1]:

    

t.var = seq(1,10, by = 1)


    

In [2]:

    

alpha = c(0, .1, .5, .9); beta = .2;
parameters = data.frame(alpha = alpha, beta = beta)


    

In [3]:

    

options(repr.plot.width = 5)
options(repr.plot.height = 5)

options(warn=-1) ## OFF 경고
suppressMessages(library(dplyr))


    

In [4]:

    

## 기억할 확률 = exp(-1*alpha*t.var)+beta
res.df = parameters %>% group_by(alpha) %>% do(res = exp(-1*.$alpha*t.var)+.$beta)


    

In [34]:

    

par(family = 'AppleGothic')
par(mfrow = c(2,2))

for (i in 1:nrow(res.df))
{
    res.df$res[[i]][res.df$res[[i]] > 1] = 1 ## 1 이상 1로
    plot(t.var, res.df$res[[i]], xlab = '시간', ylab = '기억 (기대) 확률', main = paste0('alpha = ', res.df$alpha[[i]]), type = 'l')
}


    

In [6]:

    

hist(rbeta(1000000,1,1))


    

In [7]:

    

suppressMessages(library(rstan))
#rstan_options(auto_write = TRUE)
#options(mc.cores = parallel::detectCores())


    

In [8]:

    

model <- "
// Retention With No Individual Differences
data { 
int ns;
int nt;
int k[ns - 1,nt - 1];
int t[nt];
int n;
}
parameters {
real<lower=0,upper=1> alpha;
real<lower=0,upper=1> beta;
} 
transformed parameters {
matrix<lower=0,upper=1>[ns,nt] theta;

// Retention Rate At Each Lag For Each Subject Decays Exponentially
for (i in 1:ns) {
for (j in 1:nt) {
theta[i,j] <- fmin(1, exp(-alpha * t[j]) + beta);
}
}
}
model {
// Priors
alpha ~ beta(1, 1);  // can be removed
beta ~ beta(1, 1);  // can be removed

// Observed Data
for (i in 1:(ns - 1))
for (j in 1:(nt - 1))
k[i,j] ~ binomial(n, theta[i,j]);
}
generated quantities {
int<lower=0,upper=n> predk[ns,nt];

// Predicted Data
for (i in 1:ns)
for (j in 1:nt)
predk[i,j] <- binomial_rng(n, theta[i,j]);
}"


    

In [9]:

    

## Data
t     <- c(1, 2, 4, 7, 12, 21, 35, 59, 99, 200)
nt    <- length(t)
slist <- 1:4
ns    <- length(slist)

k1 <- matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4, NA,
               17, 13,  9,  6, 4, 4, 4, 4, 4, NA,
               14, 10,  6,  4, 4, 4, 4, 4, 4, NA,
               NA, NA, NA, NA,NA,NA,NA,NA,NA, NA), nrow=ns, ncol=nt, byrow=T)

k <- k1[1:(ns - 1), 1:(nt - 1)]   # Excluding NAs (for Stan solution)

n <- 18

data <- list(k=k, n=n, t=t, ns=ns, nt=nt) # To be passed on to Stan


    

In [10]:

    

data


    

    
Out[10]:





	
$k
		


	
18
18
16
13
 9
 6
 4
 4
 4
	
17
13
 9
 6
 4
 4
 4
 4
 4
	
14
10
 6
 4
 4
 4
 4
 4
 4



	
$n
		
18
	
$t
		

	
1
	
2
	
4
	
7
	
12
	
21
	
35
	
59
	
99
	
200


	
$ns
		
4
	
$nt
		
10

In [11]:

    

myinits <- list(
    list(alpha=.5, beta=.1))

parameters <- c("alpha", "beta", "predk")  # Parameters to be monitored


    

In [12]:

    

samples <- stan(model_code=model,   
                data=data, 
                init=myinits,  # If not specified, gives random inits
                pars=parameters,
                iter=20000, 
                chains=1, 
                thin=5,
                warmup = 5000,  # Stands for burn-in; Default = iter/2
                seed = 12345  # Setting seed; Default is random seed
)


    

    




SAMPLING FOR MODEL '7d2b1c02fbf3ff8ebac7c6b5ce5600a1' NOW (CHAIN 1).

Chain 1, Iteration:     1 / 20000 [  0%]  (Warmup)
Chain 1, Iteration:  2000 / 20000 [ 10%]  (Warmup)
Chain 1, Iteration:  4000 / 20000 [ 20%]  (Warmup)
Chain 1, Iteration:  5001 / 20000 [ 25%]  (Sampling)
Chain 1, Iteration:  7000 / 20000 [ 35%]  (Sampling)
Chain 1, Iteration:  9000 / 20000 [ 45%]  (Sampling)
Chain 1, Iteration: 11000 / 20000 [ 55%]  (Sampling)
Chain 1, Iteration: 13000 / 20000 [ 65%]  (Sampling)
Chain 1, Iteration: 15000 / 20000 [ 75%]  (Sampling)
Chain 1, Iteration: 17000 / 20000 [ 85%]  (Sampling)
Chain 1, Iteration: 19000 / 20000 [ 95%]  (Sampling)
Chain 1, Iteration: 20000 / 20000 [100%]  (Sampling)
#  Elapsed Time: 0.290048 seconds (Warm-up)
#                1.04542 seconds (Sampling)
#                1.33547 seconds (Total)

In [13]:

    

##Figure 10.2

alpha <- extract(samples)$alpha
beta <- extract(samples)$beta
d.beta <- density(beta)

layout(matrix(c(1,2,3,0),2,2,byrow=T), width=c(2/3, 1/3), heights=c(2/3,1/3))
#layout.show()

par(mar=c(2,2,1,0))
plot(alpha,beta, xlab="", ylab="", pch = 3, xlim=c(0,1),ylim=c(0,1),axes=F)
box(lty=1)

par(mar=c(2,2,1,4))
plot(d.beta$y, d.beta$x, ylim=range(c(0,1)), xlim=rev(range(d.beta$y)), 
     type='l', axes=F, xlab="", ylab="")
axis(4, at=c(0,1))
mtext('Baseline', side=4,line=1, cex=1.3)
box(lty=1)

par(mar=c(6,2,0,0))
plot(density(alpha),zero.line=F, main="", ylab="", xlab="", cex.lab=1.3,
     xlim=c(0,1), axes=F)
axis(1,at=c(0,1))
mtext('Decay Rate', side=1.2,line=1, cex=1.3)
box(lty=1)


    

In [14]:

    

##Figure 10.3

layout(matrix(c(1:4),2,2,byrow=T))
#layout.show()
sc <- 3.5
jj <- numeric()
xx <- numeric()

for (i in 1:ns) 
{
    plot(-1,100,xlim=c(0,10),ylim=c(0,18), main=(paste("Subject", i)),
         xlab=("Time Lags"), ylab=("Retention Count"),cex.lab=1.3, axes=F)
    axis(1, at=c(1,2,3,4,5,6,7,8,9,10), 
         lab=c("1","2","4","7","12","21","35","59","99","200"),cex.axis=0.7)
    axis(2, at=c(0,18),lab=c("0","18"),cex.axis=0.7)
    box(lty=1)
    for (j in 1:nt) 
    {
        count <- hist(extract(samples)$predk[,i,j],c(0:n),plot=F)
        count <- count$counts
        count <- count/sum(count)
        for (x in 1:n)
        {
            if (count[x]>0)
            {
                points(j,x,pch=22, col="black",cex=sc*sqrt(count[x]))
                if (!is.na(k1[i,j]) && k1[i,j]==x)
                {
                    points(j,x,pch=22,bg="black",cex=sc*sqrt(count[x]))
                    jj <- c(jj,j)
                    xx <- c(xx,x)
                }
            }
        }
    }
    coords <- list(x=jj, y=xx)
    lines(coords,lwd=2)
    jj <- numeric()
    xx <- numeric()
}