Construction & inference (Time series) in R


source("bayesserver.R")


# In this example we programatically create a Dynamic Bayesian network (time series).
# Note that you can automatically define nodes from data using
# classes in BayesServer.Data.Discovery,
# and you can automatically learn the parameters using classes in
# BayesServer.Learning.Parameters,
# however here we build a Bayesian network from scratch.

network <- new(Network, "DBN")

cluster1 <- new(State, "Cluster1")
cluster2 <- new(State, "Cluster2")
cluster3 <- new(State, "Cluster3")
varTransition <- new(Variable, "Transition", toStateArray(cluster1, cluster2, cluster3))
nodeTransition <- new(Node, varTransition)

# make the node temporal, so that it appears in each time slice
nodeTransition$setTemporalType(TemporalType$TEMPORAL)

varObs1 <- new(Variable, "Obs1", VariableValueType$CONTINUOUS)
varObs2 <- new(Variable, "Obs2", VariableValueType$CONTINUOUS)
varObs3 <- new(Variable, "Obs3", VariableValueType$CONTINUOUS)
varObs4 <- new(Variable, "Obs4", VariableValueType$CONTINUOUS)

# observation node is a multi variable node, consisting of 4 continuous variables
nodeObservation <- new(Node, "Observation", toVariableArray(varObs1, varObs2, varObs3, varObs4))
nodeObservation$setTemporalType(TemporalType$TEMPORAL)

network$getNodes()$add(nodeTransition)
network$getNodes()$add(nodeObservation)

# link the transition node to the observation node within each time slice
network$getLinks()$add(new(Link, nodeTransition, nodeObservation))

# add a temporal link of order 1.  This links the transition node to itself in the next time slice
network$getLinks()$add(new(Link, nodeTransition, nodeTransition, 1L))

# at this point the structural specification is complete

# now complete the distributions

# because the transition node has an incoming temporal link of order 1 (from itself), we must specify
# two distributions, the first of which is specified for time = 0

cluster1Time0 <- new(StateContext, cluster1, new(Integer, "0"))
cluster2Time0 <- new(StateContext, cluster2, new(Integer, "0"))
cluster3Time0 <- new(StateContext, cluster3, new(Integer, "0"))

prior <- nodeTransition$newDistribution(0L)$getTable()
prior$set(0.2, toStateContextArray(cluster1Time0))
prior$set(0.3, toStateContextArray(cluster2Time0))
prior$set(0.5, toStateContextArray(cluster3Time0))

# NewDistribution does not assign the new distribution, so it still must be assigned
nodeTransition$setDistribution(prior)

# the second is specified for time >= 1
transition <- nodeTransition$newDistribution(1L)$getTable()

# when specifying temporal distributions, variables which belong to temporal nodes must have times associated
# NOTE: Each time is specified relative to the current point in time which is defined as zero,
# therefore the time for variables at the previous time step is -1

cluster1TimeM1 <- new(StateContext, cluster1, new(Integer, "-1"))
cluster2TimeM1 <- new(StateContext, cluster2, new(Integer, "-1"))
cluster3TimeM1 <- new(StateContext, cluster3, new(Integer, "-1"))

transition$set(0.2, toStateContextArray(cluster1TimeM1, cluster1Time0))
transition$set(0.3, toStateContextArray(cluster1TimeM1, cluster2Time0))
transition$set(0.5, toStateContextArray(cluster1TimeM1, cluster3Time0))
transition$set(0.4, toStateContextArray(cluster2TimeM1, cluster1Time0))
transition$set(0.4, toStateContextArray(cluster2TimeM1, cluster2Time0))
transition$set(0.2, toStateContextArray(cluster2TimeM1, cluster3Time0))
transition$set(0.9, toStateContextArray(cluster3TimeM1, cluster1Time0))
transition$set(0.09, toStateContextArray(cluster3TimeM1, cluster2Time0))
transition$set(0.01, toStateContextArray(cluster3TimeM1, cluster3Time0))

# an alternative would be to set values using TableIterator.CopyFrom

nodeTransition$getDistributions()$set(1L, transition)

# Node observation does not have any incoming temporal links, so
# only requires a distribution specified at time >=0
# Calling NewDistribution without specifying a time assumes time zero.
gaussian <- nodeObservation$newDistribution()

