SMC Procedure

This example illustrates how to estimate a dynamic Bayesian network (DBN). A heating, ventilation, and air conditioning (HVAC) system is used as a prototype network of DBN.

In this HVAC system, there are four state variables:

There are three observation variables:

The HVAC system can be represented by the following state space model (or the DBN model), where the subscript t denotes the time index:

Here denotes the Bernoulli distribution whose probability of success is p (success: 1; failure: 0), and denotes the normal distribution, whose mean is and variance is . Figure 34 shows the relationships within the HVAC system.

There are 16 parameters. The initial values, bounds, and actual values that the simulation data set uses for the parameters are listed in Table 23. Each of the variance parameters technically has an upper bound of , but they are constrained to have an upper bound of 0.99. This represents prior information that the variances cannot be larger than 0.99, which makes sense for this example because it is known that the true values of these parameters are all less than 0.99. Also, adding these constraints makes the model significantly easier to fit.

The following statements simulate 288 observations on the basis of the given HVAC system:

data hvacSystem;
   call streaminit(123);
   OSlag = 0;
   ISlag = 0;
   OTlag = 55;
   ITlag = 70;
   OTReadlag = 58;
   ITReadlag = 67;
   HVAClag = 3;
   pi=constant("pi");
   do t = 1 to 288*2;
      time = t;
      u = rand("uniform");
      if OSlag = 0 then OS=ifn(u<0.85,0,1);
      else OS=ifn(u<0.15,0,1);
      u = rand("uniform");
      if ISlag = 0 then IS=ifn(u<0.85,0,1);
      else IS=ifn(u<0.15,0,1);
      OT = OTlag+45*pi/287*cos(2*pi/287*(t-2)-0.5*pi)+ sqrt(0.3)*rand("normal");
      if HVAClag = 1 then IT = 0.9*ITlag + 0.2*OTlag + sqrt(0.85)*rand("normal");
      if HVAClag = 2 then IT = 0.7*ITlag + 0.2*OTlag + sqrt(0.85)*rand("normal");
      if HVAClag = 3 then IT = 0.8*ITlag + 0.2*OTlag + sqrt(0.85)*rand("normal");
      if OS = 0 then OTRead = OTReadlag + sqrt(0.4)*rand("normal");
      else OTRead = OT + sqrt(0.2)*rand("normal");
      if IS = 0 then ITRead = ITReadlag + sqrt(0.4)*rand("normal");
      else ITRead = IT + sqrt(0.2)*rand("normal");
      if IT < 68 then HVAC = 1;
      else HVAC=ifn(IT>74,2,3);
      output;
      OTReadlag = OTRead;
      ITReadlag = ITRead;
      OSlag = OS;
      ISlag = IS;
      OTlag = OT;
      ITlag = IT;
      HVAClag = HVAC;
   end;
run;

The following DATA step loads the data set hvacSystem into a data table in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.

data mylib.hvacSystem;
   set hvacSystem;
run;

In the SMC procedure, two categorical variables must be positioned before the two noncategorical ones in the state variable vector . For the four state variables, denotes OutSensorStatus; denotes InSensorStatus; denotes OutTemp; and denotes InTemp.

For the two dependent variables, OutTempRead is denoted by and InTempRead is denoted by . The independent variable, HVAC, is denoted by .

The following statements use the SUBMITMODEL statement in the SMC procedure to define the HVAC system and then use the PMCMC method to estimate the 16 parameters for the simulated data set. The median of the posterior parameter samples will be used as the parameter in the smoothing process. Moreover, in the SUBMITMODEL statement, is represented as for ; is represented by for ; is represented by for ; and is represented by . You must specify the dependent-variables in the SUBMITMODEL statement in the same order as . This example can take around 14 hours or longer to run, depending on which hardware you run it on. The TRANSFORM option for the PMCMC method in the LEARN statement significantly improves the efficiency of the sampler for this model, especially when you combine it with the constraints on the variance parameters.

