# R/plfm.R In plfm: Probabilistic Latent Feature Analysis

#### Documented in bayesplfmgendatgendatLCplfmLCplfmplfmplot.LCplfmplot.plfmplot.stepLCplfmplot.stepplfmprint.bayesplfmprint.LCplfmprint.plfmprint.stepLCplfmprint.stepplfmprint.summary.bayesplfmprint.summary.plfmprint.summary.stepLCplfmprint.summary.stepplfmstepLCplfmstepplfmsummary.bayesplfmsummary.plfmsummary.stepLCplfmsummary.stepplfm

##################################################################
### function gendat
##################################################################

gendat<-function(maprule="disj",N,objpar,attpar){
## function probability disjunctive model

pdisj<-function(F,objpar,attpar){
prob<-matrix(rep(1,J*K),nrow=J)
for (f in 1:F) {prob<-prob*(1-objpar[,f]%o%attpar[,f])}
pdisj<-1-prob
}

## function probability conjunctive model

pconj<-function(F,objpar,attpar){
prob<-matrix(rep(1,J*K),nrow=J)
for (f in 1:F) {prob<-prob*(1-(1-objpar[,f])%o%attpar[,f])}
pconj<-prob
}

prob<-matrix(rep(0,J*K),nrow=J)
for (f in 1:F) {prob<-prob+(objpar[,f]%o%attpar[,f])}
}

J<-dim(objpar)[1]
F<-dim(objpar)[2]
K<-dim(attpar)[1]

if ((dim(objpar)[2]!=dim(attpar)[2])) print("Number of columns in matrix of object parameters and matrix of attribute parameters should be the same")
else if (sum((objpar<0)|(objpar>1)|is.nan(objpar))>0) print("Object parameters should all be non-missing and between 0 and 1")
else if (sum((attpar<0)|(attpar>1)|is.nan(attpar))>0) print("Attribute parameters should all be between 0 and 1")
else if ((maprule!="disj")&(maprule!="conj")&(maprule!="add")) print("Incorrect specification of mapping rule")
else if (N<0) print("The number of replications should be positive (N>0)")
else if (N==0)
{
if (maprule=="disj") prob1<-pdisj(F,objpar,attpar)
else if (maprule=="conj") prob1<-pconj(F,objpar,attpar)
gendat<-list(call=match.call(),prob1=prob1)
}
else if (N>0){
if (maprule=="disj") prob1<-pdisj(F,objpar,attpar)
else if (maprule=="conj") prob1<-pconj(F,objpar,attpar)

freqtot<-matrix(rep(N,J*K),nrow=J)
freq1<-matrix(rbinom(J*K,N,prob1),nrow=J)
gendat<-list(call=match.call(),prob1=prob1,freq1=freq1,freqtot=freqtot)
}
}

##################################################################
### function PLFM
##################################################################

plfm<-function(data,object,attribute,rating,freq1,freqtot,F,datatype="freq",maprule="disj",M=5,emcrit1=1e-2,emcrit2=1e-10,printrun=TRUE){

#### functions

piDC<-function(J,K,F,sigma,rho){
prob<-matrix(rep(1,J*K),nrow=J)
for (f in 1:F) {prob<-prob*(1-sigma[,f]%o%rho[,f])}
piDC<-1-prob
}

prob<-matrix(rep(0,J*K),nrow=J)
for (f in 1:F) {prob<-prob+(sigma[,f]%o%rho[,f])}
}

logpost<-function(p1,sigma,rho,freq1,freq0){
logpost<-sum(log(sigma)+log(1-sigma))+sum(log(rho)+log(1-rho))+sum(freq1*log(p1)+freq0*log(1-p1))
}

loglik<-function(p1,freq1,freq0){
loglik<-sum(freq1*log(p1)+freq0*log(1-p1))
}

computeHessianDC<-function(sigma.o,rho.o,psi,nsigma.o,nrho.o,nfreq1,nfreq0,np1.o,freq1mat,pi1mat,J,K,F,NPsigma,NPtot){
## Hsigma

Hsigma<-1/sigma.o**2+1/((1-sigma.o)**2)+apply((nrho.o/(1-psi))**2*(nfreq1*((1-np1.o)/np1.o)**2+nfreq0),c(1,3),sum)

## Hrho

Hrho<-1/rho.o**2+1/((1-rho.o)**2)+apply((nsigma.o/(1-psi))**2*(nfreq1*((1-np1.o)/np1.o)**2+nfreq0),c(2,3),sum)

## Hsigma.common

temp1<-c(t(rho.o))%o%rep(1,F)
temp2<-matrix(c(aperm((rep(1,F)%o%rho.o),c(3,1,2))),nrow=F*K,byrow=TRUE)
temp3<-aperm(array(c(t(temp1*temp2)),c(F,F,K)),c(3,1,2))
temp4<-rep(1,J)%o%temp3
temp5<-c(aperm(psi,c(3,1,2)))%o%rep(1,F)
temp6<-matrix(c(aperm(psi%o%rep(1,F),c(3,4,1,2))),nrow=J*K*F,byrow=TRUE)
temp7<-(1-temp5)*(1-temp6)
temp8<-aperm(array(t(temp7),c(F,F,J,K)),c(3,4,1,2))
Hsigma.common<-apply(freq1mat*temp4*(1-pi1mat)/((pi1mat**2)*temp8),c(1,3,4),sum)

## Hrho.common

temp1<-c(t(sigma.o))%o%rep(1,F)
temp2<-matrix(c(aperm((rep(1,F)%o%sigma.o),c(3,1,2))),nrow=F*J,byrow=TRUE)
temp3<-aperm(array(c(t(temp1*temp2)),c(F,F,J)),c(3,1,2))
temp4<-rep(1,K)%o%temp3
temp4<-aperm(temp4,c(2,1,3,4))
Hrho.common<-apply(freq1mat*temp4*(1-pi1mat)/((pi1mat**2)*temp8),c(2,3,4),sum)

## Hsigmarho.common

Hsigmarho.common<-nfreq1*((1-np1.o)/(1-psi))*(-(1/np1.o)+(psi*(1-np1.o)/((1-psi)*np1.o**2)))+nfreq0/((1-psi)**2)

## Hsigmarho.distinct

temp9<-(rep(1,J)%o%rho.o)%o%rep(1,F)
temp10<-rep(1,K)%o%(aperm(array(t(c(t(sigma.o))%o%rep(1,F)),c(F,F,J)),c(3,1,2)))
temp11<-aperm(temp10,c(2,1,3,4))
temp12<-aperm(temp9*temp11,c(1,2,4,3))
temp13<-aperm(sigma.o%o%rep(1,K),c(1,3,2))
temp14<-rep(1,J)%o%rho.o
temp15<-(temp13*temp14)%o%(rep(1,F))
temp16<-aperm(temp15,c(1,2,4,3))
temp17<-(1-temp15)*(1-temp16)
Hsigmarho.distinct<-aperm((freq1mat*temp12*(1-pi1mat))/((pi1mat**2)*temp17),c(1,2,4,3))

#########################
##complete Hessian matrix
#########################

## deriv(sigma,sigma)

Hes.sigma<-array(rep(0,J*J*F*F),c(J,F,J,F))

for (j1 in 1:J){
for (f1 in 1:F){
Hes.sigma[j1,f1,j1,f1]<-Hsigma[j1,f1]}}

for (j in 1:J){
for (f1 in 1:F){
for (f2 in 1:F){

if (f1!=f2) Hes.sigma[j,f1,j,f2]<-Hsigma.common[j,f1,f2]}}}

Hesmat.sigma<-matrix(aperm(Hes.sigma,c(2,1,4,3)),nrow=J*F)

## deriv(rho,rho)

Hes.rho<-array(rep(0,K*K*F*F),c(K,F,K,F))

for (k in 1:K){
for (f1 in 1:F){
Hes.rho[k,f1,k,f1]<-Hrho[k,f1]}}

for (k in 1:K){
for (f1 in 1:F){
for (f2 in 1:F){
if (f1!=f2) Hes.rho[k,f1,k,f2]<-Hrho.common[k,f1,f2]}}}

Hesmat.rho<-matrix(aperm(Hes.rho,c(2,1,4,3)),nrow=K*F)

##deriv(sigma,rho)

Hes.sigmarho<-array(rep(0,J*F*K*F),c(J,F,K,F))

for (j in 1:J){
for (k in 1:K){
for (f1 in 1:F){
Hes.sigmarho[j,f1,k,f1]<-Hsigmarho.common[j,k,f1]}}}

for (j in 1:J){
for (k in 1:K){
for (f1 in 1:F){
for (f2 in 1:F){
if (f1!=f2) Hes.sigmarho[j,f1,k,f2]<-Hsigmarho.distinct[j,k,f1,f2]}}}}

Hesmat.sigmarho<-matrix(c(aperm(Hes.sigmarho,c(2,1,4,3))),nrow=J*F)

##deriv(rho,sigma)

Hessian<-matrix(rep(0,(NPtot)**2),nrow=NPtot)
Hessian[1:NPsigma,1:NPsigma]<-Hesmat.sigma
Hessian[1:NPsigma,(NPsigma+1):NPtot]<-Hesmat.sigmarho
Hessian[((NPsigma+1):NPtot),1:NPsigma]<-t(Hesmat.sigmarho)
Hessian[(NPsigma+1):NPtot,(NPsigma+1):NPtot]<-Hesmat.rho
computeHessianDC<-Hessian

}

## parameter estimation

if (datatype=="freq") {
if (length(freqtot)==1) freqtot<-matrix(rep(freqtot,((dim(freq1)[1])*(dim(freq1)[2]))),nrow=dim(freq1)[1])
}

if (datatype=="dataframe"){
agg<-table(data\$object,data\$attribute,data\$rating)
freq1<-agg[,,"1"]
freqtot<-agg[,,"0"]+agg[,,"1"]}

if (sum(is.nan(freq1)|is.na(freq1)|freq1<0)){print("The elements of the matrix freq1 should be non-missing and larger than or equal to 0")}
else if (sum(is.nan(freqtot)|is.na(freqtot)|freqtot<1)){print("The elements of the matrix freqtot should be non-missing and larger than 0")}
else if (sum(abs(dim(freqtot)-dim(freq1)))>0){print("The matrices freq1 and freqtot should have the same dimensionality")}
else if (sum(freq1>freqtot)>0) {print("The corresponding values of the matrix freqtot should be larger than or equal to the elements in the matrix freq1")}
else if ((maprule!="disj") & (maprule!="conj") & (maprule!="add")) {print("The mapping rule is not correctly specified")}
else if (M<1) {print("The number of requested runs should be at least one (M>=1)")}
else if (((dim(freq1)[1])*(dim(freq1)[2]))<F*sum(dim(freq1))) {print("The number of model parameters (J+K)*F should be smaller than the number of observations J*K used to fit the model")}
else if (emcrit1<0 | emcrit2<0) {print("emcrit1 and emcrit2 should be small positive numbers")}

