# Estimate Three-Parameter Normal Ogive Model. J.-P. Fox, University of Twente, Augustus 2010. ## K:: number of items ## N:: number persons ## DrawZ. Data Augmentation Z: Normal Ogive Model ## DrawS. Data Augmentation S: Knows the correct answes /guess the item correct ## DrawTheta. Draw values latent variable theta0 ## DrawBeta. Draw values difficulty parameter beta0 ## DrawAlpha. Draw values discrimination parameter alpha0 ## DrawC. Draw values guessing parameter guess0 DrawZ <- function(alpha0,beta0,theta0,S){ N <- nrow(S) K <- ncol(S) eta <- t(matrix(alpha0, ncol = N, nrow= K)) * matrix(theta0,ncol=K,nrow=N) - t(matrix(beta0, ncol = N, nrow = K)) BB <- matrix(pnorm(-eta),ncol = K, nrow = N) u <- matrix(runif(N*K), ncol = K, nrow = N) tt <- matrix( ( BB*(1-S) + (1-BB)*S )*u + BB*S, ncol = K, nrow = N) Z <- matrix(qnorm(tt),ncol = K, nrow = N) + matrix(eta,ncol = K, nrow = N) return(Z) } DrawBeta <- function(alpha0,theta0,Z) { # diffuse prior N <- nrow(Z) K <- ncol(Z) Zp <- t(matrix(alpha0,nrow= K,ncol = N)) * matrix(theta0,ncol=K,nrow=N) - matrix(Z,ncol=K,nrow=N) mu <- apply(Zp,2,sum)/N beta <- matrix(rnorm(K,mean=mu,sd=sqrt(1/N)),ncol=1,nrow=K) return(beta) } DrawAlpha <- function(beta0,theta0,Z){ # diffuse prior N <- nrow(Z) K <- ncol(Z) Zp <- (matrix(Z,ncol=K,nrow=N) + t(matrix(beta0,nrow=K,ncol=N))) * matrix(theta0,nrow=N,ncol=K) mu <- apply(Zp,2,sum)/sum(theta0^2) alpha <- matrix(rnorm(K,mean=mu,sd=sqrt(1/sum(theta0^2))),ncol=1,nrow=K) return(alpha) } DrawS <- function(alpha0,beta0,guess0,theta0,Y){ N <- nrow(Y) K <- ncol(Y) eta <- t(matrix(alpha0,ncol=N,nrow=K)) * matrix(theta0,ncol=K,nrow=N)-t(matrix(beta0,nrow=K,ncol=N)) eta <- matrix(pnorm(eta),ncol=K,nrow=N) probS <- eta/(eta + t(matrix(guess0,nrow=K,ncol=N))*(matrix(1,ncol=K,nrow=N)-eta)) S <- matrix(runif(N*K),ncol=K,nrow=N) S <- matrix(ifelse(S > probS,0,1),ncol=K,nrow=N) S <- S*Y return(S) } DrawC <- function(S,Y){ # Informative beta prior: c ~ Beta(P1,P2) P1 <- 2 P2 <- 6 N <- nrow(Y) K <- ncol(Y) Q1 <- P1 + apply((matrix(1,ncol=K,nrow=N)-S)*Y,2,sum) Q2 <- P2 + apply((matrix(1,ncol=K,nrow=N)-S),2,sum) - Q1 guess <- rbeta(K,Q1,Q2) return(guess) } THREEPNO <- function(Y,XG) { N <- nrow(Y) K <- ncol(Y) #Initialise alpha0 <- matrix(1,ncol = 1, nrow = 1) beta0 <- matrix(0,ncol=1,nrow=1) guess0 <- matrix(0,ncol=1,nrow=1) theta0 <- rnorm(N) #storage Malpha <- matrix(0,ncol=K,nrow=XG) Mbeta <- matrix(0,ncol=K,nrow=XG) Mguess <- matrix(0,ncol=K,nrow=XG) for(ii in 1:XG){ S <- DrawS(alpha0,beta0,guess0,theta0,Y) Z <- DrawZ(alpha0,beta0,theta0,S) beta0 <- DrawBeta(alpha0,theta0,Z) alpha0 <- DrawAlpha(beta0,theta0,Z) guess0 <- DrawC(S,Y) Malpha[ii,1:K] <- alpha0[1:K] Mbeta[ii,1:K] <- beta0[1:K] Mguess[ii,1:K] <- guess0[1:K] } return(list(Malpha,Mbeta,Mguess)) } ## simulate data N = 1000 K = 20 sim <- THREEPNO(N,K) Y <- matrix(sim[[1]],ncol=K,nrow=N) theta <- matrix(sim[[2]],ncol=1,nrow=N) alpha <- matrix(sim[[3]],ncol=1,nrow=K) beta <- matrix(sim[[4]],ncol=1,nrow=K) guess <- matrix(sim[[5]],ncol=1,nrow=K) out <- THREEPNO(Y,XG=5000) Malpha <- matrix(out[[1]],ncol=K) Mbeta <- matrix(out[[2]],ncol=K) Mguess <- matrix(out[[3]],ncol=K)