#install.packages(“msm”)#used for truncated normal dist when generating data
#install.packages(“rstan”)
require(rstan)
require(msm)

######
generate.data<-function(c=1, k=1, s0=10, e0=100, t=1:500){
z=c*e0^(-k/(c+k))*s0^(-c/(c+k))*t - (e0/s0)^(c/(c+k)) + (c/k)*(s0/e0)^(k/(c+k))

p=matrix(nrow=length(z))
for(i in 1:length(z)){
  f<-function(x) (x/(1-x))^(k/(c+k))- ((1-x)/x)^(c/(c+k))-z[i]
  p[i]<-uniroot(f, interval= c(0,1), tol=10^-9)$root
  p[i]<-rtnorm(n=1,mean=p[i], sd=.05, lower=0, upper=1)
}
return(p)
}

subj1<-generate.data(c=1, k=1, s0=10, e0=100, t=1:500)
subj2<-generate.data(c=.5,k=.5, s0=10, e0=100, t=1:500)
subj3<-generate.data(c=.3,k=.3, s0=10, e0=100, t=1:500)
subj4<-generate.data(c=.25,k=.25, s0=10, e0=100, t=1:500)

dat<-rbind(subj1,subj2,subj3,subj4)
dat<-cbind(rep(1:length(subj1),4),dat)
dat<-cbind(rep(1:4,each=length(subj1)),dat)
colnames(dat)<-c("Subj","Time","y")

####


Ndata=nrow(dat)
Nsubj=length(unique(dat[,"Subj"]))
subj=dat[,"Subj"]
y=dat[,"y"]
t=dat[,"Time"]
dat2<- list(Ndata=Ndata, Nsubj=Nsubj, subj=subj, y=y, t = t)


####
STAN_code <- '
data {
int<lower=0> Ndata;
int<lower=0> Nsubj;
int<lower=0> subj[Ndata];
real y[Ndata];
int<lower=0> t[Ndata];
}
parameters {
real<lower=1, upper=1000> s0;
real<lower=1, upper=1000> e0;
real<lower=0, upper=1> k[Nsubj];
}
model {
real z[Ndata];
real p[Ndata];

for(i in 1:Ndata){
z[i]<-(k[subj[i]]*t[i] +s0 -e0)/(s0*e0)^.5;
p[i]<-.5*(z[i]/((z[i])^2+4)^.5 +1);
}
y ~ normal(p , .05 );

}
'

####
fit <- stan(model_code = STAN_code, data = dat2, iter = 10000, chains = 3, warmup=500)
fit2<-extract(fit)
#plot(fit)
#traceplot(fit,ask=T)


####Make Plots
par(mar=c(4.1,4.1,2.1,2.1))
par(mfrow=c(Nsubj+1,2))
s0<-mean(fit2$s0)
e0<-mean(fit2$e0)
hist(fit2$s0, main="s0", xlab="")
hist(fit2$e0, main="e0", xlab="")
for(s in 1:Nsubj){
k<-mean(fit2$k[,s])
c=k;
t=1:500
z=c*e0^(-k/(c+k))*s0^(-c/(c+k))*t - (e0/s0)^(c/(c+k)) + (c/k)*(s0/e0)^(k/(c+k))
p=matrix(nrow=length(z))
for(i in 1:length(z)){
  f<-function(x) (x/(1-x))^(k/(c+k))- ((1-x)/x)^(c/(c+k))-z[i]
  p[i]<-uniroot(f, interval= c(0,1), tol=10^-9)$root
}

hist(fit2$k[,s], main=paste("k",s), xlab="")
plot(t,dat[which(dat[,"Subj"]==s),"y"], ylim=c(0,1),
ylab="p", xlab="t", main=paste("Subj",s))
lines(t,p, lwd=4, col="Red")
}
par(mar=c(5.1,4.1,4.1,2.1))
####