else{

if (maprule=="disj"|maprule=="conj"){

## initialisation

N<-max(freqtot)
J<-dim(freq1)[1]
K<-dim(freq1)[2]
freq0<-freqtot-freq1

if (is.null(rownames(freq1))=="FALSE") rowlabels<-rownames(freq1)
else rowlabels<-paste("R_", seq(1,J),sep="")
if (is.null(colnames(freq1))=="FALSE") columnlabels<-colnames(freq1)
else columnlabels<-paste("C_", seq(1,K),sep="")

if (maprule=="conj"){
freq0<-freq1
freq1<-freqtot-freq0
}

NPsigma<-J*F
NPrho<-K*F
NPtot<-J*F+K*F

featurelabels<-paste("F",seq(1,F),sep="")
runlabels<-paste("RUN",seq(1,M),sep="")

nfreq1<-freq1%o%rep(1,F)
nfreq0<-freq0%o%rep(1,F)
nfreqtot<-freqtot%o%rep(1,F)
freq1mat<-(freq1%o%rep(1,F))%o%rep(1,F)

margobj<-apply(freqtot,1,sum)%o%rep(1,F)
margatt<-apply(freqtot,2,sum)%o%rep(1,F)

psi<-array(rep(0,J*K*F),c(J,K,F))
ksi<-array(rep(0,J*K*F),c(J,K,F))
eps<-array(rep(0,J*K*F),c(J,K,F))
gamma<-array(rep(0,J*K*F),c(J,K,F))
alfa<-array(rep(0,J*K*F),c(J,K,F))

sigma.runs<-array(rep(0,J*F*M),c(J,F,M))
rho.runs<-array(rep(0,K*F*M),c(K,F,M))
logpost.runs<-rep(0,M)

dimnames(sigma.runs)[[1]]<-rowlabels
dimnames(sigma.runs)[[2]]<-featurelabels
dimnames(sigma.runs)[[3]]<-runlabels
dimnames(rho.runs)[[1]]<-columnlabels
dimnames(rho.runs)[[2]]<-featurelabels
dimnames(rho.runs)[[3]]<-runlabels
names(logpost.runs)<-runlabels

for (run in 1:M){

## set convergence criterion to switch from EM to EM+NR
emcritX<-emcrit1

## initialisation starting values

sigma.o<-matrix(runif(J*F),nrow=J,byrow=TRUE)
sigma.n<-matrix(runif(J*F),nrow=J,byrow=TRUE)
rho.o<-matrix(runif(K*F),nrow=K,byrow=TRUE)
rho.n<-matrix(runif(K*F),nrow=K,byrow=TRUE)
theta.o<-c(t(sigma.o),t(rho.o))
theta.n<-c(t(sigma.n),t(rho.n))
p1.o<-piDC(J,K,F,sigma.o,rho.o)
p1.n<-piDC(J,K,F,sigma.n,rho.n)
logpost.o<-logpost(p1.o,sigma.o,rho.o,freq1,freq0)
logpost.n<-logpost(p1.n,sigma.n,rho.n,freq1,freq0)

while ((abs(logpost.o-logpost.n))>emcrit2){

flagNR<-0
sigma.o<-sigma.n
rho.o<-rho.n
nrho.o<-rep(1,J)%o%rho.o
nsigma.o<-aperm(sigma.o%o%rep(1,K),c(1,3,2))
theta.o<-c(t(sigma.o),t(rho.o))

p1.o<-piDC(J,K,F,sigma.o,rho.o)
np1.o<-p1.o%o%rep(1,F)
pi1mat<-np1.o%o%rep(1,F)

logpost.o<-logpost(p1.o,sigma.o,rho.o,freq1,freq0)

for (f in 1:F){psi[,,f]<-(sigma.o[,f]%o%rho.o[,f])}
for (f in 1:F){ksi[,,f]<-(sigma.o[,f]%o%(1-rho.o[,f]))}
for (f in 1:F){eps[,,f]<-((1-sigma.o[,f])%o%rho.o[,f])}

gamma<-ksi/(1-psi)
alfa<-eps/(1-psi)
sigma.n<-(1+apply(nfreq1*(1-gamma)*psi*(1/np1.o)+nfreqtot*gamma,c(1,3),sum))/(2+margobj)
rho.n<-(1+apply(nfreq1*(1-alfa)*psi*(1/np1.o)+nfreqtot*alfa,c(2,3),sum))/(2+margatt)
p1.n<-piDC(J,K,F,sigma.n,rho.n)
logpost.n<-logpost(p1.n,sigma.n,rho.n,freq1,freq0)

while ((flagNR==0) & (abs(logpost.o-logpost.n)<emcritX) & (abs(logpost.o-logpost.n)>emcrit2)){

##compute EM step

sigma.o<-sigma.n
rho.o<-rho.n
nrho.o<-rep(1,J)%o%rho.o
nsigma.o<-aperm(sigma.o%o%rep(1,K),c(1,3,2))
theta.o<-c(t(sigma.o),t(rho.o))

p1.o<-piDC(J,K,F,sigma.o,rho.o)
np1.o<-p1.o%o%rep(1,F)
pi1mat<-np1.o%o%rep(1,F)

logpost.o<-logpost(p1.o,sigma.o,rho.o,freq1,freq0)

for (f in 1:F){psi[,,f]<-(sigma.o[,f]%o%rho.o[,f])}
for (f in 1:F){ksi[,,f]<-(sigma.o[,f]%o%(1-rho.o[,f]))}
for (f in 1:F){eps[,,f]<-((1-sigma.o[,f])%o%rho.o[,f])}

gamma<-ksi/(1-psi)
alfa<-eps/(1-psi)

sigma.EM<-(1+apply(nfreq1*(1-gamma)*psi*(1/np1.o)+nfreqtot*gamma,c(1,3),sum))/(2+margobj)
rho.EM<-(1+apply(nfreq1*(1-alfa)*psi*(1/np1.o)+nfreqtot*alfa,c(2,3),sum))/(2+margatt)
theta.EM<-c(t(sigma.EM),t(rho.EM))

##compute minus second derivative augmented posterior

ES1<-apply(nfreq1*(1/np1.o)*(psi+ksi*(np1.o-psi)/(1-psi))+nfreq0*ksi/(1-psi),c(1,3),sum)
ET1<-apply(nfreq1*(1/np1.o)*(psi+eps*(np1.o-psi)/(1-psi))+nfreq0*eps/(1-psi),c(2,3),sum)
secderivsigma<-(1+ES1)/(sigma.o**2)+(1+margobj-ES1)/((1-sigma.o)**2)
secderivrho<-(1+ET1)/(rho.o**2)+(1+margatt-ET1)/((1-rho.o)**2)
secderiv<-diag(J*F+K*F)*c(t(secderivsigma),t(secderivrho))

##compute Hessian

completeHessian<-computeHessianDC(sigma.o,rho.o,psi,nsigma.o,nrho.o,nfreq1,nfreq0,np1.o,freq1mat,pi1mat,J,K,F,NPsigma,NPtot)

##compute acceleration step

theta.n<-theta.o+solve(completeHessian)%*%(secderiv%*%(theta.EM-theta.o))
sigma.n<-matrix(theta.n[1:NPsigma],nrow=J,byrow=TRUE)
rho.n<-matrix(theta.n[(NPsigma+1):NPtot],nrow=K,byrow=TRUE)

## check convergence Newton-Rhapson step
if (sum(ifelse(((sigma.n>0) & (sigma.n<1)),0,1))+sum(ifelse(((rho.n>0) & (rho.n<1)),0,1))>0)
{
sigma.n<-sigma.EM
rho.n<-rho.EM

emcritX<-emcritX/10
flagNR<-1
##print(c("change EMcrit",emcritX,emcrit1))
}

p1.n<-piDC(J,K,F,sigma.n,rho.n)
logpost.n<-logpost(p1.n,sigma.n,rho.n,freq1,freq0)

}## end loop EM+NR

}##end loop EM

sigma.runs[,,run]<-sigma.n
rho.runs[,,run]<-rho.n
logpost.runs[run]<-logpost.n

if (printrun=="TRUE"){
if (maprule=="disj") {print(paste("DISJUNCTIVE ANALYSIS  ","F=",F,"  RUN=",run,sep=""),quote="FALSE")}
else if (maprule=="conj") {print(paste("CONJUNCTIVE ANALYSIS  ","F=",F,"  RUN=",run,sep=""),quote="FALSE")}
}

}## end M

### select best solution

## select best solution and compute fit measures
best<-order(logpost.runs,decreasing=TRUE)[1]
sigma.best<-as.matrix(sigma.runs[,,best])
rownames(sigma.best)<-rowlabels
colnames(sigma.best)<-featurelabels
rho.best<-as.matrix(rho.runs[,,best])
rownames(rho.best)<-columnlabels
colnames(rho.best)<-featurelabels
p1.best<-piDC(J,K,F,sigma.best,rho.best)
rownames(p1.best)<-rowlabels
colnames(p1.best)<-columnlabels
np1.best<-p1.best%o%rep(1,F)
pi1mat.best<-np1.best%o%rep(1,F)
for (f in 1:F){psi[,,f]<-(sigma.best[,f]%o%rho.best[,f])}
nrho.best<-rep(1,J)%o%rho.best
nsigma.best<-aperm(sigma.best%o%rep(1,K),c(1,3,2))

##compute Hessian best solution

completeHessian<-computeHessianDC(sigma.best,rho.best,psi,nsigma.best,nrho.best,nfreq1,nfreq0,np1.best,freq1mat,pi1mat.best,J,K,F,NPsigma,NPtot)
SEsigma<-1/sqrt(matrix(diag(completeHessian)[1:NPsigma],nrow=J,byrow=TRUE))
rownames(SEsigma)<-rowlabels
colnames(SEsigma)<-featurelabels

SErho<-1/sqrt(matrix(diag(completeHessian)[(NPsigma+1):NPtot],nrow=K,byrow=TRUE))
rownames(SErho)<-columnlabels
colnames(SErho)<-featurelabels

## compute fit measures

loglik.best<-loglik(p1.best,freq1,freq0)
logpost.best<-logpost(p1.best,sigma.best,rho.best,freq1,freq0)
correlation<-cor(c(p1.best*freqtot),c(freq1))
VAF<-correlation**2
deviance<--2*loglik.best
AIC<-deviance+2*NPtot
BIC<-deviance+NPtot*log(N)
chisquare.fit<-sum(((freq1-freqtot*p1.best)**2/(freqtot*p1.best*(1-p1.best))))
df<-(J*K)-(J+K)*F
pchi<-1-pchisq(chisquare.fit,df)
fitmeasures<-matrix(c(loglik.best,logpost.best,deviance,AIC,BIC,chisquare.fit,df,pchi,correlation,VAF),ncol=1)
rownames(fitmeasures)<-c("Log Likelihood","Log Posterior","Deviance","AIC","BIC","Pearson Chi-square","df","p-value","Correlation observed and expected frequencies","VAF observed frequencies")

}#end disjunctive and conjunctive rule

## initialisation

N<-max(freqtot)
J<-dim(freq1)[1]
K<-dim(freq1)[2]
freq0<-freqtot-freq1

if (is.null(rownames(freq1))=="FALSE") rowlabels<-rownames(freq1) else rowlabels<-paste("R_", seq(1,J),sep="")
if (is.null(colnames(freq1))=="FALSE") columnlabels<-colnames(freq1) else columnlabels<-paste("C_", seq(1,K),sep="")

if (maprule=="conj"){
freq0<-freq1
freq1<-freqtot-freq0
}

NPsigma<-J*F
NPrho<-K*F
NPtot<-J*F+K*F

featurelabels<-paste("F",seq(1,F),sep="")
runlabels<-paste("RUN",seq(1,M),sep="")

nfreq1<-freq1%o%rep(1,F)
nfreq0<-freq0%o%rep(1,F)
nfreqtot<-freqtot%o%rep(1,F)
freq1mat<-(freq1%o%rep(1,F))%o%rep(1,F)

margobj<-apply(freqtot,1,sum)%o%rep(1,F)
margatt<-apply(freqtot,2,sum)%o%rep(1,F)

psi<-array(rep(0,J*K*F),c(J,K,F))
ksi<-array(rep(0,J*K*F),c(J,K,F))
eps<-array(rep(0,J*K*F),c(J,K,F))

sigma.runs<-array(rep(0,J*F*M),c(J,F,M))
rho.runs<-array(rep(0,K*F*M),c(K,F,M))
logpost.runs<-rep(0,M)

dimnames(sigma.runs)[[1]]<-rowlabels
dimnames(sigma.runs)[[2]]<-featurelabels
dimnames(sigma.runs)[[3]]<-runlabels
dimnames(rho.runs)[[1]]<-columnlabels
dimnames(rho.runs)[[2]]<-featurelabels
dimnames(rho.runs)[[3]]<-runlabels
names(logpost.runs)<-runlabels

for (run in 1:M){

## initialisation starting values

sigma.o<-matrix(runif(J*F),nrow=J,byrow=TRUE)
sigma.n<-matrix(runif(J*F),nrow=J,byrow=TRUE)
rho.o<-matrix(runif(K*F),nrow=K,byrow=TRUE)
rho.n<-matrix(runif(K*F),nrow=K,byrow=TRUE)
logpost.o<-logpost(p1.o,sigma.o,rho.o,freq1,freq0)
logpost.n<-logpost(p1.n,sigma.n,rho.n,freq1,freq0)

while ((abs(logpost.o-logpost.n))>emcrit2){

sigma.o<-sigma.n
rho.o<-rho.n
p1.o<-p1.n
logpost.o<-logpost.n

psi<-array(rep(0,J*K*F),c(J,K,F))
ksi<-array(rep(0,J*K*F),c(J,K,F))
eps<-array(rep(0,J*K*F),c(J,K,F))

for (f in 1:F){psi[,,f]<-sigma.o[,f]%o%rho.o[,f]}
for (f in 1:F){ksi[,,f]<-sigma.o[,f]%o%(1-rho.o[,f])}
for (f in 1:F){eps[,,f]<-(1-sigma.o[,f])%o%(rho.o[,f])}

np1.o<-p1.o%o%rep(1,F)

px0<-array(rep(0,J*K*F),c(J,K,F))
px1<-array(rep(0,J*K*F),c(J,K,F))
py0<-array(rep(0,J*K*F),c(J,K,F))
py1<-array(rep(0,J*K*F),c(J,K,F))

px0<-((F-(1+F*np1.o-psi))*psi+(F-(F*np1.o-psi))*ksi)/(F*(1-np1.o))
px1<-((1+F*np1.o-psi)*psi+(F*np1.o-psi)*ksi)/(F*np1.o)
py0<-((F-(1+F*np1.o-psi))*psi+(F-(F*np1.o-psi))*eps)/(F*(1-np1.o))
py1<-((1+F*np1.o-psi)*psi+(F*np1.o-psi)*eps)/(F*np1.o)

sigma.n<-(1+apply(nfreq1*px1+nfreq0*px0,c(1,3),sum))/(margobj+2)
rho.n<-(1+apply(nfreq1*py1+nfreq0*py0,c(2,3),sum))/(margatt+2)

logpost.n<-logpost(p1.n,sigma.n,rho.n,freq1,freq0)

sigma.runs[,,run]<-sigma.n
rho.runs[,,run]<-rho.n
logpost.runs[run]<-logpost.n
}

}## end M

### select best solution

## select best solution and compute fit measures
best<-order(logpost.runs,decreasing=TRUE)[1]
sigma.best<-as.matrix(sigma.runs[,,best])
rownames(sigma.best)<-rowlabels
colnames(sigma.best)<-featurelabels
rho.best<-as.matrix(rho.runs[,,best])
rownames(rho.best)<-columnlabels
colnames(rho.best)<-featurelabels
rownames(p1.best)<-rowlabels
colnames(p1.best)<-columnlabels
np1.best<-p1.best%o%rep(1,F)
nrho.best<-rep(1,J)%o%rho.best
nsigma.best<-aperm(sigma.best%o%rep(1,K),c(1,3,2))

##compute standard error

SEsigma<-1/(sigma.best^2)+1/((1-sigma.best)^2)+ apply(nfreq1*(nrho.best/(F*np1.best))^2+nfreq0*(nrho.best/(F*(1-np1.best)))^2,c(1,3),sum)
SEsigma<-1/sqrt(SEsigma)
rownames(SEsigma)<-rowlabels
colnames(SEsigma)<-featurelabels

SErho<-1/(rho.best^2)+1/((1-rho.best)^2)+ apply(nfreq1*(nsigma.best/(F*np1.best))^2+nfreq0*(nsigma.best/(F*(1-np1.best)))^2,c(2,3),sum)
SErho<-1/sqrt(SErho)
rownames(SErho)<-columnlabels
colnames(SErho)<-featurelabels

## compute fit measures

loglik.best<-loglik(p1.best,freq1,freq0)
logpost.best<-logpost(p1.best,sigma.best,rho.best,freq1,freq0)
correlation<-cor(c(p1.best*freqtot),c(freq1))
VAF<-correlation**2
deviance<--2*loglik.best
AIC<-deviance+2*NPtot
BIC<-deviance+NPtot*log(N)
chisquare.fit<-sum(((freq1-freqtot*p1.best)**2/(freqtot*p1.best*(1-p1.best))))
df<-(J*K)-(J+K)*F
pchi<-1-pchisq(chisquare.fit,df)
fitmeasures<-matrix(c(loglik.best,logpost.best,deviance,AIC,BIC,chisquare.fit,df,pchi,correlation,VAF),ncol=1)
rownames(fitmeasures)<-c("Log Likelihood","Log Posterior","Deviance","AIC","BIC","Pearson Chi-square","df","p-value","Correlation observed and expected frequencies","VAF observed frequencies")

else if (maprule=="conj")

class(plfm)<-"plfm"
plfm
}
}

