# These are files for the data generation routines. ## Regression coefficients for multivariate regression. evalA <- matrix(c(.7,.3,.3,.7),2,2) ## Noise std for multivariate regression. evalSig <- c(.65,.65) ## This function takes the proficiency expressed as a vector ## (Mechanics, Fluency) and generates a random lesson value. evalLesson <- function(proficiency) { result <- proficiency %*% evalA + rnorm(length(evalSig))*evalSig result <- round(2*result)/2 result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) } ## Calculates likelihood of observation given proficiency. ## Assume Proficiency is a matrix with columns equal to the ## particles. Assume observation is a vector (Mechanics, ## Fluency) Returns a vector of likelihoods, one for each ## particle evalLike <- function (proficiency, observation) { if (!is.matrix(proficiency)) proficiency <- matrix(proficiency,1) # I'm ignoring the normalization constant and the rounding # error ker <- sweep(sweep(proficiency %*% evalA,2,observation), 2,evalSig,"/") exp(-apply(ker^2,1,sum)/2) #exp(-sum( ( (proficiency %*% evalA - observation)/evalSig)^2)/2) } ## Constant noise term in noisy-or model q0 <- c(Mechanics=.8,Fluency=.9) ## Upper and lower bounds for delta in noisy-or model threshlow <- c(Mechanics=-1,Fluency=-2) threshhigh <- c(Mechanics=2,Fluency=1) ## skill bypass parameters in noisy-or model q1 <- c(Mechanics=.8,Fluency=.7) ## NOTE: These values are overridden for Simulation 2 (used ## in the paper). # Calculates success probability of action based on current # proficiency level. # Assumes action is a vectors of the form (Mechanics,Fluency) # Assumes proficiency is a matrix with columns (Mechanics, # Fluency) returns a vector of probabilities of success for # the lesson. actionSuccessP <- function (proficiency, action) { if (!is.matrix(proficiency)) { proficiency <- matrix(proficiency,1) } ## delta <- threshlow <= action-proficiency & ## action - proficiency <= threshhigh diff <- sweep(proficiency,2,action) delta <- sweep(diff,2,threshlow,">=") & sweep(diff,2,threshhigh,"<=") ## q0 * prod(q1^(1-delta)) ## Need to transpose matrix to get rows varying fastest qq1 <- apply(q1^(1-t(delta)),2,prod) outer(qq1,q0,"*") } ## Constant term in learning rate eta <- c(Mechanics=.2,Fluency=.2) ## Zone of Proximal Development term in learning rate gamma <- c(Mechanics=.1,Fluency=.15) ## Calculates the increase rate based on current position ## and selected action. ## Assumes proficiency is a matrix and action is a vector learnRate <- function(proficiency,action) { if (!is.matrix(proficiency)) { proficiency <- matrix(proficiency,1) } # eta - gamma*(proficiency-action)^2 lam <- sweep(proficiency,2,action)^2 sweep(-sweep(lam,2,gamma,"*"),2,eta,"+") } ## Constant parameters of the learning process forgetRate <- c(Mechanics=.02,Fluency=.01) ## R(t) = (rTimediff*success + rTimeconst)*t rTimediff <- .9 rTimeconst <- .1 ## Variance of learning processes is sigma^2*t sigma <- c(Mechanics=.02,Fluency=.02) ## This function calculates what happens between two lessons ## The the Action arguments are assumed to vectors of the ## form (Mechanics, Fluency). The time argument should be a ## scalar. Success argument is optionally observed. ## The proficiency and success arguments should be matrixes ## whose rows correspond to particles and have the ## (Mechanics, Fluency) pattern. lessonEffect <- function (proficiency, action, time=1, success = runif(length(proficiency)) < actionSuccessP(proficiency,action) ){ if (!is.matrix(proficiency)) { proficiency <- matrix(proficiency,1) } if (!is.matrix(success)) { success <- matrix(success, nrow(proficiency), ncol(proficiency), byrow=TRUE) } rehearsal <- (rTimediff*success+rTimeconst)*time lambda <- learnRate(proficiency,action) sigmat <- sigma*sqrt(time) proficiency1 <- proficiency + lambda*rehearsal - forgetRate*time proficiency1 <- proficiency1 + matrix(rnorm(length(proficiency))*sigmat, nrow(proficiency),ncol(proficiency), byrow=TRUE) proficiency1 <- ifelse (proficiency1 < 0, 0, proficiency1) proficiency1 <- ifelse (proficiency1 > 6, 0, proficiency1) list(success=success,proficiency=proficiency1) } ## Mean and variance of the population. ## popA %*% t(popA) is the variance popA <- matrix(c(.6,.4,.4,.6),2,2) popMean <- c(Mechanics=1.5,Fluency=1.5) ## Generates the requested number of students ## returns a matrix with rows corresponding to students ## (particles) and columns corresponding to proficiencies. randomStudent <- function(nstudents=1) { nvar <- length(popMean) z <- matrix(rnorm(nstudents*nvar),nstudents,nvar) result <- sweep(z%*%popA,2,popMean,"+") #result <- as.vector(popA%*%rnorm(length(popMean))) + popMean result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) colnames(result)<-names(popMean) result } ## ncycles -- number of cycles to simulate ## estimator -- function used to estimate current student ## level from observables. ## policy -- function used to estimate next action from ## estimate ## evaluation -- function used to generate observables from ## current proficiencies ## advance -- function used to generate proficiencies at ## next time step from current time step. Note: assume ## that this returns a list of two values, success of action ## and new proficiencies ## initial -- function to generate random initial starting ## value for student. ## cycletime --- time between measurements (could be scaler ## or vector of length ncycles) simulation <- function(ncycles,estimator,policy, evaluation = evalLesson, advance = lessonEffect, initial = randomStudent, cycletime=1) { proficiency <- do.call(initial,list()) nvar <- length(proficiency) tnames <- paste("Time",0:ncycles) vnames <- colnames(proficiency) # We will fill these in as we go. truth <- matrix(NA,ncycles+1,nvar, dimnames=list(tnames,vnames)) obs <- truth est <- truth act <- truth suc <- truth if (length(cycletime) == 1) { cycletime <- rep(cycletime,ncycles) } ## Last cycle is stub only cycletime <- c(cycletime,0) # Now we start the simulation for (icycle in 1:ncycles) { truth[icycle,] <- proficiency observed <- do.call(evaluation,list(proficiency)) obs[icycle,] <- observed estimate <- do.call(estimator,list(obs,est,act,suc)) est[icycle,] <- estimate action <- do.call(policy,list(estimate)) act[icycle,] <- action result <- do.call(advance, list(proficiency, action, time=cycletime[icycle])) suc[icycle,] <- result$success proficiency <- result$proficiency } truth[ncycles+1,]<-proficiency observed <- do.call(evaluation,list(proficiency)) obs[ncycles+1,] <- observed sim <- data.frame(truth,obs,est,act,suc,cycletime) nsim <- paste("true",vnames,sep=".") nsim <- c(nsim,paste("obs",vnames,sep=".")) nsim <- c(nsim,paste("est",vnames,sep=".")) nsim <- c(nsim,paste("act",vnames,sep=".")) nsim <- c(nsim,paste("suc",vnames,sep=".")) nsim <- c(nsim,"cycletime") names(sim) <- nsim sim } ## Estimate current position by a running average of last ## three observables. Note that action and success history ## are passed so we can use particle filter type estimators. ## obs, est, act, and success are assumed to be matrixes ## with rows representing time points and columns ## representing variables. Values for not yet reached time ## points are assumed to be NA. raveEstimator <- function(obs,est,act,suc) { nobs <- sum(!is.na(obs[,1])) if (nobs == 1) return(obs[1,]) if (nobs < 3) { obs1 <- obs[!is.na(obs[,1]),] } else { obs1 <- obs[(nobs-2):nobs,] } result <- apply(obs1,2,mean) result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) result } expFilter <- function (obs,est,act,suc, alpha = filterAlpha, xi = meanInnovation) { nobs <- sum(!is.na(obs[,1])) if (nobs == 1) return(obs[1,]) result <- alpha*obs[nobs,] + (1-alpha)*(est[nobs-1,]+xi) result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) result } ## Generates an action from a proficiency estimate ## Simple policy which tries to match precisely with ## estimate floorPolicy <- function(estimate) { floor(2*estimate)/2 } ############################################### ## Revised Functions, taking pre-requisite effect into ## account. lessonEffectPreq <- function (proficiency,action,time=1, success = runif(length(proficiency)) < actionSuccessP(proficiency,action) ){ result <- lessonEffect(proficiency,action,time,success) prof <- result$proficiency # Prof 1 is pre-req for Prof 2 result$proficiency[,2] <- ifelse(prof[,1]>prof[,2], prof[,2],prof[,1]) result } randomStudentPreq <- function(nstudents=1) { nvar <- length(popMean) z <- matrix(rnorm(nstudents*nvar),nstudents,nvar) result <- sweep(z%*%popA,2,popMean,"+") #result <- as.vector(popA%*%rnorm(length(popMean))) + popMean result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) # Prof 1 is pre-req for Prof 2 result[,2] <- ifelse(result[,1]>result[,2], result[,2],result[,1]) colnames(result)<-names(popMean) result } ## Variation on our naive estimator which takes the pre-req ## condition into account raveEstimatorPreq <- function(obs,est,act,suc) { result <- raveEstimator(obs,est,act,suc) # Prof 1 is pre-req for Prof 2 result[2] <- ifelse(result[1]>result[2], result[2],result[1]) result } ## Variation on our naive estimator which takes the pre-req ## condition into account expFilterPreq <- function(obs,est,act,suc,alpha=filterAlpha, xi=meanInnovation) { result <- expFilter(obs,est,act,suc,alpha,xi) # Prof 1 is pre-req for Prof 2 result[2] <- ifelse(result[1]>result[2], result[2],result[1]) result } timeStep <- 1 varInnovation <- sigma*sigma*timeStep + (timeStep*rTimediff)^2*q0*(1-q0)*eta*eta filterAlpha <- evalSig^-2/(evalSig^-2+ 1/varInnovation) meanInnovation <- timeStep*rTimediff*q0*eta - timeStep*forgetRate