library(lda) Y<-as.matrix(read.table("CaltechY.dat")) # replace the diagonals with 1s diag(Y)<-1 #net<-network(Y,directed=FALSE) link.predict<- function(fit, Y){ N<- nrow(fit\$net.assignments.left) theta<- thetahat(fit,Y) EY<- matrix(0,N,N) for(i in 1:N){ for(j in (1:N)){ if(i==j) EY[i,j]<- 0 else EY[i,j]<- theta[fit\$net.assignments.right[i,j]+1 , fit\$net.assignments.left[i,j]+1] } } return(EY) } thetahat<- function(fit, Y){ G<- nrow(fit\$blocks.pos) theta.pos<- theta.neg<- matrix(0, G, G) for(g1 in 1:G){ for(g2 in 1:G){ theta.pos[g1, g2]<- sum((fit\$net.assignments.left==(g1-1))*(fit\$net.assignments.right==(g2-1))*(Y)) theta.neg[g1, g2]<- sum((fit\$net.assignments.left==(g1-1))*(fit\$net.assignments.right==(g2-1))*(1-Y)) } } theta.pos/(theta.pos + theta.neg) } fit2<- mmsb.collapsed.gibbs.sampler(network=Y,K=2, num.iterations=5100, alpha=1/2, beta.prior=list(matrix(1,2,2), matrix(1,2,2)), burnin=as.integer(50000)) fit3<- mmsb.collapsed.gibbs.sampler(network=Y,K=3, num.iterations=5100, alpha=1/3, beta.prior=list(matrix(1,3,3), matrix(1,3,3)), burnin=as.integer(50000)) fit4<- mmsb.collapsed.gibbs.sampler(network=Y,K=4, num.iterations=5100, alpha=1/4, beta.prior=list(matrix(1,4,4), matrix(1,4,4)), burnin=as.integer(50000)) fit5<- mmsb.collapsed.gibbs.sampler(network=Y,K=5, num.iterations=5100, alpha=1/5, beta.prior=list(matrix(1,5,5), matrix(1,5,5)), burnin=as.integer(50000)) ## Deliberately Overfit fit9<- mmsb.collapsed.gibbs.sampler(network=Y,K=9, num.iterations=5100, alpha=1/9, beta.prior=list(matrix(1,9,9), matrix(1,9,9)), burnin=as.integer(50000)) fit10<- mmsb.collapsed.gibbs.sampler(network=Y,K=10, num.iterations=5100, alpha=1/10, beta.prior=list(matrix(1,10,10), matrix(1,10,10)), burnin=as.integer(50000)) fit11<- mmsb.collapsed.gibbs.sampler(network=Y,K=11, num.iterations=5100, alpha=1/11, beta.prior=list(matrix(1,11,11), matrix(1,11,11)), burnin=as.integer(50000)) l2<- link.predict(fit2, Y) l3<- link.predict(fit3, Y) l4<- link.predict(fit4, Y) l5<- link.predict(fit5, Y) l9<- link.predict(fit9, Y) l10<- link.predict(fit10, Y) l11<- link.predict(fit11, Y)