############################
## function print.plfm
############################

print.plfm<-function(x, ...)
{
cat("Call:\n")
print(x\$call)

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
descript<-as.matrix(x\$fitmeasures[9:10,])
colnames(descript)<-""
print(descript,digits=3)

cat("\nESTIMATE OBJECT PARAMETERS:\n")
strout<-as.matrix(capture.output(round(x\$objpar,2)))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\nESTIMATE ATTRIBUTE PARAMETERS:\n")
strout<-as.matrix(capture.output(round(x\$attpar,2)))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")
}

#############################################################
## function plot.plfm
############################################################

plot.plfm<-function(x,feature=1,element="object",cexsymb=1,cexlabel=1,...){

if (element=="object"){
J<-dim(x\$objpar)[1]
F<-dim(x\$objpar)[2]
rowlab<-rownames(x\$objpar)
upper<-x\$objpar+1.96*x\$SE.objpar
upper<-ifelse(upper>1,1,upper)
lower<-x\$objpar-1.96*x\$SE.objpar
lower<-ifelse(lower<0,0,lower)
##par(pty="s")
plot(-1,-1,xlim=c(-1,1),ylim=c(1,J),xlab="",ylab="",yaxt="n",bty="n",xaxp=c(0,1,4),...)
for (i in 1:J) {points(x\$objpar[i,feature],i,cex=cexsymb)}
for (i in 1:J) {lines(c(lower[i,feature],upper[i,feature]),c(i,i))}
for (i in 1:J) {text(-0.8,i,rowlab[i],cex=cexlabel)}
} ## end class-independent object parameters

else if (element=="attribute"){
K<-dim(x\$attpar)[1]
F<-dim(x\$attpar)[2]
rowlab<-rownames(x\$attpar)
upper<-x\$attpar+1.96*x\$SE.attpar
upper<-ifelse(upper>1,1,upper)
lower<-x\$attpar-1.96*x\$SE.attpar
lower<-ifelse(lower<0,0,lower)
##par(pty="s")
plot(-1,-1,xlim=c(-1,1),ylim=c(1,K),xlab="",ylab="",yaxt="n",bty="n",xaxp=c(0,1,4),...)
for (i in 1:K) {points(x\$attpar[i,feature],i,cex=cexsymb)}
for (i in 1:K) {lines(c(lower[i,feature],upper[i,feature]),c(i,i))}
for (i in 1:K) {text(-0.8,i,rowlab[i],cex=cexlabel)}
} ## end class-independent attribute parameters

}

###############################
## function summary.plfm
###############################

summary.plfm <- function(object, ...)
{

crit<-as.matrix(object\$fitmeasures[1:5,])
colnames(crit)<-""
crit<-round(crit)

chi<-as.matrix(object\$fitmeasures[6:8,])
colnames(chi)<-""
chi<-round(chi,3)

descript<-as.matrix(object\$fitmeasures[9:10,])
colnames(descript)<-""
descript<-round(descript,3)

featurelabels<-colnames(object\$objpar)
F<-length(featurelabels)
estobjectpar<-round(object\$objpar,2)
colnames(estobjectpar)<-paste("Est(",featurelabels,")",sep="")
seobjectpar<-round(object\$SE.objpar,3)
colnames(seobjectpar)<-paste("SE(",featurelabels,")",sep="")

estattpar<-round(object\$attpar,2)
colnames(estattpar)<-paste("Est(",featurelabels,")",sep="")
seattpar<-round(object\$SE.attpar,3)
colnames(seattpar)<-paste("SE(",featurelabels,")",sep="")

result<-list(call=object\$call,
informationcriteria=crit,
chisquaretest=chi,
descriptivefit=descript,
objpar=estobjectpar,
SE.objpar=seobjectpar,
attpar=estattpar,
SE.attpar=seattpar)
class(result)<-"summary.plfm"
result
}

##############################################
## function print.summary.plfm
##############################################