proc smc data=mylib.hvacSystem  nthreads=32 seed=1;
   id time = t;
   learn nparticle=8000 algorithm=BF method=PMCMC(nsample=40000
        samsum samset nbi=20000 thin=5 seed=123
        outpost = mylib.outle sampler = rwm(transform ntu=700 ntunestage=48)
        statistics = all diagnostics(settings)=all) ;
   smooth nparticle = 20480 alg=bf parm=median method=fixedlag(30) out=mylib.sm;
   submitmodel OTRead ITRead = HVAC / nstate(2,2)=4 nParm = 16 name=dbn
      xlag = (state=1) ylag=(model=1) zlag=(state=1);
      subroutine dbn_sampprio(theta[*]);
         outargs theta;
         theta[1] = 0.65;
         theta[2] = 0.3;
         theta[3] = 0.65;
         theta[4] = 0.3;
         theta[5] = 0.65;
         theta[6] = 0.3;
         theta[7] = 0.7;
         theta[8] = 0.55;
         theta[9] = 0.3;
         theta[10] = 0.7;
         theta[11] = 0.6;
         theta[12] = 0.3;
         theta[13] = 0.7;
         theta[14] = 0.6;
         theta[15] = 0.6;
         theta[16] = 0.5;
      endsub;

      subroutine dbn_lowerbounds(theta[*]);
         outargs theta;
         theta[1] = 0.01;
         theta[2] = 0.01;
         theta[3] = 0.01;
         theta[4] = 0.01;
         theta[5] = 0.01;
         theta[6] = 0.01;
         theta[7] = 0.01;
         theta[8] = 0.01;
         theta[9] = 0.01;
         theta[10] = 0.01;
         theta[11] = 0.01;
         theta[12] = 0.01;
         theta[13] = 0.01;
         theta[14] = 0.01;
         theta[15] = 0.01;
         theta[16] = 0.01;
      endsub;

      subroutine dbn_upperbounds(theta[*]);
         outargs theta;
         theta[1] = 0.99;
         theta[2] = 0.99;
         theta[3] = 0.99;
         theta[4] = 0.99;
         theta[5] = 0.99;
         theta[6] = 0.99;
         theta[8] = 0.99;
         theta[9] = 0.99;
         theta[11] = 0.99;
         theta[12] = 0.99;

         theta[7] = 0.99;
         theta[10] = 0.99;
         theta[13] = 0.99;
         theta[14] = 0.99;
         theta[15] = 0.99;
         theta[16] = 0.99;
      endsub;

      subroutine dbn_sampinit(x[*],xlag[*],ylag[*],z[*],zlag[*],theta[*],_time);
         outargs x, xlag;
         a = rand('uniform');
         x[1] = ifn(a<0.5, 1, 0); /*OutStatus*/
         a = rand('uniform');
         x[2] = ifn(a<0.5, 1, 0); /*InnerStatus*/
         x[3] = 57 + 2*rand('NORMAL', 0, 1); /*OutTempRead*/
         x[4] = 66 + 2*rand('NORMAL', 0, 1); /*InTempRead*/
      endsub;

      function dbn_evalmeas(x[*],xlag[*],y[*],ylag[*],z[*],zlag[*],theta[*],_time);
         if x[1] = 0 then lp1 = LOGPDF('NORMAL',y[1],ylag[1],sqrt(theta[15]));
         else lp1 = LOGPDF('NORMAL',y[1],x[3],sqrt(theta[16]));
         if x[2] = 0 then lp2 = LOGPDF('NORMAL',y[2],ylag[2],sqrt(theta[15]));
         else lp2 = LOGPDF('NORMAL',y[2],x[4],sqrt(theta[16]));
         return (lp1+lp2);
      endsub;

      subroutine dbn_samptran(x[*],xlag[*],ylag[*],z[*],zlag[*],theta[*],_time);
         outargs x;
         pi=constant('PI');
         a = rand('uniform');
         if xlag[1] = 0 then x[1] = ifn(a<theta[1],0,1);
         if xlag[1] = 1 then x[1] = ifn(a<theta[2],0,1);
         a = rand('uniform');
         if xlag[2] = 0 then x[2] = ifn(a<theta[3],0,1);
         if xlag[2] = 1 then x[2] = ifn(a<theta[4],0,1);
         x[3] = xlag[3] + 45*pi/287*cos(2*pi/287*(_time-2)-0.5*pi)
            +sqrt(theta[14])*rand('NORMAL', 0, 1);
         if zlag[1] = 1 then x[4] = theta[5]*xlag[4] + theta[6]*xlag[3]
            +sqrt(theta[7])*rand('NORMAL', 0, 1);
         if zlag[1] = 2 then x[4] = theta[8]*xlag[4] + theta[9]*xlag[3]
            + sqrt(theta[10])*rand('NORMAL', 0, 1);
         if zlag[1] = 3 then x[4] = theta[11]*xlag[4] + theta[12]*xlag[3]
            +sqrt(theta[13])*rand('NORMAL', 0, 1);
      endsub;

      function dbn_evalprio(theta[*]);
         return (0);
      endsub;
   endsubmit;
