## Simulation of the variable slopes model.

## Structural Parameters

I <- 100                                #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
## 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
act_target <- 40
act_slope <- 2
## policy funciton:  probability dose will be given.
act_pol <- function(obs,thetahat=NULL) {
  1-pnorm(obs,act_target,act_slope)
}

## 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,])*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)
}