print.summary.plfm<-function(x, ...)
{
cat("Call:\n")
print(x\$call)

cat("\nINFORMATION CRITERIA:\n")
print(x\$informationcriteria)

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
print(x\$chisquaretest)

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
print(x\$descriptivefit)

cat("\nESTIMATE OBJECT PARAMETERS:\n")
strout<-as.matrix(capture.output(x\$objpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\nSTANDARD ERROR OBJECT PARAMETERS:\n")
strout<-as.matrix(capture.output(x\$SE.objpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$SE.objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\nESTIMATE ATTRIBUTE PARAMETERS:\n")
strout<-as.matrix(capture.output(x\$attpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\nSTANDARD ERROR ATTRIBUTE PARAMETERS:\n")
strout<-as.matrix(capture.output(x\$SE.attpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$SE.attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")
}

##############################
## function stepplfm
##############################

stepplfm<-function(minF,maxF,data,object,attribute,rating,freq1,freqtot,datatype="freq",maprule="disj",M=5,emcrit1=1e-2,emcrit2=1e-10,printrun=TRUE){

tempcall<-match.call()
tempcall<-tempcall[c(-2,-3)]
if (maprule=="disj") tempcall\$maprule<-"disj"
else if (maprule=="conj") tempcall\$maprule<-"conj"

if (maprule=="disj"){
stepplfm<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
stepplfm[[f-minF+1]]<-plfm(maprule="disj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepplfm[[f-minF+1]]\$call<-tempcall
}
}

else if (maprule=="conj"){
stepplfm<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
stepplfm[[f-minF+1]]<-plfm(maprule="conj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepplfm[[f-minF+1]]\$call<-tempcall

}
}

stepplfm<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
tempcall\$F<-as.numeric(f)
stepplfm[[f-minF+1]]\$call<-tempcall
}
}

else if (maprule=="disj/conj"){

stepdisj<-vector(mode="list",maxF-minF+1)
stepconj<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
stepdisj[[f-minF+1]]<-plfm(maprule="disj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepdisj[[f-minF+1]]\$call<-tempcall
}
for (f in (minF:maxF)){
stepconj[[f-minF+1]]<-plfm(maprule="conj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepconj[[f-minF+1]]\$call<-tempcall
}
stepplfm<-list(disj=stepdisj,conj=stepconj)
}

stepdisj<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
stepdisj[[f-minF+1]]<-plfm(maprule="disj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepdisj[[f-minF+1]]\$call<-tempcall
}
for (f in (minF:maxF)){
tempcall\$F<-as.numeric(f)
}
}

stepconj<-vector(mode="list",maxF-minF+1)
for (f in (minF:maxF)){
tempcall\$F<-as.numeric(f)
}
for (f in (minF:maxF)){
stepconj[[f-minF+1]]<-plfm(maprule="conj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepconj[[f-minF+1]]\$call<-tempcall
}
}

stepdisj<-vector(mode="list",maxF-minF+1)
stepconj<-vector(mode="list",maxF-minF+1)

for (f in (minF:maxF)){
stepdisj[[f-minF+1]]<-plfm(maprule="disj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepdisj[[f-minF+1]]\$call<-tempcall
}
for (f in (minF:maxF)){
stepconj[[f-minF+1]]<-plfm(maprule="conj",data=data,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=f,M=M,emcrit1=emcrit1,emcrit2=emcrit2,printrun=printrun)
tempcall\$F<-as.numeric(f)
stepconj[[f-minF+1]]\$call<-tempcall
}
for (f in (minF:maxF)){
tempcall\$F<-as.numeric(f)
}
}

class(stepplfm)<-"stepplfm"
stepplfm
}

############################
## function print.stepplfm
############################

print.stepplfm<-function(x, ...)
{
nfeat<-length(x)
minfeat<-dim(x[[1]]\$objpar)[2]
maxfeat<-dim(x[[nfeat]]\$objpar)[2]
fitm<-t(sapply(x,function(obj) obj\$fitmeasures))
colnames(fitm)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitm)<-paste("F=",seq(minfeat,maxfeat),sep="")

if (x[[1]]\$call\$maprule=="disj"){
cat("\n*****************************")
cat("\nDISJUNCTIVE MODEL\n")
cat("******************************\n")
}
else if (x[[1]]\$call\$maprule=="conj"){
cat("\n*********************************")
cat("\nCONJUNCTIVE MODEL\n")
cat("*********************************\n")
}

cat("\n*********************************")
cat("*********************************\n")
}

cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitm[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitm[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitm[,c(9,10)],3))
}

nfeat<-length(x\$disj)
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")

cat("\n*****************************")
cat("\nDISJUNCTIVE MODEL\n")
cat("******************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitdisj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(9,10)],3))

cat("\n*********************************")
cat("\nCONJUNCTIVE MODEL\n")
cat("*********************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitconj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(9,10)],3))
}

nfeat<-length(x\$disj)
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")

cat("\n*****************************")
cat("\nDISJUNCTIVE MODEL\n")
cat("******************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitdisj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(9,10)],3))

cat("\n*********************************")
cat("*********************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
}

fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")

cat("\n*********************************")
cat("\nCONJUNCTIVE MODEL\n")
cat("*********************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitconj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(9,10)],3))

cat("\n*****************************")
cat("******************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
}

nfeat<-length(x\$disj)
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")

cat("\n*****************************")
cat("\nDISJUNCTIVE MODEL\n")
cat("******************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitdisj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitdisj[,c(9,10)],3))

cat("\n*********************************")
cat("\nCONJUNCTIVE MODEL\n")
cat("*********************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")
print(round(fitconj[,c(1:5)],0))

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(6:8)],3))

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
print(round(fitconj[,c(9,10)],3))

cat("\n*****************************")
cat("******************************\n")
cat("\nINFORMATION CRITERIA:\n")
cat("\n")

cat("\nPEARSON CHI SQUARE TEST OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
cat("\n")
}

}

################################
## function plot.stepplfm
################################

plot.stepplfm<-function(x,which="BIC",...){

nfeat<-length(x)
if (nfeat==1){print("The fit of the model is only displayed if F>1",quote="FALSE")}
else{
minfeat<-dim(x[[1]]\$objpar)[2]
maxfeat<-dim(x[[nfeat]]\$objpar)[2]
fitm<-t(sapply(x,function(obj) obj\$fitmeasures))
colnames(fitm)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitm)<-paste("F=",seq(minfeat,maxfeat),sep="")
plot(seq(minfeat,maxfeat),fitm[,which],xlab="Number of features",ylab=which,type="b",xaxp=c(minfeat,maxfeat,maxfeat-minfeat),...)

}
}

nfeat<-length(x\$disj)
if (nfeat==1){print("The fit of the model is only displayed if F>1",quote="FALSE")}
else{
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")
comb<-rbind(fitdisj,fitconj)
minstat<-apply(comb,2,min)
maxstat<-apply(comb,2,max)

plot(seq(minfeat,maxfeat),fitdisj[,which],xlab="Number of features",ylab=which,type="b",xaxp=c(minfeat,maxfeat,maxfeat-minfeat),ylim=c(minstat[which],maxstat[which]),...)
lines(seq(minfeat,maxfeat),fitconj[,which],lty=2,type="b")
if ((which=="BIC")|(which=="AIC")|(which=="Deviance")|(which=="Chisquare")){
legend("topright",c("Disjunctive","Conjunctive"),lty=c(1,2),border=NULL,bty="n")}
else if ((which=="Correlation")|(which=="VAF")){
legend("bottomright",c("Disjunctive","Conjunctive"),lty=c(1,2),border=NULL,bty="n")}

}
}

nfeat<-length(x\$disj)
if (nfeat==1){print("The fit of the model is only displayed if F>1",quote="FALSE")}
else{
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")
minstat<-apply(comb,2,min)
maxstat<-apply(comb,2,max)

plot(seq(minfeat,maxfeat),fitdisj[,which],xlab="Number of features",ylab=which,type="b",xaxp=c(minfeat,maxfeat,maxfeat-minfeat),ylim=c(minstat[which],maxstat[which]),...)
if ((which=="BIC")|(which=="AIC")|(which=="Deviance")|(which=="Chisquare")){
else if ((which=="Correlation")|(which=="VAF")){

}
}

if (nfeat==1){print("The fit of the model is only displayed if F>1",quote="FALSE")}
else{
fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")
minstat<-apply(comb,2,min)
maxstat<-apply(comb,2,max)

lines(seq(minfeat,maxfeat),fitconj[,which],lty=2,type="b")
if ((which=="BIC")|(which=="AIC")|(which=="Deviance")|(which=="Chisquare")){
else if ((which=="Correlation")|(which=="VAF")){

}
}

nfeat<-length(x\$disj)
if (nfeat==1){print("The fit of the model is only displayed if F>1",quote="FALSE")}
else{
minfeat<-dim(x\$disj[[1]]\$objpar)[2]
maxfeat<-dim(x\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(x\$disj,function(obj) obj\$fitmeasures))
fitconj<-t(sapply(x\$conj,function(obj) obj\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("F=",seq(minfeat,maxfeat),sep="")
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("F=",seq(minfeat,maxfeat),sep="")
minstat<-apply(comb,2,min)
maxstat<-apply(comb,2,max)

plot(seq(minfeat,maxfeat),fitdisj[,which],xlab="Number of features",ylab=which,type="b",xaxp=c(minfeat,maxfeat,maxfeat-minfeat),ylim=c(minstat[which],maxstat[which]),...)
lines(seq(minfeat,maxfeat),fitconj[,which],lty=2,type="b")
if ((which=="BIC")|(which=="AIC")|(which=="Deviance")|(which=="Chisquare")){
else if ((which=="Correlation")|(which=="VAF")){

}
}

}

############################
## function summary.stepplfm
############################

summary.stepplfm<-function(object, ...)
{
nfeat<-length(object)
minfeat<-dim(object[[1]]\$objpar)[2]
maxfeat<-dim(object[[nfeat]]\$objpar)[2]
fitm<-t(sapply(object,function(x) x\$fitmeasures))
colnames(fitm)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")

if (object[[1]]\$call\$maprule=="disj"){
rownames(fitm)<-paste("DISJUNCTIVE F=",seq(minfeat,maxfeat),sep="")
}
else if (object[[1]]\$call\$maprule=="conj"){
rownames(fitm)<-paste("CONJUNCTIVE F=",seq(minfeat,maxfeat),sep="")
}
}

result<-fitm
class(result)<-"summary.stepplfm"
result
}

nfeat<-length(object\$disj)
minfeat<-dim(object\$disj[[1]]\$objpar)[2]
maxfeat<-dim(object\$disj[[nfeat]]\$objpar)[2]
fitdisj<-t(sapply(object\$disj,function(x) x\$fitmeasures))
fitconj<-t(sapply(object\$conj,function(x) x\$fitmeasures))
colnames(fitdisj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitdisj)<-paste("DISJUNCTIVE F=",seq(minfeat,maxfeat),sep="")
colnames(fitconj)<-c("LogLik","LogPost","Deviance","AIC","BIC","Chisquare","df","p-value","Correlation","VAF")
rownames(fitconj)<-paste("CONJUNCTIVE F=",seq(minfeat,maxfeat),sep="")
class(result)<-"summary.stepplfm"
result

}
}

#####################################
## function print.summary.stepplfm
#####################################

print.summary.stepplfm<-function(x,digits=2,...)
{
print(x[],digits=digits,...)
}

##################################################################
### function bayesplfm
##################################################################

bayesplfm<-function(data,object,attribute,rating,freq1,freqtot,F,Nchains=2,Nburnin=0,maxNiter=4000,Nstep=1000,Rhatcrit=1.2,maprule="disj",datatype="freq",start.bayes="best",fitted.plfm=NULL){

## functions

pi<-function(J,K,F,sigma,rho){
prob<-matrix(rep(1,J*K),nrow=J)
for (f in 1:F) {prob<-prob*(1-sigma[,f]%o%rho[,f])}
pi<-1-prob
}

post95<-function(x){post95<-quantile(x,c(0.025,0.5,0.975))}

logit<-function(x){logit<-log(x/(1-x))}

Rhat<-function(x){
nrow<-dim(x)[1]
ncol<-dim(x)[2]
gsample<-apply(x,2,mean)
varsample<-apply(x,2,var)
between<-nrow*var(gsample)
within<-mean(varsample)
varW<-var(varsample)/ncol
varB<-2*(between^2)/(ncol-1)
cov2<-cov(varsample,gsample)
muhat<-mean(gsample)
covWB<-(cov1-2*muhat*cov2)*(nrow/ncol)
postvar<-((nrow-1)/nrow)*within+(1/nrow)*(1+(1/ncol))*between
varpostvar<-((nrow-1)^2*varW +(1+(1/ncol))^2*varB +2*(nrow-1)*(1+(1/ncol))*covWB)/(nrow^2)
df<-2*(postvar^2)/varpostvar
Rhat<-sqrt((postvar/within)*((df+3)/(df+1)))
}

computedraw<-function(J,K,F,C,binmat,sigma,rho){

## compute P(X_jki=x_jki|d_ijk,theta)

probd.x0<-matrix(rep(1,C*K),nrow=C)
for (f in 1:F){probd.x0<-probd.x0*(1-binmat[,f]%o%rho[,f])}
nprobd.x0<-rep(1,J)%o%t(probd.x0)
nprobd.x1<-1-nprobd.x0

margprob.sigma<-exp(binmat%*%t(log(sigma))+(1-binmat)%*%log(1-t(sigma)))
nmargprob.sigma<-aperm(margprob.sigma%o%rep(1,K),c(2,3,1))

probx.d1<-nprobd.x1*nmargprob.sigma
probx.d0<-nprobd.x0*nmargprob.sigma

## normalize
norm1<-(apply(probx.d1,c(1,2),sum))%o%rep(1,C)
norm0<-(apply(probx.d0,c(1,2),sum))%o%rep(1,C)
probx.d1<-probx.d1/norm1
probx.d0<-probx.d0/norm0

## draw latent data objects
latx1<-array(rep(0,J*K*C),c(J,K,C))
latx0<-array(rep(0,J*K*C),c(J,K,C))
for (j in 1:J){
for (k in 1:K){
if (freq1[j,k]>0) latx1[j,k,]<-rmultinom(1,freq1[j,k],probx.d1[j,k,])
if (freq0[j,k]>0) latx0[j,k,]<-rmultinom(1,freq0[j,k],probx.d0[j,k,])
}}

statobj1<-apply(latx1+latx0,c(1,3),sum)%*%binmat
statobj0<-apply(latx1+latx0,c(1,3),sum)%*%(1-binmat)

## draw objectparameters
sigma.n<-matrix(rbeta(J*F,1+statobj1,1+statobj0),nrow=J)

## compute P(Y_kji=y_kji|X_jki=x_jki,d_ijk,theta)

probd.xy0<-matrix(rep(1,C*C),nrow=C)
for (f in 1:F){probd.xy0<-probd.xy0*(1-binmat[,f]%o%binmat[,f])}
nprobd.xy0<-rep(1,K)%o%probd.xy0
nprobd.xy1<-1-nprobd.xy0

margprob.rho<-t(exp(binmat%*%t(log(rho))+(1-binmat)%*%log(1-t(rho))))
nmargprob.rho<-aperm(margprob.rho%o%rep(1,C),c(1,3,2))

proby.xd1<-nprobd.xy1*nmargprob.rho
proby.xd0<-nprobd.xy0*nmargprob.rho

## normalize

norm1<-(apply(proby.xd1,c(1,2),sum))%o%rep(1,C)
norm0<-(apply(proby.xd0,c(1,2),sum))%o%rep(1,C)

proby.xd1<-proby.xd1/norm1
proby.xd0<-proby.xd0/norm0

proby.xd1[is.na(proby.xd1)]<-0
proby.xd0[is.na(proby.xd0)]<-0

Slatx1<-apply(latx1,c(2,3),sum)
Slatx0<-apply(latx0,c(2,3),sum)

## draw latent data attributes

laty1<-array(rep(0,K*C*C),c(K,C,C))
laty0<-array(rep(0,K*C*C),c(K,C,C))
for (k in 1:K){
for (c in 1:C){
if (Slatx1[k,c]>0) laty1[k,c,]<-rmultinom(1,Slatx1[k,c],proby.xd1[k,c,])
if (Slatx0[k,c]>0) laty0[k,c,]<-rmultinom(1,Slatx0[k,c],proby.xd0[k,c,])
}}

statatt1<-(apply(laty1+laty0,c(1,3),sum))%*%binmat
statatt0<-(apply(laty1+laty0,c(1,3),sum))%*%(1-binmat)

## draw attribute parameters
rho.n<-matrix(rbeta(K*F,1+statatt1,1+statatt0),nrow=K)

computedraw<-list(sigma.n=sigma.n,rho.n=rho.n)
}

## derive freq1 and freqtot  data from input

if (datatype=="freq") {
if (length(freqtot)==1) freqtot<-matrix(rep(freqtot,((dim(freq1)[1])*(dim(freq1)[2]))),nrow=dim(freq1)[1])
}

if (datatype=="dataframe"){
agg<-table(data\$object,data\$attribute,data\$rating)
freq1<-agg[,,"1"]
freqtot<-agg[,,"0"]+agg[,,"1"]}

## check input parameters, define starting values

if (sum(is.nan(freq1)|is.na(freq1)|freq1<0)){print("The elements of the matrix freq1 should be non-missing and larger than or equal to 0")}
else if (sum(is.nan(freqtot)|is.na(freqtot)|freqtot<1)){print("The elements of the matrix freqtot should be non-missing and larger than 0")}
else if (sum(abs(dim(freqtot)-dim(freq1)))>0){print("The matrices freq1 and freqtot should have the same dimensionality")}
else if (sum(freq1>freqtot)>0) {print("The corresponding values of the matrix freqtot should be larger or equal to the elements in the matrix freq1")}
else if ((maprule!="disj") & (maprule!="conj")) {print("The mapping rule is not correctly specified")}
else if (((dim(freq1)[1])*(dim(freq1)[2]))<F*sum(dim(freq1))) {print("The number of model parameters (J+K)*F should be smaller than the number of observations J*K used to fit the model")}
else if (F<1){print("The requested number of features should be larger than zero")}
else if (Nburnin<0){print("The requested number of burn-in iterations should be larger than or equal to zero")}
else if (Nchains<1){print("The requested number of chains should be larger than or equal to 1")}
else if (Nstep>maxNiter) {print("maxNiter should be larger than Nstep")}
else if (Rhatcrit<1.05){print("The convergence criterion Rhatcrit should not be smaller than 1.05")}

else {

J<-dim(freq1)[1]
K<-dim(freq1)[2]
C<-2**F
featurelabels<-paste("F",seq(1,F),sep="")
if (is.null(rownames(freq1))=="FALSE") rowlabels<-rownames(freq1)
else rowlabels<-paste("R_", seq(1,J),sep="")
if (is.null(colnames(freq1))=="FALSE") columnlabels<-colnames(freq1)
else columnlabels<-paste("C_", seq(1,K),sep="")

## define start values

if (!(is.null(fitted.plfm)) & (start.bayes=="fitted.plfm")){
if  ( (!(is.list(fitted.plfm)))) {print("fitted.plfm should be a list with components objpar and attpar")}
else if  ( (is.list(fitted.plfm)) & (length(fitted.plfm\$objpar)==0)) {print("The list fitted.plfm should contain the component objpar")}
else if  ( (is.list(fitted.plfm)) & (length(fitted.plfm\$attpar)==0)) {print("The list fitted.plfm should contain the component attpar")}
else if  ( (is.list(fitted.plfm)) & (sum(is.nan(fitted.plfm\$objpar)|is.na(fitted.plfm\$objpar)|fitted.plfm\$objpar<0|fitted.plfm\$objpar>1)>0)) {print("The values in the matrix fitted.plfm\$objpar should be non-missing and between 0 and 1")}
else if  ( (is.list(fitted.plfm)) & (sum(is.nan(fitted.plfm\$attpar)|is.na(fitted.plfm\$attpar)|fitted.plfm\$attpar<0|fitted.plfm\$attpar>1)>0)) {print("The values in the matrix fitted.plfm\$attpar should be non-missing and between 0 and 1")}
else if  ( (is.list(fitted.plfm)) & (dim(fitted.plfm\$objpar)[1]!=dim(freq1)[1])){print("The matrix fitted.plfm\$objpar should have the same number of rows as the matrix freq1")}
else if  ( (is.list(fitted.plfm)) & (dim(fitted.plfm\$attpar)[1]!=dim(freq1)[2])){print("The matrix number of rows of the matrix fitted.plfm\$attpar should equal the number of columns of the matrix freq1")}
else if  ( (is.list(fitted.plfm)) & (dim(fitted.plfm\$objpar)[2]!=F)){print(paste("The matrix fitted.plfm\$objpar should have ", F," columns",sep=""))}
else if ( (is.list(fitted.plfm)) & (dim(fitted.plfm\$attpar)[2]!=F)){print(paste("The matrix fitted.plfm\$attpar should have ", F," columns",sep=""))}
else {
start.objpar<-fitted.plfm\$objpar
start.attpar<-fitted.plfm\$attpar}
}

else if (start.bayes=="best")
{
postmode<-plfm(maprule=maprule,object=object,attribute=attribute,rating=rating,freq1=freq1,freqtot=freqtot,datatype=datatype,F=F,M=20,emcrit1=1e-2,emcrit2=1e-10,printrun="FALSE")
start.objpar<-postmode\$objpar
start.attpar<-postmode\$attpar
}

else if (start.bayes=="random")
{
start.objpar<-matrix(runif(J*F),nrow=J)
start.attpar<-matrix(runif(K*F),nrow=K)
}

initial.objpar<-array(rep(0,J,F,Nchains),c(J,F,Nchains))
for (chain in 1:Nchains){initial.objpar[,,chain]<-start.objpar}
initial.attpar<-array(rep(0,K,F,Nchains),c(K,F,Nchains))
for (chain in 1:Nchains){initial.attpar[,,chain]<-start.attpar}

freq0<-freqtot-freq1
if (maprule=="conj"){
freq0<-freq1
freq1<-freqtot-freq0
initial.objpar<-1-initial.objpar
}

## generate binary matrices
binmat<-matrix(rep(0,C*F),nrow=C)
for (i in 1:(C-1)){binmat[i+1,]<-c(digitsBase(i,base=2,F))}

## define output objects

sample.sigma<-array(rep(0,J*F*maxNiter*Nchains),c(J,F,maxNiter,Nchains))
dimnames(sample.sigma)[[1]]<-rowlabels
dimnames(sample.sigma)[[2]]<-featurelabels

sample.rho<-array(rep(0,K*F*maxNiter*Nchains),c(K,F,maxNiter,Nchains))
dimnames(sample.rho)[[1]]<-columnlabels
dimnames(sample.rho)[[2]]<-featurelabels

for (chain in 1:Nchains){

## define initial values

sigma.o<-as.matrix(initial.objpar[,,chain])
rho.o<-as.matrix(initial.attpar[,,chain])

for (iter in 1:Nburnin){

##if (iter%%100==0) print(c("iteration ",iter),justify="left",quote="FALSE")
newdraw<-computedraw(J,K,F,C,binmat,sigma.o,rho.o)
sigma.o<-newdraw\$sigma.n
rho.o<-newdraw\$rho.n
}

for (iter in 1:Nstep){
##if (iter%%1000==0) print(c("chain ",chain,"iteration ",iter),justify="left",quote="FALSE")
newdraw<-computedraw(J,K,F,C,binmat,sigma.o,rho.o)
sigma.o<-newdraw\$sigma.n
rho.o<-newdraw\$rho.n
sample.sigma[,,iter,chain]<-newdraw\$sigma.n
sample.rho[,,iter,chain]<-newdraw\$rho.n
}
}

## permute features to avoid label switching

if ((Nchains>1) & (F>1)){
temp<-abind(sample.sigma[,,(1:Nstep),],sample.rho[,,(1:Nstep),],along=1)
gtemp<-apply(temp,c(1,2,4),mean)
new<-temp
for (chain in 2:Nchains){
new[,,,chain]<-temp[,c(apply(cor(gtemp[,,1],gtemp[,,chain]),1,order)[F,]),,chain]}
sample.sigma[,,(1:Nstep),]<-new[(1:J),,,]
sample.rho[,,(1:Nstep),]<-new[((J+1):(J+K)),,,]
}

## compute convergence

if (Nchains==1)
{convstat<-(J+K)*F
print(c("chainlength",Nstep),justify="left",quote="FALSE")}

if (Nchains>1){
Rhat.sigma<-matrix(rep(0,J*F),nrow=J)
Rhat.rho<-matrix(rep(0,K*F),nrow=K)
rownames(Rhat.sigma)<-rowlabels
colnames(Rhat.sigma)<-featurelabels
rownames(Rhat.rho)<-columnlabels
colnames(Rhat.rho)<-featurelabels

for (j in 1:J){
for (f in 1:F){
mat<-logit(sample.sigma[j,f,(1:Nstep),])
Rhat.sigma[j,f]<-Rhat(mat)
}}

for (k in 1:K){
for (f in 1:F){
mat<-logit(sample.rho[k,f,(1:Nstep),])
Rhat.rho[k,f]<-Rhat(mat)
}}

convstat<-sum(Rhat.sigma>Rhatcrit)+sum(Rhat.rho>Rhatcrit)
print(c("chainlength",Nstep,"# parameters to converge",convstat),justify="left",quote="FALSE")
}

stepnumber<-1

while ((convstat>0) & (((stepnumber+1)*Nstep)<=maxNiter)){

## define output step

for (chain in 1:Nchains){

## define initial values

sigma.o<-sample.sigma[,,(stepnumber*Nstep),chain]
rho.o<-sample.rho[,,(stepnumber*Nstep),chain]

for (iter in 1:Nstep){
newdraw<-computedraw(J,K,F,C,binmat,sigma.o,rho.o)
sigma.o<-newdraw\$sigma.n
rho.o<-newdraw\$rho.n
sample.sigma[,,(stepnumber*Nstep+iter),chain]<-newdraw\$sigma.n
sample.rho[,,(stepnumber*Nstep+iter),chain]<-newdraw\$rho.n
}
}

stepnumber<-stepnumber+1

## permute features to avoid label switching

if ((Nchains>1) & (F>1)){
temp<-abind(sample.sigma[,,(1:(stepnumber*Nstep)),],sample.rho[,,(1:(stepnumber*Nstep)),],along=1)
gtemp<-apply(temp,c(1,2,4),mean)
new<-temp
for (chain in 2:Nchains){
new[,,,chain]<-temp[,c(apply(cor(gtemp[,,1],gtemp[,,chain]),1,order)[F,]),,chain]}
sample.sigma[,,(1:(stepnumber*Nstep)),]<-new[(1:J),,,]
sample.rho[,,(1:(stepnumber*Nstep)),]<-new[((J+1):(J+K)),,,]
}

## compute convergence

if (Nchains==1)
{print(c("chainlength",(stepnumber*Nstep)),justify="left",quote="FALSE")}

if (Nchains>1){
Rhat.sigma<-matrix(rep(0,J*F),nrow=J)
Rhat.rho<-matrix(rep(0,K*F),nrow=K)
rownames(Rhat.sigma)<-rowlabels
colnames(Rhat.sigma)<-featurelabels
rownames(Rhat.rho)<-columnlabels
colnames(Rhat.rho)<-featurelabels

for (j in 1:J){
for (f in 1:F){
mat<-logit(sample.sigma[j,f,(1:(stepnumber*Nstep)),])
Rhat.sigma[j,f]<-Rhat(mat)
}}

for (k in 1:K){
for (f in 1:F){
mat<-logit(sample.rho[k,f,(1:(stepnumber*Nstep)),])
Rhat.rho[k,f]<-Rhat(mat)
}}

convstat<-sum(Rhat.sigma>Rhatcrit)+sum(Rhat.rho>Rhatcrit)
print(c("chainlength",(stepnumber*Nstep),"# parameters to converge",convstat),justify="left",quote="FALSE")

}
}

##compute posterior mean and posterior intervals
if (F==1)
{
pmean.sigma<-as.matrix(apply(sample.sigma[,,(1:(stepnumber*Nstep)),],1,mean))
pmean.rho<-as.matrix(apply(sample.rho[,,(1:(stepnumber*Nstep)),],1,mean))
p95.sigma<-as.matrix(apply(sample.sigma[,,(1:(stepnumber*Nstep)),],1,post95))
p95.rho<-as.matrix(apply(sample.rho[,,(1:(stepnumber*Nstep)),],1,post95))
}

if (F>1)
{
pmean.sigma<-apply(sample.sigma[,,(1:(stepnumber*Nstep)),],c(1,2),mean)
pmean.rho<-apply(sample.rho[,,(1:(stepnumber*Nstep)),],c(1,2),mean)
p95.sigma<-apply(sample.sigma[,,(1:(stepnumber*Nstep)),],c(1,2),post95)
p95.rho<-apply(sample.rho[,,(1:(stepnumber*Nstep)),],c(1,2),post95)
}

## compute fitmeasures

p1.mean<-pi(J,K,F,pmean.sigma,pmean.rho)
correlation<-cor(c(p1.mean*freqtot),c(freq1))
VAF<-correlation**2
fitmeasures<-as.matrix(c(correlation,VAF))
rownames(fitmeasures)<-c("Correlation observed and expected frequencies","VAF observed frequencies")
colnames(fitmeasures)<-""

if ((maprule=="disj") & (Nchains>1))
bayesplfm<-list(call=match.call(),sample.objpar=sample.sigma[,,(1:(stepnumber*Nstep)),],sample.attpar=sample.rho[,,(1:(stepnumber*Nstep)),],pmean.objpar=pmean.sigma,pmean.attpar=pmean.rho,p95.objpar=p95.sigma,p95.attpar=p95.rho,Rhat.objpar=Rhat.sigma,Rhat.attpar=Rhat.rho,fitmeasures=fitmeasures,convstat=convstat)
else if ((maprule=="disj") & (Nchains==1))
bayesplfm<-list(call=match.call(),sample.objpar=sample.sigma[,,(1:(stepnumber*Nstep)),],sample.attpar=sample.rho[,,(1:(stepnumber*Nstep)),],pmean.objpar=pmean.sigma,pmean.attpar=pmean.rho,p95.objpar=p95.sigma,p95.attpar=p95.rho,fitmeasures=fitmeasures)
else if ((maprule=="conj") & (Nchains>1))
bayesplfm<-list(call=match.call(),sample.objpar=1-sample.sigma[,,(1:(stepnumber*Nstep)),],sample.attpar=sample.rho[,,(1:(stepnumber*Nstep)),],pmean.objpar=1-pmean.sigma,pmean.attpar=pmean.rho,p95.objpar=1-p95.sigma,p95.attpar=p95.rho,Rhat.objpar=Rhat.sigma,Rhat.attpar=Rhat.rho,fitmeasures=fitmeasures,convstat=convstat)
else if ((maprule=="conj") & (Nchains==1))
bayesplfm<-list(call=match.call(),sample.objpar=1-sample.sigma[,,(1:(stepnumber*Nstep)),],sample.attpar=sample.rho[,,(1:(stepnumber*Nstep)),],pmean.objpar=1-pmean.sigma,pmean.attpar=pmean.rho,p95.objpar=1-p95.sigma,p95.attpar=p95.rho,fitmeasures=fitmeasures)

class(bayesplfm)<-"bayesplfm"
bayesplfm
}
}

############################
## function print.bayesplfm
############################

print.bayesplfm <- function(x, ...)
{
J<-dim(x\$pmean.objpar)[1]
K<-dim(x\$pmean.attpar)[1]
F<-dim(x\$pmean.objpar)[2]
Nchains<-dim(x\$sample.objpar)[4]
Npar<-(J+K)*F

cat("CALL:\n")
print(x\$call)
if (is.na(Nchains)=="FALSE"){
cat("\nNUMBER OF PARAMETERS THAT DO NOT MEET CONVERGENCE CRITERION:\n")
convstat<-matrix(c(Npar,x\$convstat),ncol=1)
colnames(convstat)<-""
rownames(convstat)<-c("total number of parameters","number of parameters without convergence")
print(convstat)
}

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
print(round(x\$fitmeasures,3))

cat("\nPOSTERIOR MEAN OBJECTPARAMETERS:\n")
strout<-as.matrix(capture.output(round(x\$pmean.objpar,2)))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$pmean.objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\nPOSTERIOR MEAN ATTRIBUTEPARAMETERS:\n")
strout<-as.matrix(capture.output(round(x\$pmean.attpar,2)))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$pmean.attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")
}

###############################
## function summary.bayesplfm
###############################

summary.bayesplfm <- function(object, ...)
{
J<-dim(object\$pmean.objpar)[1]
K<-dim(object\$pmean.attpar)[1]
F<-dim(object\$pmean.objpar)[2]
Nchains<-dim(object\$sample.objpar)[4]

descript<-round(object\$fitmeasures,3)
estobjpar<-round(object\$pmean.objpar,2)
estattpar<-round(object\$pmean.attpar,2)
p95objpar<-matrix(rep("",J*F),nrow=J)
rownames(p95objpar)<-rownames(object\$pmean.objpar)
colnames(p95objpar)<-colnames(object\$pmean.objpar)
for (f in 1:F) p95objpar[,f]<-paste("[",round(object\$p95.objpar[1,,f],3),";",round(object\$p95.objpar[3,,f],3),"]",sep="")
p95attpar<-matrix(rep("",K*F),nrow=K)
rownames(p95attpar)<-rownames(object\$pmean.attpar)
colnames(p95attpar)<-colnames(object\$pmean.attpar)
for (f in 1:F) p95attpar[,f]<-paste("[",round(object\$p95.attpar[1,,f],3),";",round(object\$p95.attpar[3,,f],3),"]",sep="")
if (is.na(Nchains)){
result<-list(call=object\$call,
descriptivefit=descript,
objpar=estobjpar,
attpar=estattpar,
p95objpar=p95objpar,
p95attpar=p95attpar)
}
else if (Nchains>1){
Rhatobjpar<-round(object\$Rhat.objpar,3)
Rhatattpar<-round(object\$Rhat.attpar,3)
result<-list(call=object\$call,
descriptivefit=descript,
objpar=estobjpar,
attpar=estattpar,
p95objpar=p95objpar,
p95attpar=p95attpar,
Rhatobjpar=Rhatobjpar,
Rhatattpar=Rhatattpar)
}
class(result)<-"summary.bayesplfm"
result
}

####################################
## function print.summary.bayesplfm
####################################

print.summary.bayesplfm<-function(x, ...)
{
cat("Call:\n")
print(x\$call)

cat("\nDESCRIPTIVE FIT OBJECT X ATTRIBUTE TABLE:\n")
print(x\$descriptivefit)

cat("\nPOSTERIOR MEAN OBJECTPARAMETERS:\n")
strout<-as.matrix(capture.output(x\$objpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\n95% POSTERIOR INTERVAL OBJECTPARAMETERS:\n")
strout<-as.matrix(capture.output(x\$p95objpar))
strout<-gsub("\"", "", strout,fixed="TRUE")
strout<-gsub("[0.", "[.", strout,fixed="TRUE")
strout<-gsub(";0.", ";.", strout,fixed="TRUE")
rownames(strout)<-rep("",(1+dim(x\$p95objpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

if (is.null(x\$Rhatobjpar)=="FALSE")
{cat("\n RHAT CONVERGENCE OBJECTPARAMETERS:\n")
print(x\$Rhatobjpar)}

cat("\nPOSTERIOR MEAN ATTRIBUTEPARAMETERS:\n")
strout<-as.matrix(capture.output(x\$attpar))
strout<- gsub(" 0", "  ", strout)
rownames(strout)<-rep("",(1+dim(x\$attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

cat("\n95% POSTERIOR INTERVAL ATTRIBUTEPARAMETERS:\n")
strout<-as.matrix(capture.output(x\$p95attpar))
strout<-gsub("\"", "", strout,fixed="TRUE")
strout<-gsub("[0.", "[.", strout,fixed="TRUE")
strout<-gsub(";0.", ";.", strout,fixed="TRUE")
rownames(strout)<-rep("",(1+dim(x\$p95attpar)[1]))
colnames(strout)<-rep("",1)
print(strout,quote="FALSE")

if (is.null(x\$Rhatattpar)=="FALSE")
{cat("\n RHAT CONVERGENCE ATTRIBUTEPARAMETERS:\n")
print(x\$Rhatattpar)}
}

#############################
## function gendatLCplfm
############################

gendatLCplfm<-function(N,objpar,attpar,sizepar,maprule="disj",model=1){

T<-length(sizepar)
J<-dim(objpar)[1]
K<-dim(attpar)[1]
F<-dim(objpar)[2]

if (sum((objpar<0)|(objpar>1)|is.nan(objpar))>0) print("Object parameters should all be non-missing and between 0 and 1")
else if (sum((attpar<0)|(attpar>1)|is.nan(attpar))>0) print("Attribute parameters should all be between 0 and 1")
else if ((maprule!="disj")&(maprule!="conj")&(maprule!="add")) print("Incorrect specification of mapping rule")
else if ((model!=1) & (model!=2)  & (model!=3) & (model!=4) & (model!=5) & (model!=6)) print("model should be specified as 1, 2, 3, 4, 5 or 6")
else if (N<=0) print("The number of replications should be positive (N>0)")
else if (N>0)
{

## draw latent class membership
lc<-rmultinom(N,1,sizepar)
lcnum<-t(lc)%*%c(1:T)

## create empty output data
data<-array(rep(0,N*J*K),c(N,J,K))

if (model==1){
for (i in 1:N){
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,,lcnum[i]])
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f, t] %o% attpar[, f])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f, t]) %o% attpar[, f])}}
}

condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f, t] %o% attpar[, f])}}
condprob.JKT<-condprob.JKT/F
}
} ## end model==1

