# Differences Between Inference Engines in Regression Model

• ### Question

• Hello Infer.net Team,

The code below is designed to perform a regression analysis to infer the statistics of hidden variables based on a large number of measurements over a set of selected observable values. Some questions:

1. The model appears to run well under Gibbs Sampling and Variational Message passing, but blows up under EP. Is there any way to address this?

2. Although the mean values for the results between GS and VMP appear to be relatively consistent, the error terms are different. Gibbs appears to yield the higher precision on average, although this is not uniform. Which is to be believed?

3. I am thinking of adding some additional models, including some based on Bernoulli distributions and possibly probit models. Is there a recommended way to proceed to ensure that all models will come together in stable form?

Thanks for any help.

Philip

namespace test
{
class Program
{
static void Main(string[] args)
{
// Define observed data

double[] PS01D_03S_PR_obs = { -10.0, -50.0, 0.0, 0.0, 30.0, 22.0, -5.0, 10.0, 0.0, 30.0 };
double[] PS03D_04S_PR_obs = { 510.0, 250.0, 3.0, 0.0, 0.0, -92.0, -50.0, 0.0, 0.0, 201.0 };
double[] PS04D_05S_PR_obs = { -3010.0, -7050.0, -1109.0, -10.0, 0.0, 207.0, 2850.0, 3900.0, 1690.0, 372.0 };
double[] PS01D_03S_VB_obs = { -3.0, -5.0, 2.0, -8.0, 3.0, 8.0, -1.0, 4.0, 0.0, 8.0 };
double[] PS03D_04S_VB_obs = { -60.0, -10.0, 10.0, 12.0, 10.0, 8.0, -5.0, 12.0, 0.0, 62.0 };
double[] PS04D_05S_VB_obs = { 300.0, 559.0, 360.0, 20.0, -40.0, 12.0, -105.0, -310.0, -462.0, -330.0 };
double[] PS03_VB_obs = { 220.0, -349.0, -90.0, 160.0, 30.0, -112.0, 75.0, -94.0, 62.0, 130.0 };
double[] PS04_VB_obs = { 30.0, 25.0, -5.0, -72.0, -12.0, 42.0, 19.0, 44.0, -37.0, 12.0 };

int n = PS01D_03S_PR_obs.Length;    // snapshot array length

// Priors on statistical distributions

double sigsqBeta = 1.0; // prior on variance for means
double A = 1.1;     // Prior on shape factor for Gamma used in calculation of variances
double B = .9;      // Prior on scale factor for Gamma used in calculation of variances

Range index = new Range(n).Named("index");      // range based on snapshot length

//Variable<bool> model1 = Variable.Bernoulli(0.5).Named("model1");

// Meter priors

Variable<double> PS01D_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS01D_mean");
Variable<double> PS01D_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS01D_precision");
Variable<double> PS03S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS03S_mean");
Variable<double> PS03S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS03S_precision");
Variable<double> PS03T_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS03T_mean");
Variable<double> PS03T_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS03T_precision");
Variable<double> PS03D_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS03D_mean");
Variable<double> PS03D_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS03D_precision");
Variable<double> PS04S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS04S_mean");
Variable<double> PS04S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS04S_precision");
Variable<double> PS04T_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS04T_mean");
Variable<double> PS04T_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS04T_precision");
Variable<double> PS04D_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS04D_mean");
Variable<double> PS04D_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS04D_precision");
Variable<double> PS05S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS05S_mean");
Variable<double> PS05S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS05S_precision");

// Rate priors:

Variable<double> PS01D_PS03S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS01D_PS03S_mean");
Variable<double> PS01D_PS03S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS01D_PS03S_precision");
Variable<double> PS03D_PS04S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS03D_PS04S_mean");
Variable<double> PS03D_PS04S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS03D_PS04S_precision");
Variable<double> PS04D_PS05S_mean = Variable.GaussianFromMeanAndVariance(0.0, sigsqBeta).Named("PS04D_PS05S_mean");
Variable<double> PS04D_PS05S_precision = Variable.GammaFromShapeAndScale(A, B).Named("PS04D_PS05S_precision");

// Meter errors

VariableArray<double> PS01D_Flo_err = Variable.Array<double>(index).Named("PS01D_Flo_err");
VariableArray<double> PS03S_Flo_err = Variable.Array<double>(index).Named("PS03S_Flo_err");
VariableArray<double> PS03T_Flo_err = Variable.Array<double>(index).Named("PS03T_Flo_err");
VariableArray<double> PS03D_Flo_err = Variable.Array<double>(index).Named("PS03D_Flo_err");
VariableArray<double> PS04S_Flo_err = Variable.Array<double>(index).Named("PS04S_Flo_err");
VariableArray<double> PS04T_Flo_err = Variable.Array<double>(index).Named("PS04T_Flo_err");
VariableArray<double> PS04D_Flo_err = Variable.Array<double>(index).Named("PS04D_Flo_err");
VariableArray<double> PS05S_Flo_err = Variable.Array<double>(index).Named("PS05S_Flo_err");

// Rate error parameters

VariableArray<double> PS01D_PS03S_PR_err = Variable.Array<double>(index).Named("PS01D_PS03S_PR_err");
VariableArray<double> PS03D_PS04S_PR_err = Variable.Array<double>(index).Named("PS03D_PS04S_PR_err");
VariableArray<double> PS04D_PS05S_PR_err = Variable.Array<double>(index).Named("PS04D_PS05S_PR_err");

// VB errors

VariableArray<double> PS01D_PS03S_VB = Variable.Array<double>(index).Named("PS01D_PS03S_VB");
VariableArray<double> PS03D_PS04S_VB = Variable.Array<double>(index).Named("PS03D_PS04S_VB");
VariableArray<double> PS04D_PS05S_VB = Variable.Array<double>(index).Named("PS04D_PS05S_VB");

VariableArray<double> PS03_VB = Variable.Array<double>(index).Named("PS03_VB");
VariableArray<double> PS04_VB = Variable.Array<double>(index).Named("PS04_VB");

// Rates

VariableArray<double> PS01D_PS03S_PR = Variable.Array<double>(index).Named("PS01D_PS03S_PR");
VariableArray<double> PS03D_PS04S_PR = Variable.Array<double>(index).Named("PS03D_PS04S_PR");
VariableArray<double> PS04D_PS05S_PR = Variable.Array<double>(index).Named("PS04D_PS05S_PR");

// Error sum parameters

VariableArray<double> PS01D_PS03S_mu = Variable.Array<double>(index).Named("PS01D_PS03S_mu");
VariableArray<double> PS03D_PS04S_mu = Variable.Array<double>(index).Named("PS03D_PS04S_mu");
VariableArray<double> PS04D_PS05S_mu = Variable.Array<double>(index).Named("PS04D_PS05S_mu");

VariableArray<double> PS03_mu = Variable.Array<double>(index).Named("PS03_mu");
VariableArray<double> PS04_mu = Variable.Array<double>(index).Named("PS04_mu");

// Flow errors drawn from Gaussian based on meter priors

PS01D_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS01D_mean, PS01D_precision).ForEach(index);
PS03S_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS03S_mean, PS03S_precision).ForEach(index);
PS03T_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS03T_mean, PS03T_precision).ForEach(index);
PS03D_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS03D_mean, PS03D_precision).ForEach(index);
PS04S_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS04S_mean, PS04S_precision).ForEach(index);
PS04T_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS04T_mean, PS04T_precision).ForEach(index);
PS04D_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS04D_mean, PS04D_precision).ForEach(index);
PS05S_Flo_err[index] = Variable.GaussianFromMeanAndPrecision(PS05S_mean, PS05S_precision).ForEach(index);

