module HasGP.Classification.EP.ClassificationEP
(
EPValue(eValue,siteState,count),
EPConvergenceTest,
EPSiteState,
EPState,
SiteOrder,
generateRandomSiteOrder,
generateFixedSiteOrder,
gpClassifierEPEvidence,
gpClassifierEPLearn,
gpClassifierEPPredict,
gpClassifierEPLogEvidence,
gpClassifierEPLogEvidenceList,
gpClassifierEPLogEvidenceVec
) where
import Numeric.LinearAlgebra
import Control.Monad.State
import System.Random
import HasGP.Types.MainTypes
import HasGP.Support.MatrixFunction
import HasGP.Support.Linear
import HasGP.Support.Functions
import HasGP.Support.Solve
import HasGP.Support.Iterate
import HasGP.Support.Random
import HasGP.Covariance.Basic
import HasGP.Likelihood.Basic
data EPValue = EPValue {
eValue::Double,
siteState::EPSiteState,
count::Int
}
type EPConvergenceTest = (EPValue -> EPValue -> Bool)
data EPSiteState = EPSiteState {
var::DMatrix,
tauTilde::DVector,
mu::DVector,
nuTilde::DVector,
tauMinus::DVector,
muMinus::DVector
}
type EPState = (EPSiteState, StdGen, Int)
type SiteOrder = State EPState [Int]
generateInitialSiteState :: CovarianceMatrix
-> Int
-> EPSiteState
generateInitialSiteState k n = EPSiteState k z1 z2 z3 z4 z5
where
z1 = constant 0.0 n
z2 = constant 0.0 n
z3 = constant 0.0 n
z4 = constant 0.0 n
z5 = constant 0.0 n
cavityParameters :: Double
-> Double
-> Double
-> Double
-> (Double,Double)
cavityParameters varI tauTildeI muI nuTildeI = (tauMinusI, nuMinusI)
where
varIInv = 1 / varI
tauMinusI = varIInv tauTildeI
nuMinusI = (varIInv * muI) nuTildeI
marginalMoments :: Double
-> Double
-> Double
-> (Double,Double)
marginalMoments muMinusI tI varMinusI = (muHatI, varHatI)
where
zI = (tI * muMinusI) / (sqrt (1 + varMinusI))
nopz = nOverPhi zI
muHatI = muMinusI + (((tI * varMinusI) / (sqrt (1 + varMinusI))) * nopz)
varHatI = varMinusI
((((square varMinusI) / (1 + varMinusI)) * nopz) * (zI + nopz))
siteParameters :: Double
-> Double
-> Double
-> Double
-> Double
-> (Double,Double,Double)
siteParameters tauTildeI varHatI tauMinusI muHatI nuMinusI =
(deltaTauTilde, tauTildeINew, nuTildeINew)
where
deltaTauTilde = (1 / varHatI) tauMinusI tauTildeI
tauTildeINew = tauTildeI + deltaTauTilde
nuTildeINew = ((1 / varHatI) * muHatI) nuMinusI
updateOneSite :: Targets
-> Int
-> Int
-> EPSiteState
-> EPSiteState
updateOneSite t n i (EPSiteState var tauTilde mu nuTilde oldTauMinus
oldMuMinus) =
EPSiteState newVar newTauTilde newMu newNuTilde newTauMinus newMuMinus
where
i2 = i1
tauTildeI = (tauTilde @> i2)
varII = var @@> (i2,i2)
(tauMinusI, nuMinusI) =
cavityParameters varII tauTildeI (mu @> i2) (nuTilde @> i2)
varMinusI = (1 / tauMinusI)
muMinusI = (nuMinusI / tauMinusI)
(muHatI, varHatI) = marginalMoments muMinusI (t @> i2) varMinusI
(deltaTauTilde, tauTildeINew, nuTildeINew) =
siteParameters tauTildeI varHatI tauMinusI muHatI nuMinusI
newTauTilde = replaceInVector tauTilde (i2+1) tauTildeINew
newNuTilde = replaceInVector nuTilde (i2+1) nuTildeINew
newTauMinus = replaceInVector oldTauMinus (i2+1) tauMinusI
newMuMinus = replaceInVector oldMuMinus (i2+1) muMinusI
sI = subMatrix (0,i2) (n,1) var
newVar = var (scale (deltaTauTilde / (1 + (deltaTauTilde * varII)))
(sI <> (trans sI)))
newMu = flatten $ (newVar <> (asColumn newNuTilde))
randomPermutation :: StdGen
-> Int
-> (StdGen, [Int])
randomPermutation g n = rP g (n1) [1..n] []
where
rP g' n' [] result = (g', result)
rP g' n' [x] result = (g', x:result)
rP g' n' xs result = rP newG (n' 1) newXs (m:result)
where
(r, newG) = randomR (0, n') g'
m = xs !! r
newXs = filter (\x -> x /= m) xs
generateRandomSiteOrder :: SiteOrder
generateRandomSiteOrder = do
(state, g, n) <- get
let (newG, p) = randomPermutation g (dim $ tauTilde state)
put (state, newG, n)
return p
generateFixedSiteOrder :: SiteOrder
generateFixedSiteOrder = do
(state, g, n) <- get
return [1..(dim $ tauTilde state)]
updateAllSites :: Targets
-> Int
-> [Int]
-> EPSiteState
-> EPSiteState
updateAllSites t n [] state = state
updateAllSites t n (s:ss) state = updateAllSites t n ss newState
where
newState = updateOneSite t n s state
recomputeApproximation :: CovarianceMatrix
-> Int
-> DVector
-> DVector
-> (DMatrix, DMatrix, DVector)
recomputeApproximation k n tauTilde nuTilde =
(l, finalVar, finalMu)
where
rootSTilde = mapVector sqrt tauTilde
l = trans $ chol ((ident n) + (abaDiagDiag rootSTilde k))
v = generalSolve lowerSolve l (preMultiply rootSTilde k)
finalVar = k ((trans v) <> v)
finalMu = flatten $ finalVar <> (asColumn nuTilde)
gpClassifierEPEvidence :: CovarianceMatrix
-> Targets
-> DMatrix
-> EPSiteState
-> Double
gpClassifierEPEvidence k t l state =
(terms1and4 + term3 + terms2and5)
where
tT = tauTilde state
nT = nuTilde state
tM = tauMinus state
mM = muMinus state
oneOverTauMinus = mapVector (1/) $ tM
sumLog = sum $ toList $ mapVector log $ takeDiag l
terms1and4 =
(0.5 * (sum $ toList $ mapVector (log . (1+))
(zipVectorWith (*) tT oneOverTauMinus))) sumLog
term3 =
sum $ toList $
mapVector logPhi (zipVectorWith (/)
(zipVectorWith (*) t mM)
(mapVector (sqrt . (1+)) oneOverTauMinus))
rootSTilde = diag $ mapVector sqrt tT
vM = generalSolve lowerSolve l
(preMultiply (mapVector sqrt tT) k)
invTPlusSTilde =
diag $ mapVector (1/) $ tT + tM
part1 =
0.5 * ((asRow nT) <>
(k ((trans vM) <> vM) invTPlusSTilde) <>
(asColumn nT))
part2 =
0.5 * ((asRow mM) <> (diag tM) <> invTPlusSTilde <>
((asColumn (zipVectorWith (*) tT mM))
(asColumn $ (scale 2 nT))))
terms2and5 = (part1 + part2) @@> (0,0)
doOneUpdate :: CovarianceMatrix
-> Targets
-> SiteOrder
-> State EPState EPValue
doOneUpdate k t siteOrder = do
order <- siteOrder
(state,g,i) <- get
let sites = dim $ tauTilde state
let state' = updateAllSites t sites order state
let finalTT = tauTilde state'
let finalNT = nuTilde state'
let finalTM = tauMinus state'
let finalMM = muMinus state'
let (l, finalVar, finalMu) = recomputeApproximation k sites finalTT finalNT
let state'' = EPSiteState finalVar finalTT finalMu finalNT finalTM finalMM
put (state'', g, i+1)
let logML = gpClassifierEPEvidence k t l state''
return $ EPValue logML state'' (i+1)
gpClassifierEPLearn :: CovarianceMatrix
-> Targets
-> SiteOrder
-> EPConvergenceTest
-> (EPValue, EPState)
gpClassifierEPLearn k t siteOrder converged =
runState (iterateToConvergence'' doOnce converged) start
where
doOnce = doOneUpdate k t siteOrder
start = ((generateInitialSiteState k (dim t)), mkStdGen 0, 0)
gpClassifierEPPredict :: (CovarianceFunction c) => EPSiteState
-> Inputs
-> Targets
-> CovarianceMatrix
-> c
-> Inputs
-> DVector
gpClassifierEPPredict state i t k c xStars
= fromList $ map phiIntegral (zipWith (/) fStar (map (sqrt . (1+)) vfStar))
where
nT = nuTilde state
tT = tauTilde state
d = dim t
rootSTildeV = mapVector sqrt tT
rootSTilde = diag rootSTildeV
l = trans $ chol ((ident d) + (abaDiagDiag rootSTildeV k))
z = rootSTilde
<> (asColumn
(upperSolve (trans l)
(lowerSolve l (flatten $
(rootSTilde <> k <> (asColumn nT))))))
xStarsRows = toRows xStars
covarianceWithTestInputs =
fromRows [covarianceWithPoint c i xStar | xStar <- xStarsRows]
fStar =
toList $ flatten $
(covarianceWithTestInputs <> ((asColumn nT) z))
v = [lowerSolve l (rootSTilde <> kxStar) |
kxStar <- (toColumns $ trans $ covarianceWithTestInputs)]
vTv = zipWith (<.>) v v
kxStarxStar = zipWith (covariance c) xStarsRows xStarsRows
vfStar = zipWith () kxStarxStar vTv
gpClassifierEPLogEvidence :: (CovarianceFunction c) => c
-> Inputs
-> Targets
-> SiteOrder
-> EPConvergenceTest
-> (Double, DVector)
gpClassifierEPLogEvidence c i t siteOrder converged
= (logEvidence, (zipVectorWith (*) (trueHyper c) (fromList dLogEvidence)))
where
d = dim t
k = covarianceMatrix c i
(value, s) = gpClassifierEPLearn k t siteOrder converged
nT = nuTilde $ siteState value
tT = tauTilde $ siteState value
logEvidence = eValue value
rootSTildeV = mapVector sqrt tT
rootSTilde = diag rootSTildeV
l = trans $ chol ((ident d) + (abaDiagDiag rootSTildeV k))
b = (asColumn nT)
(rootSTilde <>
(asColumn $ upperSolve (trans l)
(lowerSolve l (flatten $
rootSTilde <> k <> (asColumn nT)))))
r = (b <> (trans b))
(rootSTilde <>
(generalSolve upperSolve (trans l)
(generalSolve lowerSolve l rootSTilde)))
dLogEvidence = map ((0.5 *) . sum . toList . (abDiagOnly r))
(makeMatricesFromPairs (dCovarianceDParameters c) i)
gpClassifierEPLogEvidenceList :: (CovarianceFunction c) => Inputs
-> Targets
-> c
-> SiteOrder
-> EPConvergenceTest
-> [Double]
-> (Double, DVector)
gpClassifierEPLogEvidenceList i t cov siteOrder converged hyper =
gpClassifierEPLogEvidence cov2 i t siteOrder converged
where
cov2 = makeCovarianceFromList cov hyper
gpClassifierEPLogEvidenceVec :: (CovarianceFunction c) => Inputs
-> Targets
-> c
-> SiteOrder
-> EPConvergenceTest
-> DVector
-> (Double, DVector)
gpClassifierEPLogEvidenceVec i t cov siteOrder converged hyper =
gpClassifierEPLogEvidence cov2 i t siteOrder converged
where
cov2 = makeCovarianceFromList cov (toList hyper)