if (model==2){
for (i in 1:N){
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,])
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,,lcnum[i]])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f] %o% attpar[, f, t])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f]) %o% attpar[, f, t])}}
}

condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f] %o% attpar[, f, t])}}
condprob.JKT<-condprob.JKT/F
}

} ## end model==2

if (model==3){
for (i in 1:N){
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,,lcnum[i]])
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,,lcnum[i]])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f, t] %o% attpar[, f, t])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f, t]) %o% attpar[, f, t])}}
}
condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f, t] %o% attpar[, f, t])}}
condprob.JKT<-condprob.JKT/F
}
} ## end model==3

if (model==4){
for (i in 1:N){
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,])
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,,lcnum[i]])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f, t] %o% attpar[, f])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f, t]) %o% attpar[, f])}}
}

condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f, t] %o% attpar[, f])}}
condprob.JKT<-condprob.JKT/F
}

} ## end model==4

if (model==5){
for (i in 1:N){
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,,lcnum[i]])
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f] %o% attpar[, f, t])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f]) %o% attpar[, f, t])}}
}

condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f] %o% attpar[, f, t])}}
condprob.JKT<-condprob.JKT/F
}

} ## end model==5

if (model==6){
for (i in 1:N){
for (k in 1:K){
laty<-rbinom(F,1,attpar[k,,lcnum[i]])
for (j in 1:J){
latx<-rbinom(F,1,objpar[j,,lcnum[i]])
if (maprule=="disj"){data[i,j,k]<-ifelse(sum(latx*laty)>0,1,0)}
else if (maprule=="conj") {data[i,j,k]<-ifelse(min(latx-laty)<0,0,1)}
else if (maprule=="add"){if (runif(1)<=mean(latx*laty)) data[i,j,k]<-1 else data[i,j,k]<-0}
}}}

## compute conditional probabilities
if (maprule=="disj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1 - objpar[, f, t] %o% attpar[, f, t])}}
condprob.JKT <- 1 - condprob.JKT
}
else if (maprule=="conj"){
condprob.JKT <- array(rep(1, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] *
(1-(1 - objpar[, f, t]) %o% attpar[, f, t])}}
}