// Rate parameter index errors based on priors:

PS01D_PS03S_PR_err[index] = Variable.GaussianFromMeanAndPrecision(PS01D_PS03S_mean, PS01D_PS03S_precision).ForEach(index);
PS03D_PS04S_PR_err[index] = Variable.GaussianFromMeanAndPrecision(PS03D_PS04S_mean, PS03D_PS04S_precision).ForEach(index);
PS04D_PS05S_PR_err[index] = Variable.GaussianFromMeanAndPrecision(PS04D_PS05S_mean, PS04D_PS05S_precision).ForEach(index);

// Residual error terms for VBs

var PS01D_PS03S_tau = Variable.GammaFromShapeAndScale(A, 1 / B).Named("PS01D_PS03S_tau");
var PS03D_PS04S_tau = Variable.GammaFromShapeAndScale(A, 1 / B).Named("PS03D_PS04S_tau");
var PS04D_PS05S_tau = Variable.GammaFromShapeAndScale(A, 1 / B).Named("PS04D_PS05S_tau");

var PS03_tau = Variable.GammaFromShapeAndScale(A, 1 / B).Named("PS03_tau");
var PS04_tau = Variable.GammaFromShapeAndScale(A, 1 / B).Named("PS04_tau");

// Model

PS01D_PS03S_mu[index] = PS01D_Flo_err[index] + PS03S_Flo_err[index] + PS01D_PS03S_PR_err[index] * PS01D_PS03S_PR[index];
PS03D_PS04S_mu[index] = PS03D_Flo_err[index] + PS04S_Flo_err[index] + PS03D_PS04S_PR_err[index] * PS03D_PS04S_PR[index];
PS04D_PS05S_mu[index] = PS04D_Flo_err[index] + PS05S_Flo_err[index] + PS04D_PS05S_PR_err[index] * PS04D_PS05S_PR[index];
PS03_mu[index] = PS03S_Flo_err[index] + PS03D_Flo_err[index] + PS03T_Flo_err[index];
PS04_mu[index] = PS04S_Flo_err[index] + PS04D_Flo_err[index] + PS04T_Flo_err[index];

