#  Estimate Two-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 
##	DrawTheta. Draw values latent variable theta
## DrawBeta. Draw values difficulty parameter beta0
## DrawAlpha. Draw values discrimination parameter alpha0

DrawZ 	<- function(alpha0,beta0,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, 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-Y) + (1-BB)*Y )*u + BB*Y, ncol = K, nrow = N)
	Z		<- matrix(qnorm(tt),ncol = K, nrow = N) + matrix(eta,ncol = K, nrow = N) 
  return(Z)
}

DrawTheta <- function(alpha0,beta0,Z) {
#prior theta N(0,1)
	N		<- nrow(Z)
	K 		<- ncol(Z) 
	pvar		<- (sum(alpha0^2)+1)
	thetahat	<- (((Z + t(matrix(beta0,ncol= N, nrow = K))) %*% matrix(alpha0,ncol = 1,nrow = K)))
	mu		<- thetahat/pvar
	theta		<- rnorm(N,mean=mu,sd=sqrt(1/pvar))
	theta		<- (theta - mean(theta))#/sqrt(var(theta))
return(theta)
}


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)
}
 

TWOPNO <- function(Y,XG)
{
	N 		<- nrow(Y)
	K 		<- ncol(Y) 

#Initialise
	alpha0 	<- matrix(1,ncol = 1, nrow = 1)
	beta0 		<- matrix(0,ncol=1,nrow=1)
	theta0 	<- matrix(rnorm(N),ncol=1,nrow=N)
	EAPtheta 	<- matrix(0,ncol=1,nrow=N)
	SQtheta 	<- matrix(0,ncol=1,nrow=N)
	Malpha 	<- matrix(0,ncol=K,nrow=XG)
	Mbeta 		<- matrix(0,ncol=K,nrow=XG)

for(ii in 1:XG){

	Z					<- DrawZ(alpha0,beta0,theta0,Y) 
	theta0				<- DrawTheta(alpha0,beta0,Z)
	beta0				<- DrawBeta(alpha0,theta0,Z) 
	alpha0 			<- DrawAlpha(beta0,theta0,Z)	
	EAPtheta 			<- ((ii - 1)*EAPtheta + theta0)/ii
	SQtheta			<- ((ii - 1)*SQtheta + theta0^2)/ii
	Malpha[ii,1:K]	<- alpha0[1:K]
	Mbeta[ii,1:K] 	<- beta0[1:K]
}	
	
	return(list(Malpha,Mbeta,EAPtheta,SQtheta))
}


N	<- 1000
K	<- 20
sim		<-	sim2PNO(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)

out <- TWOPNO(Y,XG=1000)
Malpha <- matrix(out[[1]],ncol=K)
Mbeta <- matrix(out[[2]],ncol =K)
EAPtheta <- out[[3]]
SDtheta <- sqrt(out[[4]] - out[[3]]^2) ### estimate standard deviation theta

plot(theta[1:N],EAPtheta[1:N])
abline(0,1)

### Use estMLIRT program
S <- c(0,0,0,0) 
out2 <- estMLIRT(Y, S, nll=N, XG=1000) 
MLIRTout(500,out2)