condprob.JKT <- array(rep(0, J * K * T), c(J, K, T))
for (t in 1:T) {
for (f in 1:F) {
condprob.JKT[, , t] <- condprob.JKT[, , t] +
(objpar[, f, t] %o% attpar[, f, t])}}
condprob.JKT<-condprob.JKT/F
}

} ## end model==6

## compute marginal association probabilities
margprob.JK <- matrix(rep(1, J * K), nrow = J)
weight <- rep(1, J) %o% rep(1, K) %o% sizepar
margprob.JK <- apply(weight * condprob.JKT, c(1, 2),sum)

gendatLCplfm<-list(call=match.call(),data=data,class=lcnum,condprob.JKT=condprob.JKT,margprob.JK=margprob.JK)
}
}

####################
## function LCplfm
####################

LCplfm<-function(data,F=2,T=2,M=5,maprule="disj",emcrit1=1e-3,emcrit2=1e-8,model=1,
start.objectparameters=NULL,start.attributeparameters=NULL,
start.sizeparameters=NULL,delta=0.0001,printrun=FALSE,update.objectparameters=NULL,update.attributeparameters=NULL, Nbootstrap=2000)
{

##functions
post95<-function(x){post95<-quantile(x,c(0.025,0.5,0.975))}
post99<-function(x){post99<-quantile(x,c(0.005,0.5,0.995))}

# check input and analysis specifications
# dimensions data

NI<-dim(data)[1]
NJ<-dim(data)[2]
NK<-dim(data)[3]

STOBJ<-rep(0,4)
if (!is.null(start.objectparameters)){
for (i in (1:length(dim(start.objectparameters)))){STOBJ[i]<-dim(start.objectparameters)[i]}
}

STATT<-rep(0,4)
if (!is.null(start.attributeparameters)){
for (i in (1:length(dim(start.attributeparameters)))){STATT[i]<-dim(start.attributeparameters)[i]}
}

STSIZE<-rep(0,2)
if (!is.null(start.sizeparameters)){
for (i in (1:length(dim(start.sizeparameters)))){STSIZE[i]<-dim(start.sizeparameters)[i]}
}

UPOBJ<-rep(0,3)
if (!is.null(update.objectparameters)){
for (i in (1:length(dim(update.objectparameters)))){UPOBJ[i]<-dim(update.objectparameters)[i]}
}

UPATT<-rep(0,3)
if (!is.null(update.attributeparameters)){
for (i in (1:length(dim(update.attributeparameters)))){UPATT[i]<-dim(update.attributeparameters)[i]}
}

if (length(dim(data))!=3|(NI<5)) {print("The array 'data' should be a three-way array with sufficient replications (I>5)")}
else if (sum((data!=0) & (data!=1))>0){print("The elements of the array 'data' should equal 0 or 1")}
else if (F<1){print("The number of latent features should be at least 1 (F>=1)")}
else if (T<1){print("The number of latent classes should be at least 1 (T>=1)")}
else if ((maprule!="disj") & (maprule!="conj") & (maprule!="add")) {print("The mapping rule is not correctly specified")}
else if (M<1) {print("The number of requested runs should be at least one (M>=1)")}
else if (emcrit1<0 | emcrit2<0) {print("emcrit1 and emcrit2 should be small positive numbers")}

else if  ((model==1| model==4| model==3| model==6) & (!is.null(start.objectparameters)) & (STOBJ[1] != NJ| STOBJ[2]!=F | STOBJ[3]!=T | STOBJ[4]!=M | sum(start.objectparameters==0|start.objectparameters==1)>0))
{print("When model=1,4,3,6 start.objectparameters should be either NULL or a J x F x T x M array. The boundary values 0,1 cannot be used as starting values.")}

else if  ((model==2|model==5) & (!is.null(start.objectparameters)) & (STOBJ[1] != NJ | STOBJ[2]!=F | STOBJ[3]!=M | sum(start.objectparameters==0|start.objectparameters==1)>0))
{print("When model=2,5 start.objectparameters should be either NULL or a J x F x M array. The boundary values 0,1 cannot be used as starting values.")}

else if  ((model==2|model==5|model==3|model==6) & (!is.null(start.attributeparameters)) & (STATT[1] != NK| STATT[2]!=F | STATT[3]!=T | STATT[4]!=M | sum(start.attributeparameters==0|start.attributeparameters==1)>0))
{print("When model=2,5,3,6 start.attributeparameters should be either NULL or a K x F x T x M array. The boundary values 0,1 cannot be used as starting values.")}

else if  ((model==1|model==4) & (!is.null(start.attributeparameters)) & (STATT[1] != NK | STATT[2]!=F | STATT[3]!=M | sum(start.attributeparameters==0|start.attributeparameters==1)>0))
{print("When model=1,4 start.attributeparameters should be either NULL or a K x F x M array. The boundary values 0,1 cannot be used as starting values.")}

else if  ((!is.null(start.sizeparameters)) & (STSIZE[1] != T | STSIZE[2]!=M))
{print("Start.sizeparameters should be either NULL or a T x M matrix")}

else if  ((model==1|model==4|model==3|model==6) & (!is.null(update.objectparameters)) & (UPOBJ[1] != NJ| UPOBJ[2]!=F | UPOBJ[3]!=T | sum(update.objectparameters!=0 & update.objectparameters!=1)>0))
{print("When model=1,4,3,6 update.objectparameters should be either NULL or a binary valued (0/1) J x F x T array")}

else if  ((model==2|model==5) & (!is.null(update.objectparameters)) & (UPOBJ[1] != NJ | UPOBJ[2]!=F | UPOBJ[3]!=M | sum(update.objectparameters!=0 & update.objectparameters!=1)>0))
{print("When model=2,5 update.objectparameters should be either NULL or a binary valued (0/1) J x F matrix")}

else if  ((model==2|model==5|model==3|model==6) & (!is.null(update.attributeparameters)) & (UPATT[1] != NK | UPATT[2]!=F | UPATT[3]!=T | sum(update.attributeparameters!=0 & update.attributeparameters!=1)>0))
{print("When model=2,5,3,6 update.attributeparameters should be either NULL or a binary valued (0/1) K x F x T array")}

else if  ((model==1|model==4) & (!is.null(update.attributeparameters)) & (UPATT[1] != NK | UPATT[2]!=F | sum(update.attributeparameters!=0 & update.attributeparameters!=1)>0))
{print("When model=1,4 update.attributeparameters should be either NULL or a binary valued (0/1) K x F matrix")}

## estimate model M_1
## LCPLFM with constant classification objects and heterogeneity object parameters

else if (model==1){

if (maprule=="conj") {data<-1-data}

## data should be an IxJxK array
I<-dim(data)[1]
J<-dim(data)[2]
K<-dim(data)[3]

# compute number of object/attribute pairs
Nattpair<-K*(K-1)/2
Nobjpair<-J*(J-1)/2

## define labels

if (is.null(dimnames(data)[[1]])=="FALSE") raterlabels<-dimnames(data)[[1]]
else raterlabels<-paste("RATER_", seq(1,I),sep="")

if (is.null(dimnames(data)[[2]])=="FALSE") objectlabels<-dimnames(data)[[2]]
else objectlabels<-paste("O_", seq(1,J),sep="")

if (is.null(dimnames(data)[[3]])=="FALSE") attributelabels<-dimnames(data)[[3]]
else attributelabels<-paste("A_", seq(1,K),sep="")

featurelabels<-paste("F",seq(1,F),sep="")
classlabels<-paste("T",seq(1,T),sep="")
runlabels<-paste("RUN",seq(1,M),sep="")

##transform data to be used by c++ function
ndata<-c(aperm(data,c(1,3,2)))

## define constants
S<-2^F
##Npar<-J*F*T+K*F+T-1

## generate feature pattern matrix
binmat<-matrix(rep(0,S*F),nrow=S)
for (i in 1:(S-1)){binmat[i+1,]<-c(digitsBase(i,base=2,F))}
Pat <- c(t(binmat))

## define objects for output storage

sigma.runs <- array(rep(0, J * F * T* M), c(J, F, T, M))
rho.runs <- array(rep(0, K * F * M), c(K, F, M))
gamma.runs<-array(rep(0, T * M), c(T, M))
logpost.runs <- rep(0, M)
names(logpost.runs)<-runlabels
loglik.runs <-rep(0, M)
names(loglik.runs)<-runlabels

## generate random start values for M runs

start.obj<-array(runif(J*F*T*M),c(J,F,T,M))
start.att<-array(runif(K*F*M),c(K,F,M))
start.size<-array(rep(1/T,T*M),c(T,M))

for (run in 1:M){

## define starting values for each run of algorithm

if (is.null(start.objectparameters)){sigma.n<-c(start.obj[,,,run])}
else {sigma.n<-c(aperm(start.objectparameters[,,,run,drop="FALSE"],c(4,3,2,1)))}

if (is.null(start.attributeparameters)){rho.n<-c(start.att[,,run])}
else {rho.n<-c(t(start.attributeparameters[,,run]))}

if (is.null(start.sizeparameters)){gamma.n<-start.size[,run]}
else {gamma.n<-start.sizeparameters[,run]}

## define parameters to update

if (is.null(update.objectparameters)){
sigma.update<-c(array(rep(1,J*F*T),c(J,F,T)))}
else {sigma.update<-c(aperm(update.objectparameters,c(3,2,1)))}

if (is.null(update.attributeparameters)){
rho.update<-c(array(rep(1,K*F),c(K,F)))}
else {rho.update<-c(t(update.attributeparameters))}

## compute number of parameters
Npar<-sum(sigma.update)+sum(rho.update)+T-1

SEsigma <- double(J*F*T)
SErho <- double (K*F)
SEgamma <- double (T)

logpost.n <- 0.0
loglik.n <- 0.0

postprob <- double (T*I)

OR.commonatt.rep<-double(Nbootstrap*J*Nattpair)
OR.commonobj.rep<-double(Nbootstrap*K*Nobjpair)
report.commonatt<-double(J*Nattpair*5)
report.commonobj<-double(K*Nobjpair*5)

## call C function to estimate candidate model
## use flag==0 to skip computation of gradient and standard errors
## use convergence criterion emcrit1

if (maprule=="disj"|maprule=="conj"){
result <- .C("PlFm_XZ_Y_DC", as.integer(ndata), as.integer(I),as.integer(J),as.integer(K),
as.integer(F),as.integer(T),as.integer(Pat),as.double(emcrit1),as.double(sigma.n),
as.double(delta),as.double(logpost.n),as.double(loglik.n),as.double(postprob),as.integer(0),as.integer(Nbootstrap),
as.double(OR.commonatt.rep),as.double(OR.commonobj.rep),as.double(report.commonatt),as.double(report.commonobj))
}
as.integer(F),as.integer(T),as.integer(Pat),as.double(emcrit1),as.double(sigma.n),
as.double(delta),as.double(logpost.n),as.double(loglik.n),as.double(postprob),as.integer(0),as.integer(Nbootstrap),
as.double(OR.commonatt.rep),as.double(OR.commonobj.rep),as.double(report.commonatt),as.double(report.commonobj))
}

## take output from result
sigma.runs[,,,run]<-aperm(array(result[[9]],c(T,F,J)),c(3,2,1))
rho.runs[,,run]<-matrix(result[[10]][1:(K*F)],nrow=K,byrow=T)
gamma.runs[,run]<-result[[11]]
logpost.runs[run]<-result[[21]]
loglik.runs[run]<-result[[22]]

if (printrun=="TRUE"){
if (maprule=="disj") {print(paste("DISJUNCTIVE ANALYSIS  ","F=",F,"  T=",T,"  RUN=",run,sep=""),quote="FALSE")}
else if (maprule=="conj") {print(paste("CONJUNCTIVE ANALYSIS  ","F=",F,"  T=",T,"  RUN=",run,sep=""),quote="FALSE")}

}

} ##end run

