## Simulation of the variable slopes model.

## Structural Parameters

I <- 400                                #Number of students
Tmax <- 10                              #Maximum time periods

## Experimental Design:  Time and dosage patterns.
## For the moment assume all have the same number of measurements.
Time <- matrix(1/Tmax,Tmax-1,I)         #Equally spaced measurements
cumTime <- rbind(0,apply(Time,2,cumsum))
## vector[I] Time[Tmax-1];
Dose <- matrix(0,Tmax-1,I)              #Start with 0, we will fill
## vector[I] Dose[Tmax-1];
                                        #this in later.
T <- rep(Tmax,I)                        #Number of measurement points
                                        #for each person
## Place holder for data
Obs <- matrix(NA,Tmax,I)
## vector[I] Obs[Tmax-1];

## Hyperparameters
slope_mu_mu <- 1                        #Mean slope
slope_mu_std <- .1

## Evidence Model priors
obs_rel <- .8
obs_mean_t1 <- 50
obs_std_t1 <- 10
res_std <- obs_std_t1*sqrt(1-obs_rel)
obs_slope <- obs_std_t1*sqrt(obs_rel)
obs_int <- obs_mean_t1
## Action assignment function parameters
## Two changes here.  First, now a hard cut rather than a soft one.
## Second, we assume that the teacher who has perfect insight
## overrides a certain percentage of the time.
## Third, we assume the target moves over time.

act_target <- -1
act_target_slope <- 1
override_p <- 0
## policy funciton:  probability dose will be given.
act_pol <- function(obs,theta,cumtime) {
  target <- act_target+act_target_slope*cumtime
  thetahat <- (obs-obs_int)/obs_slope
  thetahat <- ifelse(runif(I) < override_p, theta, thetahat)
  as.integer(thetahat < target)
}

## Parameters:  Brownian Motion
var_innov <- .1                          #Variance of the innovations
treat_eff <- .25                        #Treatment effect
slope_std <- .33
slope_r2 <- .36
## Level 2 parameters
slope_mu <- rnorm(1,slope_mu_mu,slope_mu_std)

## Level 1 parameters
slope_res <- rnorm(I,0.0,1.0)
T0_val <- rnorm(I,0,1.0)
innov <- matrix(NA,Tmax-1,I)
## vector[I] innov[Tmax-1];
slope <- slope_mu +
  slope_std*(sqrt(1-slope_r2)*slope_res + slope_r2*T0_val)



theta <- matrix(NA,Tmax,I)
## vector[I] theta[Tmax];
theta[1,] <- T0_val
Obs[1,] <- rnorm(I,obs_int+obs_slope*theta[1,],res_std)
for (t in 2:Tmax) {
  ## Dosage based on observation.
  Dose[t-1,] <-
    rbinom(I,1,act_pol(Obs[t-1,],theta[t-1,],cumTime[t-1,]))*
    Time[t-1,]
  innov[t-1,] <- rnorm(I,0,sqrt(var_innov*Time[t-1,]))
  theta[t,] <- theta[t-1,] + slope*Time[t-1,] +
    treat_eff*Dose[t-1,]+innov[t-1,]

  Obs[t,] <- rnorm(I,obs_int+obs_slope*theta[t,],res_std)
}


override_p <- .5
theta5 <- matrix(NA,Tmax,I)
Obs5 <- matrix(NA,Tmax,I)
Dose5 <- matrix(0,Tmax-1,I)              #Start with 0, we will fill

theta5[1,] <- T0_val
Obs5[1,] <- rnorm(I,obs_int+obs_slope*theta5[1,],res_std)
for (t in 2:Tmax) {
  ## Dosage based on observation.
  Dose5[t-1,] <-
    rbinom(I,1,act_pol(Obs5[t-1,],theta5[t-1,],cumTime[t-1,]))*
    Time[t-1,]
  ## Reuse innovations
  ## innov[t-1,] <- rnorm(I,0,sqrt(var_innov*Time[t-1,]))
  theta5[t,] <- theta5[t-1,] + slope*Time[t-1,] +
    treat_eff*Dose5[t-1,]+innov[t-1,]

  Obs5[t,] <- rnorm(I,obs_int+obs_slope*theta5[t,],res_std)
}