# set the Gaussian parameters corresponding to the state "Cluster1" of variable "transition"

varObs1Time0 = new(VariableContext, varObs1, new(Integer, "0"), HeadTail$HEAD)
varObs2Time0 = new(VariableContext, varObs2, new(Integer, "0"), HeadTail$HEAD)
varObs3Time0 = new(VariableContext, varObs3, new(Integer, "0"), HeadTail$HEAD)
varObs4Time0 = new(VariableContext, varObs4, new(Integer, "0"), HeadTail$HEAD)

gaussian$setMean(varObs1Time0, 3.2, toStateContextArray(cluster1Time0))
gaussian$setMean(varObs2Time0, 2.4, toStateContextArray(cluster1Time0))
gaussian$setMean(varObs3Time0, -1.7, toStateContextArray(cluster1Time0))
gaussian$setMean(varObs4Time0, 6.2, toStateContextArray(cluster1Time0))

gaussian$setVariance(varObs1Time0, 2.3, toStateContextArray(cluster1Time0))
gaussian$setVariance(varObs2Time0, 2.1, toStateContextArray(cluster1Time0))
gaussian$setVariance(varObs3Time0, 3.2, toStateContextArray(cluster1Time0))
gaussian$setVariance(varObs4Time0, 1.4, toStateContextArray(cluster1Time0))

gaussian$setCovariance(varObs1Time0, varObs2Time0, -0.3, toStateContextArray(cluster1Time0))
gaussian$setCovariance(varObs1Time0, varObs3Time0, 0.5, toStateContextArray(cluster1Time0))
gaussian$setCovariance(varObs1Time0, varObs4Time0, 0.35, toStateContextArray(cluster1Time0))
gaussian$setCovariance(varObs2Time0, varObs3Time0, 0.12, toStateContextArray(cluster1Time0))
gaussian$setCovariance(varObs2Time0, varObs4Time0, 0.1, toStateContextArray(cluster1Time0))
gaussian$setCovariance(varObs3Time0, varObs4Time0, 0.23, toStateContextArray(cluster1Time0))

# set the Gaussian parameters corresponding to the state "Cluster2" of variable "transition"
gaussian$setMean(varObs1Time0, 3.0, toStateContextArray(cluster2Time0))
gaussian$setMean(varObs2Time0, 2.8, toStateContextArray(cluster2Time0))
gaussian$setMean(varObs3Time0, -2.5, toStateContextArray(cluster2Time0))
gaussian$setMean(varObs4Time0, 6.9, toStateContextArray(cluster2Time0))

gaussian$setVariance(varObs1Time0, 2.1, toStateContextArray(cluster2Time0))
gaussian$setVariance(varObs2Time0, 2.2, toStateContextArray(cluster2Time0))
gaussian$setVariance(varObs3Time0, 3.3, toStateContextArray(cluster2Time0))
gaussian$setVariance(varObs4Time0, 1.5, toStateContextArray(cluster2Time0))

gaussian$setCovariance(varObs1Time0, varObs2Time0, -0.4, toStateContextArray(cluster2Time0))
gaussian$setCovariance(varObs1Time0, varObs3Time0, 0.5, toStateContextArray(cluster2Time0))
gaussian$setCovariance(varObs1Time0, varObs4Time0, 0.45, toStateContextArray(cluster2Time0))
gaussian$setCovariance(varObs2Time0, varObs3Time0, 0.22, toStateContextArray(cluster2Time0))
gaussian$setCovariance(varObs2Time0, varObs4Time0, 0.15, toStateContextArray(cluster2Time0))
gaussian$setCovariance(varObs3Time0, varObs4Time0, 0.24, toStateContextArray(cluster2Time0))

# set the Gaussian parameters corresponding to the state "Cluster3" of variable "transition"

gaussian$setMean(varObs1Time0, 3.8, toStateContextArray(cluster3Time0))
gaussian$setMean(varObs2Time0, 2.0, toStateContextArray(cluster3Time0))
gaussian$setMean(varObs3Time0, -1.9, toStateContextArray(cluster3Time0))
gaussian$setMean(varObs4Time0, 6.25, toStateContextArray(cluster3Time0))