## run final analysis with parameters of best solution as starting point

best<-order(logpost.runs,decreasing=T)[1]
if (F==1) {sigma.n<-array(sigma.runs[,,,best],c(J,F,T))} else if (F>1) {sigma.n<-sigma.runs[,,,best]}
if (F==1){rho.n<-as.matrix(rho.runs[,,best])} else if (F>1) {rho.n<-rho.runs[,,best]}
gamma.n<-gamma.runs[,best]

if ((F==1) & (T==1) & (J>1)) {sigma.n<-c(sigma.n)}
if ((F==1) & (T>1) & (J>1)) {sigma.n<-c(aperm(sigma.n,c(3,2,1)))}
if ((F>1) & (T==1) & (J>1)) {sigma.n<-c(aperm(sigma.n,c(2,1)))}
if ((F>1) & (T>1) & (J>1)) {sigma.n<-c(aperm(sigma.n,c(3,2,1)))}

if ((F==1) & (T==1) & (J==1)) {sigma.n<-c(sigma.n)}
if ((F==1) & (T>1) & (J==1)) {sigma.n<-c(sigma.n)}
if ((F>1) & (T==1) & (J==1)) {sigma.n<-c(sigma.n)}
if ((F>1) & (T>1) & (J==1)) {sigma.n<-c(aperm(sigma.n,c(2,1)))}

