# Data
# ----
yt     = matrix(scan("pollutiondata-NO3.txt"),342,22,byrow=T)
subset = c(1,2,4,13,14,16,17,19,20)
y      = yt[,subset]
y      = scale(y)
n      = nrow(y)
p      = ncol(y)

par(mfrow=c(3,3))
for (i in 1:p)
  ts.plot(y[,i],ylab="",main=paste("Location: ",subset[i],sep=""))

# Classical factor analysis
# -------------------------
k   = 1
fac = factanal(y,k,scores="regression")

fac

fac$load

par(mfrow=c(1,2))
ts.plot(fac$scores[,1])
acf(fac$scores[,1],lag.max=50)


# Bayesian factor analyis
# -----------------------
source("factormodel-routines.R")

# Prior specification
# -------------------
m0    =  0.0
C0    =  1.0
c0    =  sqrt(C0)
iC0   =  1/C0
m0c0  =  m0/c0
m0C0  =  m0*iC0
E     = rep(1,p)
D     = 0.25
a     = (E/D)^2+2
b     = (a-1)*E

set.seed(23733)
M0 = 2000        # burn-in length
M  = 2000        # length of the chain, after burn-in)

# Inivial values
beta     = matrix(fac$load[1:p,1],p,1)
beta     = beta/beta[1,1]
sigma2   = fac$uniqueness

# Storage variables
betas    = matrix(0,M,p)
sigma2s  = matrix(0,M,p)
mfactor  = 0.0
VD1      = array(0,c(M,p,2))

# Algorithm
for (iter in 1:(M0+M)){
  f      = sampling.f(beta,sigma2)
  sigma2 = sampling.sigma2(f,beta)
  beta   = sampling.beta(f,sigma2)
  if (iter>M0){
    mfactor = mfactor + f
    sigma2s[iter-M0,] = sigma2
    betas[iter-M0,] = beta[,1]
    V = diag(beta%*%t(beta))
    O = V + sigma2
    VD1[iter-M0,,1] = V/O
    VD1[iter-M0,,2] = sigma2/O
  }
}
mfactor = mfactor/M

par(mfrow=c(1,1))
ts.plot(fac$scores[,1])
lines(mfactor,col=2)

par(mfrow=c(3,3))
for (i in 1:p){
  hist(betas[,i],prob=TRUE,main=paste("beta(",i,")",sep=""),xlab="")
  abline(v=fac$load[i,1],col=2,lwd=2)
}

par(mfrow=c(3,3))
for (i in 1:p){
  hist(sigma2s[,i],prob=TRUE,main=paste("sigma2(",i,")",sep=""),xlab="")
  abline(v=fac$uniq[i],col=2,lwd=2)
}

par(mfrow=c(3,3))
for (i in 1:p){
  hist(VD1[,i,1],prob=TRUE,main=paste("y(",i,")",sep=""),xlab="",xlim=c(0,1))
  abline(v=0.5,lty=2)
}

par(mfrow=c(1,1))
boxplot(VD1[,,1],ylim=c(0,1),outline=FALSE)
abline(h=0.5,lty=2)