--- title: "MusicTutorExample" author: "Russell Almond" date: "10/20/2021" output: html_document runtime: shiny --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) library(shiny) library(tidyverse) library(plotly) library(DT) ``` ## Parameter Definitions ### Initial Skill Distribution ```{r pm, echo=FALSE} ## Mean and variance of the population. ## popA %*% t(popA) is the variance inputPanel( sliderInput("Mech", label = "Mechanics Mean", min=0, max=10, step=.05, value = 1.5), sliderInput("Flue",label = "Fluency Mean", min=0, max=10, step=.05, value = 1.5), sliderInput("s.mech", label= "Mechanics SD", min=0, max = 2, step=.05, value=.75), sliderInput("s.flue", label= "Fluency SD", min=0, max = 2, step=.05, value=.75), sliderInput("r.mf", label= "Mechanics--Fluency Correlation", min=-1, max = 1, step=.05, value=.9) ) popA <- reactive({diag(c(input$s.mech,input$s.flue))%*%chol(matrix(c(1,input$r.mf,input$r.mf,1),2,2))}) popMean <- reactive({ c(Mechanics=input$Mech,Fluency=input$Flue)}) ``` ```{r initialFunction} ## Generates the requested number of students ## returns a matrix with rows corresponding to students ## (particles) and columns corresponding to proficiencies. randomStudent <- function(nstudents=1,pMean=popMean(),pA=popA()) { nvar <- length(pMean) z <- matrix(rnorm(nstudents*nvar),nstudents,nvar) result <- sweep(z%*%pA,2,pMean,"+") #result <- as.vector(popA%*%rnorm(length(popMean))) + popMean result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) colnames(result)<-names(pMean) result } ``` ### Action Success Model ```{r q0, echo=FALSE} ## Constant noise term in noisy-or model #q0 <- c(Mechanics=.8,Fluency=.9) inputPanel( sliderInput("q0.Mech", label = "Mechanics Noise", min=0, max=1, step=.05, value = .8), sliderInput("q0.Flue",label = "Fluency Noise", min=0, max=1, step=.05, value = .9) ) q0 <- reactive(c(Mechanics=input$q0.Mech, Fluency=input$q0.Flue)) ``` ```{r highlow, echo=FALSE} ## Upper and lower bounds for delta in noisy-or model #threshlow <- c(Mechanics=-1,Fluency=-2) inputPanel( sliderInput("low.Mech", label = "Mechanics Delta Lower Bound", min=-6, max=0, step=.25, value = -1), sliderInput("low.Flue",label = "Fluency Delta Lower Bound", min=-6, max=0, step=.05, value = -2), #threshhigh <- c(Mechanics=2,Fluency=1) sliderInput("high.Mech", label = "Mechanics Delta Upper Bound", min=0, max=6, step=.25, value = 2), sliderInput("high.Flue",label = "Fluency Delta Upper Bound", min=0, max=6, step=.05, value = 1), ) ## Upper and lower bounds for delta in noisy-or model threshlow <- reactive(c(Mechanics=input$low.Mech, Fluency=input$low.Flue)) threshhigh <- reactive(c(Mechanics=input$high.Mech, Fluency=input$high.Mech)) ``` ```{r q1, echo=FALSE} ## skill bypass parameters in noisy-or model #q1 <- c(Mechanics=.8,Fluency=.7) inputPanel( sliderInput("q1.Mech", label = "Mechanics Bypass", min=0, max=1, step=.05, value = .8), sliderInput("q1.Flue",label = "Fluency Noise", min=0, max=1, step=.05, value = .7) ) q1 <- reactive(c(Mechanics=input$q1.Mech, Fluency=input$q1.Flue)) ## NOTE: These values are overridden for Simulation 2 (used ## in the paper). ``` ```{r action success} # 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, qq0=q0(), qq1=q1(), tlow=threshlow(), thigh=threshhigh()) { if (!is.matrix(proficiency)) { proficiency <- matrix(proficiency,1) } ## delta <- threshlow <= action-proficiency & ## action - proficiency <= threshhigh diff <- sweep(proficiency,2,action) delta <- sweep(diff,2,tlow,">=") & sweep(diff,2,thigh,"<=") ## q0 * prod(q1^(1-delta)) ## Need to transpose matrix to get rows varying fastest print(qq1^(1-t(delta))) qq1 <- apply(qq1^(1-t(delta)),2,prod) outer(qq1,qq0,"*") } ``` ## Growth Model ```{r Growth Model,echo=FALSE} inputPanel( sliderInput("eta.mech", label = "Learning Rate Mechanics", min=0, max=1, step=.01, value = .2), sliderInput("eta.flue",label = "Learning Rate Fluency", min=0, max=1, step=.01, value = .2), sliderInput("gamma.mech", label= "ZPD weight Mechanics", min=0, max = .5, step=.01, value=.1), sliderInput("gamma.flue", label= "ZPD weight Fluency", min=0, max = .5, step=.01, value=.15), sliderInput("forget.mech", label= "Forget Rate Mechanics", min=-.5, max = .5, step=.005, value=.02), sliderInput("forget.flue", label= "Forget Rate Fluency", min=-.5, max = .5, step=.005, value=.01) ) ## Constant term in learning rate eta <- reactive(c(Mechanics=input$eta.mech, Fluency=input$eta.flue)) ## Zone of Proximal Development term in learning rate gamma <- reactive(c(Mechanics=input$gamma.mech, Fluency=input$gamma.flue)) ## Constant parameters of the learning process forgetRate <- reactive(c(Mechanics=input$forget.mech, Fluency=input$forget.flue)) ``` ``` R(t) = (rTimediff*success + rTimeconst)*t ``` ```{r rates, echo=FALSE} ## R(t) = (rTimediff*success + rTimeconst)*t #rTimediff <- reactive(.9) #rTimeconst <- reactive(.1) inputPanel( sliderInput("rtd", label = "Time difference in learning rate", min=0, max=1, step=.01, value = .9), sliderInput("rtc",label = "Time constant learning rate", min=0, max=1, step=.01, value = .1) ) rTimediff <- reactive(input$rtd) rTimeconst <- reactive(input$rtc) ``` ## Variance of learning processes is sigma^2*t ```{r sigma, echo=FALSE} inputPanel( sliderInput("sig.Mech", label = "Learning Variance Mechanics", min=0, max=.25, step=.005, value = .02), sliderInput("sig.Flue",label = "Fluency Noise", min=0, max=.25, step=.005, value = .02) ) sigma <- reactive(c(Mechanics=input$sig.Mech, Fluency=input$sig.Flue)) ``` ```{r growth functions} ## 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, Eta=eta(), Gamma=gamma()) { 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,"+") } ## 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), rtdiff=rTimediff(), rtconst=rTimeconst(), sig=sigma(), fRate=forgetRate() ){ print("Do Lesson") if (!is.matrix(proficiency)) { proficiency <- matrix(proficiency,1) } if (!is.matrix(success)) { success <- matrix(success, nrow(success), ncol(success), byrow=TRUE) } print(success) rehearsal <- (rtdiff*success+rtconst)*time print(rehearsal) lambda <- learnRate(proficiency,action) print(lambda) sigmat <- sig*sqrt(time) proficiency1 <- proficiency + lambda*rehearsal - fRate*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) } 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 } ``` ### Evidence Model ```{r em} ## Regression coefficients for multivariate regression. evalA <- reactive(matrix(c(.7,.3,.3,.7),2,2)) ## Noise std for multivariate regression. evalSig <- reactive(c(.65,.65)) ``` ```{r em functions} ## This function takes the proficiency expressed as a vector ## (Mechanics, Fluency) and generates a random lesson value. evalLesson <- function(proficiency, eA=evalA(), eSig=evalSig()) { result <- proficiency %*% eA + rnorm(length(eSig))*eSig result <- round(2*result)/2 result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) } ``` ## Policies ## Estimation Functions ```{r estfun} ## 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. runlen <- reactive(3) raveEstimator <- function(obs,est,act,suc,runlen=runlen()) { nobs <- sum(!is.na(obs[,1])) if (nobs == 1) return(obs[1,]) if (nobs < runlen) { obs1 <- obs[!is.na(obs[,1]),] } else { obs1 <- obs[(nobs-runlen+1):nobs,] } result <- apply(obs1,2,mean) result <- ifelse((result < 0),0,result) result <- ifelse((result > 6),6,result) result } # This filter has two parameters filterAlpha <- reactive(.8) meanInnovation <- reactive (0) 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. 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 } ``` ## Simulation # The Simulation ## Simulation Parameters ```{r sim param, echo=FALSE} inputPanel( selectInput("N", label = "Number of Cycles:", choices = (1:10)*5, selected = 10), numericInput("seed", label = "Random number Seed (integer)", min = 0, max = .Machine$integer.max, value = floor(runif(1)*.Machine$integer.max), step = 157) ) ``` ```{r simulation} ## 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) print(tnames) 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 print(ncycles) 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) { print(icycle) truth[icycle,] <- proficiency print("Observed") observed <- do.call(evaluation,list(proficiency)) obs[icycle,] <- observed print("Estimate") estimate <- do.call(estimator,list(obs,est,act,suc)) est[icycle,] <- estimate print("Action") action <- do.call(policy,list(estimate)) act[icycle,] <- action print("Result") result <- do.call(advance, list(proficiency, action, time=cycletime[icycle])) str(result) suc[icycle,] <- result$success proficiency <- result$proficiency } print("Last Obs") truth[ncycles+1,]<-proficiency observed <- do.call(evaluation,list(proficiency)) obs[ncycles+1,] <- observed sim <- data.frame(cycle=0:ncycles,truth,obs,est,act,suc,cycletime) nsim <- c("cycle",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 } simdat <- reactive({input$seed;simulation(as.numeric(input$N), expFilter,floorPolicy)}) DT::renderDataTable(simdat()) ``` ```{r plot} cols <- c("true"="red","estimate"="purple","observed"="blue") value_type <- factor(c("true","estimate","observed")) shapes <- c("true"="asterisk", "estimate"="cross","observed"="triangle") renderPlotly({ plot <- ggplot(simdat(),aes(x=true.Mechanics,y=true.Fluency,frame=cycle)) + geom_point(aes(x=true.Mechanics,y=true.Fluency, shape=value_type[1], colour = value_type[1], size=5)) + geom_point(aes(x=obs.Mechanics,y=obs.Fluency, shape=value_type[3], colour = value_type[3], size=5)) + geom_point(aes(x=est.Mechanics,y=est.Fluency, shape=value_type[2], colour = value_type[2], size=5)) + scale_shape_manual(values=shapes) + scale_color_manual(values=cols) ggplotly(plot) }) ```