rho.n<-c(t(rho.n))
gamma.n<-c(gamma.n)

SEsigma <- double(J*F*T)
SErho <- double (K*F)
SEgamma <- double (T)

logpost.n <- 0.0
loglik.n <- 0.0

postprob <- double (T*I)

## call C function to estimate final model
## use flag==1 to compute gradient and standard errors,
## use emcrit2 as convergence criterion

if (maprule=="disj"|maprule=="conj"){
result <- .C("PlFm_XZ_Y_DC", as.integer(ndata), as.integer(I),as.integer(J),as.integer(K),
as.integer(F),as.integer(T),as.integer(Pat),as.double(emcrit2),as.double(sigma.n),
as.double(delta),as.double(logpost.n),as.double(loglik.n),as.double(postprob),as.integer(1),as.integer(Nbootstrap),
as.double(OR.commonatt.rep),as.double(OR.commonobj.rep),as.double(report.commonatt),as.double(report.commonobj))
}

as.integer(F),as.integer(T),as.integer(Pat),as.double(emcrit2),as.double(sigma.n),
as.double(delta),as.double(logpost.n),as.double(loglik.n),as.double(postprob),as.integer(1),as.integer(Nbootstrap),
as.double(OR.commonatt.rep),as.double(OR.commonobj.rep),as.double(report.commonatt),as.double(report.commonobj))
}

## take output from result
sigma.best<-aperm(array(result[[9]],c(T,F,J)),c(3,2,1))
rho.best<-matrix(result[[10]][1:(K*F)],nrow=K,byrow=T)
gamma.best<-result[[11]]

dimnames(sigma.best)[[1]]<-as.list(objectlabels)
dimnames(sigma.best)[[2]]<-as.list(featurelabels)
dimnames(sigma.best)[[3]]<-as.list(classlabels)
dimnames(rho.best)[[1]]<-as.list(attributelabels)
dimnames(rho.best)[[2]]<-as.list(featurelabels)
names(gamma.best)<-classlabels

SEsigma.best<-aperm(array(result[[17]],c(T,F,J)),c(3,2,1))
SEsigma.best[update.objectparameters==0]<-0

SErho.best<-matrix(result[[18]][1:(K*F)],nrow=K,byrow=T)
SErho.best[update.attributeparameters==0]<-0

if (T==1){SEgamma.best<-0} else if (T>1) {SEgamma.best<-result[[19]]}

dimnames(SEsigma.best)[[1]]<-as.list(objectlabels)
dimnames(SEsigma.best)[[2]]<-as.list(featurelabels)
dimnames(SEsigma.best)[[3]]<-as.list(classlabels)
dimnames(SErho.best)[[1]]<-as.list(attributelabels)
dimnames(SErho.best)[[2]]<-as.list(featurelabels)
names(SEgamma.best)<-classlabels

logpost.best<-result[[21]]
loglik.best<-result[[22]]

postprob.best<-matrix(result[[23]],nrow=I)
colnames(postprob.best)<-classlabels

OR.comatt.rep<-aperm(array(result[[26]],c(Nattpair,J,Nbootstrap)),c(3,2,1))
OR.comobj.rep<-aperm(array(result[[27]],c(Nobjpair,K,Nbootstrap)),c(3,2,1))

report.comatt<-matrix(result[[28]],nrow=J*Nattpair,byrow=TRUE)
report.comobj<-matrix(result[[29]],nrow=K*Nobjpair,byrow=TRUE)

#report bootstrap

p95.ORcomatt<-apply(OR.comatt.rep,c(2,3),post95)
p99.ORcomatt<-apply(OR.comatt.rep,c(2,3),post99)

p95.ORcomobj<-apply(OR.comobj.rep,c(2,3),post95)
p99.ORcomobj<-apply(OR.comobj.rep,c(2,3),post99)

temp<-cbind(c(t(p95.ORcomatt[1,,])),c(t(p95.ORcomatt[3,,])),c(t(p99.ORcomatt[1,,])),c(t(p99.ORcomatt[3,,])))
report.OR.attpair<-cbind(report.comatt,temp)
dimnames(report.OR.attpair)[[2]]<-c("object j","attribute k1","attribute k2","OR.obs","OR.mean","OR.p025","OR.p975","OR.p005","OR.p995")

temp<-cbind(c(t(p95.ORcomobj[1,,])),c(t(p95.ORcomobj[3,,])),c(t(p99.ORcomobj[1,,])),c(t(p99.ORcomobj[3,,])))
report.OR.objpair<-cbind(report.comobj,temp)
dimnames(report.OR.objpair)[[2]]<-c("attribute k","object j1","object j2","OR.obs","OR.mean","OR.p025","OR.p975","OR.p005","OR.p995")

pOR.attpair.inCI95<-mean(ifelse((report.OR.attpair[,4]>= report.OR.attpair[,6])& (report.OR.attpair[,4]<= report.OR.attpair[,7]),1,0))
pOR.attpair.inCI99<-mean(ifelse((report.OR.attpair[,4]>= report.OR.attpair[,8])& (report.OR.attpair[,4]<= report.OR.attpair[,9]),1,0))

pOR.objpair.inCI95<-mean(ifelse((report.OR.objpair[,4]>= report.OR.objpair[,6])& (report.OR.objpair[,4]<= report.OR.objpair[,7]),1,0))
pOR.objpair.inCI99<-mean(ifelse((report.OR.objpair[,4]>= report.OR.objpair[,8])& (report.OR.objpair[,4]<= report.OR.objpair[,9]),1,0))

##########################################################
## compute conditional and marginal probabilities
###########################################################

# disjunctive, conjunctive
if (maprule=="disj"|maprule=="conj"){
condprob.JKT<-array(rep(1,J*K*T),c(J,K,T))
for (t in 1:T){
for (f in 1:F){
condprob.JKT[,,t]<-condprob.JKT[,,t]*(1-sigma.best[,f,t]%o%rho.best[,f])}}
condprob.JKT<-1-condprob.JKT
}
condprob.JKT<-array(rep(0,J*K*T),c(J,K,T))
for (t in 1:T){
for (f in 1:F){
condprob.JKT[,,t]<-condprob.JKT[,,t]+(sigma.best[,f,t]%o%rho.best[,f])}}
condprob.JKT<-condprob.JKT/F
}

margprob.JK<-matrix(rep(1,J*K),nrow=J)
weight<-rep(1,J)%o%rep(1,K)%o%gamma.best
margprob.JK<-apply(weight*condprob.JKT,c(1,2),sum)
correl<-cor(c(apply(data,c(2,3),mean)),c(margprob.JK))
VAF<-correl^2

rownames(margprob.JK)<-objectlabels
colnames(margprob.JK)<-attributelabels
dimnames(condprob.JKT)[[1]]<-as.list(objectlabels)
dimnames(condprob.JKT)[[2]]<-as.list(attributelabels)
dimnames(condprob.JKT)[[3]]<-as.list(classlabels)

#########################################
## compute fit measures
########################################

deviance<--2*loglik.best
AIC<-deviance+2*Npar
BIC<-deviance+Npar*log(I)
fitmeasures<-matrix(c(loglik.best,logpost.best,deviance,AIC,BIC,correl,VAF,pOR.attpair.inCI95,pOR.objpair.inCI95,pOR.attpair.inCI99,pOR.objpair.inCI99),ncol=1)
rownames(fitmeasures)<-c("Log Likelihood","Log Posterior","Deviance","AIC","BIC","Correlation observed and expected frequencies J X K table","VAF observed frequencies J X K table",
"Proportion OR attribute pairs in 95% CI","Proportion OR object pairs in 95% CI","Proportion OR attribute pairs in 99% CI","Proportion OR object pairs in 99% CI")

tempcall<-match.call()
if (maprule=="disj") tempcall\$maprule<-"disj"
if (maprule=="conj") tempcall\$maprule<-"conj"

tempcall\$model<-1

LCplfm<-list(call=tempcall,logpost.runs=logpost.runs,best=best,
objpar=sigma.best,attpar=rho.best,sizepar=gamma.best,
SE.objpar=SEsigma.best,SE.attpar=SErho.best,SE.sizepar=SEgamma.best,
fitmeasures=fitmeasures,postprob=postprob.best,margprob.JK=margprob.JK,condprob.JKT=condprob.JKT,
report.OR.attpair=report.OR.attpair,report.OR.objpair=report.OR.objpair)
}
else if (maprule=="conj"){
LCplfm<-list(call=tempcall,logpost.runs=logpost.runs,best=best,
objpar=1-sigma.best,attpar=rho.best,sizepar=gamma.best,
SE.objpar=SEsigma.best,SE.attpar=SErho.best,SE.sizepar=SEgamma.best,
fitmeasures=fitmeasures,postprob=postprob.best,margprob.JK=1-margprob.JK,condprob.JKT=1-condprob.JKT,
report.OR.attpair=report.OR.attpair,report.OR.objpair=report.OR.objpair)

}
class(LCplfm)<-"LCplfm"
LCplfm

}## end model 1

## estimate model M_2
## LCPLFM with constant classification objects and heterogeneity attribute parameters

else