PS01D_PS03S_VB[index] = Variable.GaussianFromMeanAndPrecision(PS01D_PS03S_mu[index], PS01D_PS03S_tau);
PS03D_PS04S_VB[index] = Variable.GaussianFromMeanAndPrecision(PS03D_PS04S_mu[index], PS03D_PS04S_tau);
PS04D_PS05S_VB[index] = Variable.GaussianFromMeanAndPrecision(PS04D_PS05S_mu[index], PS04D_PS05S_tau);
PS03_VB[index] = Variable.GaussianFromMeanAndPrecision(PS03_mu[index], PS03_tau);
PS04_VB[index] = Variable.GaussianFromMeanAndPrecision(PS04_mu[index], PS04_tau);

// Observations:

PS01D_PS03S_PR.ObservedValue = PS01D_03S_PR_obs;
PS03D_PS04S_PR.ObservedValue = PS03D_04S_PR_obs;
PS04D_PS05S_PR.ObservedValue = PS04D_05S_PR_obs;

PS01D_PS03S_VB.ObservedValue = PS01D_03S_VB_obs;
PS03D_PS04S_VB.ObservedValue = PS03D_04S_VB_obs;
PS04D_PS05S_VB.ObservedValue = PS04D_05S_VB_obs;
PS03_VB.ObservedValue = PS03_VB_obs;
PS04_VB.ObservedValue = PS04_VB_obs;

var engine = new InferenceEngine(new GibbsSampling());
//var engine = new InferenceEngine(new VariationalMessagePassing());
//var engine = new InferenceEngine();
//Console.WriteLine(engine.Infer(model1));

//Console.WriteLine("beta0=" + engine.Infer(beta0));
Console.WriteLine("PS01D_mean =" + engine.Infer(PS01D_mean));
Console.WriteLine("PS01D_precision =" + engine.Infer(PS01D_precision));
//Console.WriteLine("beta1=" + engine.Infer(beta1));
Console.WriteLine("PS03S_mean =" + engine.Infer(PS03S_mean));
Console.WriteLine("PS03S_precision =" + engine.Infer(PS03S_precision));
//Console.WriteLine("beta2=" + engine.Infer(beta2));
Console.WriteLine("PS03T_mean =" + engine.Infer(PS03T_mean));
Console.WriteLine("PS03T_precision =" + engine.Infer(PS03T_precision));

Console.WriteLine("PS03D_mean =" + engine.Infer(PS03D_mean));
Console.WriteLine("PS03D_precision =" + engine.Infer(PS03D_precision));

Console.WriteLine("PS04S_mean =" + engine.Infer(PS04S_mean));
Console.WriteLine("PS04S_precision =" + engine.Infer(PS04S_precision));

Console.WriteLine("PS04T_mean =" + engine.Infer(PS04T_mean));
Console.WriteLine("PS04T_precision =" + engine.Infer(PS04T_precision));

Console.WriteLine("PS04D_mean =" + engine.Infer(PS04D_mean));
Console.WriteLine("PS04D_precision =" + engine.Infer(PS04D_precision));

Console.WriteLine("PS05S_mean =" + engine.Infer(PS05S_mean));
Console.WriteLine("PS05S_precision =" + engine.Infer(PS05S_precision));

Console.WriteLine("PS01D_PS03S_mean =" + engine.Infer(PS01D_PS03S_mean));
Console.WriteLine("PS01D_PS03S_precision =" + engine.Infer(PS01D_PS03S_precision));

Console.WriteLine("PS03D_PS04S_mean =" + engine.Infer(PS03D_PS04S_mean));
Console.WriteLine("PS03D_PS04S_precision =" + engine.Infer(PS03D_PS04S_precision));

Console.WriteLine("PS04D_PS05S_mean =" + engine.Infer(PS04D_PS05S_mean));
Console.WriteLine("PS04D_PS05S_precision =" + engine.Infer(PS04D_PS05S_precision));

Console.WriteLine("PS01D_PS03S_tau =" + engine.Infer(PS01D_PS03S_tau));
Console.WriteLine("PS03D_PS04S_tau =" + engine.Infer(PS03D_PS04S_tau));
Console.WriteLine("PS04D_PS05S_tau =" + engine.Infer(PS04D_PS05S_tau));
Console.WriteLine("PS03_tau =" + engine.Infer(PS03_tau));
Console.WriteLine("PS04_tau =" + engine.Infer(PS04_tau));

}
}
}

Tuesday, June 12, 2012 10:42 PM

### All replies

• 1. This is probably due to your priors, which have been set very far from your data.  Based on your data, you'd get better convergence by setting sigsqBeta much higher (say 100*100) and precisions to have a mean of 1/100 with shape parameter A around 10.  I notice that in some places you have GammaFromShapeAndScale(A, B) versus GammaFromShapeAndScale(A, 1 / B).   You may find it more intuitive to use GammaFromMeanAndVariance.

2. Both algorithms are approximate so the only way to be sure is to test the estimates on a holdout set.

3. Don't understand the question.

Thursday, June 14, 2012 7:26 PM