run;

The summary tables of the model information are shown in Figure 35.

The summary tables of the posterior samples are shown in Figure 36.

The correlation and covariance matrices of the parameter samples that the PMCMC method generates are shown in Figure 37.

The following statements load the smoothing estimates for each state variable into the Work library and then sort each data set by the time ID:

data smo_s1; set mylib.sm(where=(variableID=1)); run;
proc sort data=smo_s1; by timeId; run;

data smo_s2; set mylib.sm(where=(variableID=2)); run;
proc sort data=smo_s2; by timeId; run;

data smo_s3; set mylib.sm(where=(variableID=3)); run;
proc sort data=smo_s3; by timeId; run;

data smo_s4; set mylib.sm(where=(variableID=4)); run;
proc sort data=smo_s4; by timeId; run;

The following statements create a relevant data table, combine_s3, in the Work library to record the smoothing estimates and actual values of OutTemp. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Figure 38, where the shadow band indicates the range from the 2.5th to the 97.5th percentiles of the smoothing estimates.

data combine_s3;
   set smo_s3;
   set hvacSystem;
run;

proc sgplot data=combine_s3;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=Time y=OT /legendlabel='Actual Value' lineattrs=(color=red);
   series x=Time y=Q2/ legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
run;

The following statements create a relevant data table, combine_s4, in the Work library to record the smoothing estimates and actual values of InTemp. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot is shown in Figure 39, where the shadow band indicates the range from the 2.5th to the 97.5th percentiles of the smoothing estimates.

data combine_s4;
   set smo_s4;
   set hvacSystem;
run;

proc sgplot data=combine_s4;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=Time y=IT /legendlabel='Actual Value' lineattrs=(color=red);
   series x=Time y=Q2/ legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
run;

Figure 38 and Figure 39 show both the outer temperature and inner temperature are accurately estimated by the model.

For the category variables, OutSensorStatus and InSensorStatus, the smoothing estimates of the frequency of the work status (WORK: 1) are used to compare the actual values. The following statements create a relevant data table, combine_s1, in the Work library to record the smoothing estimates and actual values of the category variable OutSensorStatus. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot, which is shown in Figure 40, indicates that the hidden outer sensor status is successfully decoded by the model.

data combine_s1;
   set smo_s1;
   set hvacSystem;
run;

proc sgplot data=combine_s1;
   series x=Time y=OS;
   series x=Time y=cat_2;
run;

The following statements create a relevant data table, combine_s2, in the Work library to record the smoothing estimates and actual values of the category variable InSensorStatus. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The plot, which is shown in Figure 41, indicates that the hidden inner sensor status is successfully decoded by the model.

data combine_s2;
   set smo_s2;
   set hvacSystem;
run;

proc sgplot data=combine_s2;
   series x=Time y=IS;
   series x=Time y=cat_2;
run;