gaussian$setVariance(varObs1Time0, 2.34, toStateContextArray(cluster3Time0))
gaussian$setVariance(varObs2Time0, 2.11, toStateContextArray(cluster3Time0))
gaussian$setVariance(varObs3Time0, 3.22, toStateContextArray(cluster3Time0))
gaussian$setVariance(varObs4Time0, 1.43, toStateContextArray(cluster3Time0))

gaussian$setCovariance(varObs1Time0, varObs2Time0, -0.31, toStateContextArray(cluster3Time0))
gaussian$setCovariance(varObs1Time0, varObs3Time0, 0.52, toStateContextArray(cluster3Time0))
gaussian$setCovariance(varObs1Time0, varObs4Time0, 0.353, toStateContextArray(cluster3Time0))
gaussian$setCovariance(varObs2Time0, varObs3Time0, 0.124, toStateContextArray(cluster3Time0))
gaussian$setCovariance(varObs2Time0, varObs4Time0, 0.15, toStateContextArray(cluster3Time0))
gaussian$setCovariance(varObs3Time0, varObs4Time0, 0.236, toStateContextArray(cluster3Time0))

nodeObservation$setDistribution(gaussian)

# optional check to validate network
network$validate(new(ValidationOptions))


# at this point the network has been fully specified

# we will now perform some queries on the network

factory <- new(RelevanceTreeInferenceFactory)
inference <- factory$createInferenceEngine(network)
queryOptions <- factory$createQueryOptions()
queryOutput <- factory$createQueryOutput()

# set some temporal evidence

inference$getEvidence()$set(varObs1, toDoubleArray(2.2, 2.4, 2.6, 2.9), 0L, 0L, 4L)
inference$getEvidence()$set(varObs2, toDoubleArray(NA, 4.0, 4.1, 4.88), 0L, 0L, 4L)
inference$getEvidence()$set(varObs3, toDoubleArray(-2.5, -2.3, NA, -4.0), 0L, 0L, 4L)
inference$getEvidence()$set(varObs4, toDoubleArray(4.0, 6.5, 4.9, 4.4), 0L, 0L, 4L)

queryOptions$setLogLikelihood(TRUE) # only ask for this if you really need it

# predict the observation variables one time step in the future
predictTime = new(Integer, 4L)

gaussianFuture <- lapply(nodeObservation$getVariables(), function(v){
  
  g <- new(CLGaussian, v, predictTime)
  inference$getQueryDistributions()$add(g)
  return(g)
  
})

# we will also demonstrate querying a joint distribution

jointFuture <- new(CLGaussian, toList(varObs1, varObs2), predictTime)
inference$getQueryDistributions()$add(jointFuture)


inference$query(queryOptions, queryOutput) # note that this can raise an exception (see help for details)

print(sprintf("LogLikelihood: %f", queryOutput$getLogLikelihood()))

for (h in 1:(length(gaussianFuture))) {
  variableH <- nodeObservation$getVariables()$get(h-1L)
  mean <- gaussianFuture[[h]]$getMean(variableH, predictTime)
  print(sprintf("P(%s(t=4)|evidence)=%s", variableH$getName(), mean))
}

print(sprintf("P(%s,%s|evidence)=", varObs1$getName(), varObs2$getName()))
print(sprintf("%f  %f", jointFuture$getMean(varObs1, predictTime), jointFuture$getMean(varObs2, predictTime)))
print(sprintf("%f  %f", jointFuture$getVariance(varObs1, predictTime), jointFuture$getCovariance(varObs1, predictTime, varObs2, predictTime)))
print(sprintf("%f  %f", jointFuture$getCovariance(varObs2, predictTime, varObs1, predictTime), jointFuture$getVariance(varObs2, predictTime)))

# Expected output...

# LogLikelihood: -26.3688322999762

# P(Obs1(t=4)|evidence)=3.33914912825023
# P(Obs2(t=4)|evidence)=2.38039739886759
# P(Obs3(t=4)|evidence)=-1.98416436694525
# P(Obs4(t=4)|evidence)=6.40822262492584

# P(Obs1,Obs2|evidence)=
# 3.33914912825023        2.38039739886759
# 2.36608725717058        -0.427500059391733
# -0.427500059391733      2.22592296205311