locked
Differences Between Inference Engines in Regression Model RRS feed

  • 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));

                Console.Read();
            }
        }
    }

    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
    Owner