Construction & inference (Time series) in C#

// --------------------------------------------------------------------------------------------------------------------
// <copyright file="DbnExample.cs" company="Bayes Server">
//   Copyright (C) Bayes Server.  All rights reserved.
// </copyright>
// --------------------------------------------------------------------------------------------------------------------

namespace BayesServer.HelpSamples
{
    using System;

    using BayesServer.Inference.RelevanceTree;

    public static class DbnExample
    {
        public static void Main()
        {
            // 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.

            var network = new Network("DBN");

            var cluster1 = new State("Cluster1");
            var cluster2 = new State("Cluster2");
            var cluster3 = new State("Cluster3");
            var varTransition = new Variable("Transition", cluster1, cluster2, cluster3);
            var nodeTransition = new Node(varTransition);

            // make the node temporal, so that it appears in each time slice
            nodeTransition.TemporalType = TemporalType.Temporal;  

            var varObs1 = new Variable("Obs1", VariableValueType.Continuous);
            var varObs2 = new Variable("Obs2", VariableValueType.Continuous);
            var varObs3 = new Variable("Obs3", VariableValueType.Continuous);
            var varObs4 = new Variable("Obs4", VariableValueType.Continuous);

            // observation node is a multi variable node, consisting of 4 continuous variables
            var nodeObservation = new Node("Observation", new Variable[] { varObs1, varObs2, varObs3, varObs4 });
            nodeObservation.TemporalType = TemporalType.Temporal;

            network.Nodes.Add(nodeTransition);
            network.Nodes.Add(nodeObservation);

            // link the transition node to the observation node within each time slice
            network.Links.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.Links.Add(new Link(nodeTransition, nodeTransition, 1));

            // 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

            var cluster1Time0 = new StateContext(cluster1, 0);
            var cluster2Time0 = new StateContext(cluster2, 0);
            var cluster3Time0 = new StateContext(cluster3, 0);

            var prior = nodeTransition.NewDistribution(0).Table;
            prior[cluster1Time0] = 0.2;
            prior[cluster2Time0] = 0.3;
            prior[cluster3Time0] = 0.5;

            // NewDistribution does not assign the new distribution, so it still must be assigned
            nodeTransition.Distribution = prior;

            // the second is specified for time >= 1
            var transition = nodeTransition.NewDistribution(1).Table;

            // 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

            var cluster1TimeM1 = new StateContext(cluster1, -1);
            var cluster2TimeM1 = new StateContext(cluster2, -1);
            var cluster3TimeM1 = new StateContext(cluster3, -1);

            transition[cluster1TimeM1, cluster1Time0] = 0.2;
            transition[cluster1TimeM1, cluster2Time0] = 0.3;
            transition[cluster1TimeM1, cluster3Time0] = 0.5;
            transition[cluster2TimeM1, cluster1Time0] = 0.4;
            transition[cluster2TimeM1, cluster2Time0] = 0.4;
            transition[cluster2TimeM1, cluster3Time0] = 0.2;
            transition[cluster3TimeM1, cluster1Time0] = 0.9;
            transition[cluster3TimeM1, cluster2Time0] = 0.09;
            transition[cluster3TimeM1, cluster3Time0] = 0.01;

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

            //new TableIterator(transition, new Variable[] { varTransition, varTransition }, new int?[] { -1, 0 }).CopyFrom(new double[]
            //    {
            //        0.2, 0.3, 0.5, 0.4, 0.4, 0.2, 0.9, 0.09, 0.01
            //    });

            nodeTransition.Distributions[1] = 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.
            var gaussian = (CLGaussian)nodeObservation.NewDistribution();

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

            var varObs1Time0 = new VariableContext(varObs1, 0, HeadTail.Head);
            var varObs2Time0 = new VariableContext(varObs2, 0, HeadTail.Head);
            var varObs3Time0 = new VariableContext(varObs3, 0, HeadTail.Head);
            var varObs4Time0 = new VariableContext(varObs4, 0, HeadTail.Head);

            gaussian.SetMean(varObs1Time0, 3.2, cluster1Time0);
            gaussian.SetMean(varObs2Time0, 2.4, cluster1Time0);
            gaussian.SetMean(varObs3Time0, -1.7, cluster1Time0);
            gaussian.SetMean(varObs4Time0, 6.2, cluster1Time0);

            gaussian.SetVariance(varObs1Time0, 2.3, cluster1Time0);
            gaussian.SetVariance(varObs2Time0, 2.1, cluster1Time0);
            gaussian.SetVariance(varObs3Time0, 3.2, cluster1Time0);
            gaussian.SetVariance(varObs4Time0, 1.4, cluster1Time0);

            gaussian.SetCovariance(varObs1Time0, varObs2Time0, -0.3, cluster1Time0);
            gaussian.SetCovariance(varObs1Time0, varObs3Time0, 0.5, cluster1Time0);
            gaussian.SetCovariance(varObs1Time0, varObs4Time0, 0.35, cluster1Time0);
            gaussian.SetCovariance(varObs2Time0, varObs3Time0, 0.12, cluster1Time0);
            gaussian.SetCovariance(varObs2Time0, varObs4Time0, 0.1, cluster1Time0);
            gaussian.SetCovariance(varObs3Time0, varObs4Time0, 0.23, cluster1Time0);

            // set the Gaussian parameters corresponding to the state "Cluster2" of variable "transition"
            gaussian.SetMean(varObs1Time0, 3.0, cluster2Time0);
            gaussian.SetMean(varObs2Time0, 2.8, cluster2Time0);
            gaussian.SetMean(varObs3Time0, -2.5, cluster2Time0);
            gaussian.SetMean(varObs4Time0, 6.9, cluster2Time0);

            gaussian.SetVariance(varObs1Time0, 2.1, cluster2Time0);
            gaussian.SetVariance(varObs2Time0, 2.2, cluster2Time0);
            gaussian.SetVariance(varObs3Time0, 3.3, cluster2Time0);
            gaussian.SetVariance(varObs4Time0, 1.5, cluster2Time0);

            gaussian.SetCovariance(varObs1Time0, varObs2Time0, - 0.4, cluster2Time0);
            gaussian.SetCovariance(varObs1Time0, varObs3Time0, 0.5, cluster2Time0);
            gaussian.SetCovariance(varObs1Time0, varObs4Time0, 0.45, cluster2Time0);
            gaussian.SetCovariance(varObs2Time0, varObs3Time0, 0.22, cluster2Time0);
            gaussian.SetCovariance(varObs2Time0, varObs4Time0, 0.15, cluster2Time0);
            gaussian.SetCovariance(varObs3Time0, varObs4Time0, 0.24, cluster2Time0);

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

            gaussian.SetMean(varObs1Time0, 3.8, cluster3Time0);
            gaussian.SetMean(varObs2Time0, 2.0, cluster3Time0);
            gaussian.SetMean(varObs3Time0, -1.9, cluster3Time0);
            gaussian.SetMean(varObs4Time0, 6.25, cluster3Time0);

            gaussian.SetVariance(varObs1Time0, 2.34, cluster3Time0);
            gaussian.SetVariance(varObs2Time0, 2.11, cluster3Time0);
            gaussian.SetVariance(varObs3Time0, 3.22, cluster3Time0);
            gaussian.SetVariance(varObs4Time0, 1.43, cluster3Time0);

            gaussian.SetCovariance(varObs1Time0, varObs2Time0, -0.31, cluster3Time0);
            gaussian.SetCovariance(varObs1Time0, varObs3Time0, 0.52, cluster3Time0);
            gaussian.SetCovariance(varObs1Time0, varObs4Time0, 0.353, cluster3Time0);
            gaussian.SetCovariance(varObs2Time0, varObs3Time0, 0.124, cluster3Time0);
            gaussian.SetCovariance(varObs2Time0, varObs4Time0, 0.15, cluster3Time0);
            gaussian.SetCovariance(varObs3Time0, varObs4Time0, 0.236, cluster3Time0);

            nodeObservation.Distribution = 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

            var inference = new RelevanceTreeInference(network);
            var queryOptions = new RelevanceTreeQueryOptions();
            var queryOutput = new RelevanceTreeQueryOutput();

            // set some temporal evidence

            inference.Evidence.Set(varObs1, new double?[] { 2.2, 2.4, 2.6, 2.9}, 0, 0, 4);
            inference.Evidence.Set(varObs2, new double?[] { null, 4.0, 4.1, 4.88}, 0, 0, 4);
            inference.Evidence.Set(varObs3, new double?[] {-2.5, -2.3, null, -4.0 }, 0, 0, 4);
            inference.Evidence.Set(varObs4, new double?[] { 4.0, 6.5, 4.9, 4.4}, 0, 0, 4);

            queryOptions.LogLikelihood = true; // only ask for this if you really need it

            // predict the observation variables one time step in the future
            var predictTime = 4;

            var gaussianFuture = new CLGaussian[nodeObservation.Variables.Count];

            for(int i = 0; i < gaussianFuture.Length; i++)
            {
                gaussianFuture[i] = new CLGaussian(nodeObservation.Variables[i], predictTime);
                inference.QueryDistributions.Add(gaussianFuture[i]);
            }

            // we will also demonstrate querying a joint distribution

            var jointFuture = new CLGaussian(new Variable[] { varObs1, varObs2 }, predictTime);
            inference.QueryDistributions.Add(jointFuture);


            inference.Query(queryOptions, queryOutput); // note that this can raise an exception (see help for details)

            Console.WriteLine("LogLikelihood: " + queryOutput.LogLikelihood.Value);
            Console.WriteLine();

            for (int h = 0; h < gaussianFuture.Length; h++)
            {
                var variableH = nodeObservation.Variables[h];
                Console.WriteLine("Mean({0}(t=4)|evidence)={1}", variableH.Name, gaussianFuture[h].GetMean(variableH, predictTime));
            }

            Console.WriteLine();
            Console.WriteLine("{0},{1}|evidence:", varObs1.Name, varObs2.Name);
            Console.WriteLine(jointFuture.GetMean(varObs1, predictTime) + "\t" + jointFuture.GetMean(varObs2, predictTime));
            Console.WriteLine(jointFuture.GetVariance(varObs1, predictTime) + "\t" + jointFuture.GetCovariance(varObs1, predictTime, varObs2, predictTime));
            Console.WriteLine(jointFuture.GetCovariance(varObs2, predictTime, varObs1, predictTime) + "\t" + 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

        }
    }
}