locked
"Quadrature found zero variance" when dimensions of input data is increased (ex: time or replicates added) RRS feed

  • Question

  • Hello everyone.

    I have written a code to model a network. The data is for network nodes in different times and also each set of data has several replicates.The value of each node is a function of other node values in previous time step. I have currently generated my data (function generateData()) in a way that all node values are a function of the first three node values. My model is a linear regression as you can see and I have added one feature for intercept term.My problem is that it only works for (numReplicate=1 and numtime=2; i.e. for no replicate option and just one time interval). It works perfect with that. But When I increase the number of replicates or times, it won't work anymore. Before implementing, I was thinking having more replicates and time points will help the model inference in thoery. I don't understand why it won't work when I increase them, and why I get zero probabilty after running. I tried to play with my prior on precision parameteres (like: generalizationNoise and slabPrecisionPrior) but it still result in "Quadrature found zero variance" error. Is it something wrong with my implementation? Or I should continue with trying to tune the parameters? Can anyone give me a clue? Any help is appreciated.

    Thank you.

    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using MicrosoftResearch.Infer.Models;
    using MicrosoftResearch.Infer;
    using MicrosoftResearch.Infer.Distributions;
    using MicrosoftResearch.Infer.Maths;
    
    namespace sparseModelTest
    {
        public class ComplexModelOneIntRep
        {
            public ComplexModelOneIntRep()
            {
            }
    
            static int dataSize = 10;
            static int featureSize = 10;
            static int numReplicate = 1;
            static int numtime = 2;
    
            public static void MainTest()
            {            
                runComplexModel();            
            }
    
    		// in network inference, a prior network can be an input to the model
    		// here a fully connected prior network is generated
            static int[][] getFullPriorNetwork(int numNodes)
            {
                int[][] ppiIndices = new int[numNodes][];
                for (int n1 = 0; n1 < numNodes - 1; n1++)
                {
                    ppiIndices[n1] = new int[numNodes];
                    for (int j = 0; j < numNodes; j++)
                    {
                        ppiIndices[n1][j] = j;
                    }
                }
                ppiIndices[numNodes - 1] = new int[0];
    
                return ppiIndices;
            }
    
    		// creating the model
            public static void runComplexModel()
            {
    			// defining ranges
                Range featureRange = new Range(featureSize+1).Named("featureRange");//featureSize+1
                Range dataRange = new Range(dataSize).Named("dataRange");
                Range replicateRange = new Range(numReplicate).Named("replicateRange");
                //////////// 
    			// W shows weight on inferred network edges
    			// ppiIndices will show which w indices should be present in the model to be learned (prior network structure - shows which edges should be inferred)
    			// wIndices for each node, contains the indices of nodes which have an outgoing edge to this node
                int[][] ppiIndices = getFullPriorNetwork(featureRange.SizeAsInt);
                int[] wsizes = new int[ppiIndices.Length];
                for (int i = 0; i < ppiIndices.Length; i++)
                {
                    wsizes[i] = ppiIndices.ElementAt(i).Length;
                }
                VariableArray<int> wsizesVar = Variable.Constant(wsizes, featureRange).Named("wsizes");
                Range wRange = new Range(wsizesVar[featureRange]).Named("wRange");
                var wIndices = Variable.Array(Variable.Array<int>(wRange), featureRange).Named("wIndices");
                wIndices.ObservedValue = ppiIndices;
                ////////////////////       
                // creating seperate observed variables for each dimension's prior on Bernouli distribution
    			// here the prior on the probabilty of the existance of each node can be entered (currentl all eges have the same probabilty of existance)
                double[][] RhoPriorObserved = new double[featureRange.SizeAsInt][];
                for (int i = 0; i < RhoPriorObserved.Length; i++)
                {
                    RhoPriorObserved[i] = new double[ppiIndices[i].Length];
                    for (int j = 0; j < RhoPriorObserved[i].Length; j++)
                    {
                        if (j < 3) RhoPriorObserved[i][j] = 0.3;
                        else RhoPriorObserved[i][j] = 0.3;
                    }
                }
                //VariableArray<double> rhoPriors = Variable.Observed(RhoPriorObserved, featureRange,);
                var rhoPriors = Variable.Array(Variable.Array<double>(wRange), featureRange);
                rhoPriors.ObservedValue = RhoPriorObserved;
                ///////////////////
    
                VariableArray<VariableArray<double>, double[][]> w = Variable.Array(Variable.Array<double>(wRange), featureRange).Named("w");
                VariableArray<double> slabMean = Variable.Array<double>(featureRange).Named("slabMean");
                VariableArray<double> slabPrecision = Variable.Array<double>(featureRange).Named("slabPrecision");
                Gaussian slabMeanPrior = Gaussian.FromMeanAndPrecision(2, 0.1);
                Gamma slabPrecisionPrior = Gamma.FromShapeAndScale(1, 899999999);// Gamma slabNoisePrior = Gamma.FromShapeAndScale(1, 99999);//150000
                using (Variable.ForEach(featureRange))
                {
                    slabMean[featureRange] = Variable.Random(slabMeanPrior);
                    slabPrecision[featureRange] = Variable.Random(slabPrecisionPrior);
                }
                VariableArray<double>[] output = new VariableArray<double>[numtime];//Variable.Array<double>(dataRange).Named("output");
                for (int t = 0; t < numtime; t++)
                {
                    output[t] = Variable.Array<double>(dataRange).Named("output_" + t);
                }
                VariableArray<double> alpha = Variable.Array<double>(dataRange).Named("alpha");
                VariableArray<double> basee = Variable.Array<double>(dataRange).Named("basee");
                VariableArray<VariableArray<VariableArray<double>, double[][]>, double[][][]>[] input = new VariableArray<VariableArray<VariableArray<double>, double[][]>, double[][][]>[numtime];
                for (int t = 0; t < numtime; t++)
                {
                    input[t] = Variable.Array(Variable.Array(Variable.Array<double>(featureRange), replicateRange), dataRange).Named("input_" + t);
                }
                VariableArray<VariableArray<bool>, bool[][]> gamma = Variable.Array(Variable.Array<bool>(wRange), featureRange).Named("gamma");
    
                // define generalization ter (error)
                Variable<double> averageGeneralization = Variable.Observed<double>(0);
                Variable<double> generalizationNoise = Variable.GammaFromShapeAndScale(1, 899999999);
    
                VariableArray<VariableArray<double>, double[][]>[] generalizationTerm = new VariableArray<VariableArray<double>, double[][]>[numtime];//Variable.Array(Variable.Array<double>(featureRange), dataRange).Named("input");
                for (int t = 0; t < numtime; t++)
                {
                    generalizationTerm[t] = Variable.Array(Variable.Array<double>(featureRange), dataRange).Named("generalizationTerm_" + t);
                    using (Variable.ForEach(dataRange))
                    using (Variable.ForEach(featureRange))
                        generalizationTerm[t][dataRange][featureRange] = Variable.GaussianFromMeanAndPrecision(averageGeneralization, generalizationNoise);
                }
    
                gamma[featureRange][wRange] = Variable.Bernoulli(rhoPriors[featureRange][wRange]);
    
                Gaussian alphaPrior = Gaussian.FromMeanAndVariance(0, 0.5);
                using (Variable.ForEach(dataRange)) alpha[dataRange] = Variable<double>.Random(alphaPrior);
                Variable.ConstrainPositive(alpha[dataRange]);
    
                using (Variable.ForEach(featureRange))
                {
                    using (Variable.ForEach(wRange))
                    {
                        using (Variable.If(gamma[featureRange][wRange]))
                        {
                            w[featureRange][wRange] = Variable.GaussianFromMeanAndPrecision(slabMean[featureRange], slabPrecision[featureRange]);//slabPrecision[featureRange]//0.01
                        }
                        using (Variable.IfNot(gamma[featureRange][wRange]))
                        {
                            w[featureRange][wRange] = Variable.GaussianFromMeanAndVariance(0, 0.0000001);//Variable<double>.Random(spike[featureRange])
                        }
                    }
                }
    
                VariableArray<double> nodeSubarray = Variable.Array<double>(wRange).Named("ndoeSubarray");
                for (int t = 1; t < numtime; t++)
                {
                    using (Variable.ForEach(dataRange))
                    {
                        using (Variable.ForEach(replicateRange))
                        {
                            using (Variable.ForEach(featureRange))
                            {
                                VariableArray<double> weightSumArray = Variable.Array<double>(wRange).Named("weightSum_" + t);
                                nodeSubarray = Variable.Subarray(input[t - 1][dataRange][replicateRange], wIndices[featureRange]).Named("nodeSubarray_" + t);
                                weightSumArray[wRange] = (w[featureRange][wRange] * nodeSubarray[wRange]).Named("weightSumArray_" + t);
                                Variable<double> linearEffect = Variable<double>.Sum(weightSumArray).Named("sum_" + t);//old: Variable<double>.Sum(tmp).Named("sum_" + t);
                                input[t][dataRange][replicateRange][featureRange] = linearEffect + generalizationTerm[t][dataRange][featureRange];
                            }
                        }//dataRange
                    }
                }
    
                GenerateData();
                for (int t = 0; t < numtime; t++)
                {
                    input[t].ObservedValue = inputData[t];
                }
                basee.ObservedValue = baseData;
    
    
                InferenceEngine ie = new InferenceEngine();
                //ie.ShowFactorGraph = true;
                ie.NumberOfIterations = 50;
                Bernoulli[][] gammaPosterior = ie.Infer<Bernoulli[][]>(gamma);
                Gaussian[][] wPosterior = ie.Infer<Gaussian[][]>(w);
                Gamma generalizationNoisePosterior = ie.Infer<Gamma>(generalizationNoise);
                Gaussian generalizationMeanPosterior = ie.Infer<Gaussian>(averageGeneralization);
    
                System.IO.StreamWriter file = new System.IO.StreamWriter(@".\sparseModel_0.3.txt");
    
                file.WriteLine("********* gamma posterior ************");
                for (int i = 0; i < gammaPosterior.Length; i++)
                {
                    for (int j = 0; j < gammaPosterior[i].Length; j++)
                    {
                        file.WriteLine("gammaPosterior" + "(" + i + "," + j + ")=" + gammaPosterior[i][j]);
                    }
                }
                file.WriteLine("********* w posterior ************");
                for (int i = 0; i < wPosterior.Length; i++)
                {
                    for (int j = 0; j < wPosterior[i].Length; j++)
                    {
                        file.WriteLine("wPosterior" + "(" + i + "," + j + ")=" + wPosterior[i][j]);
                    }
                }
                file.WriteLine("********* generalization posterior ************");
                file.WriteLine("generalizationMeanPosterior=" + generalizationMeanPosterior);
                file.WriteLine("generalizationNoisePosterior=" + generalizationNoisePosterior);
    
                file.Close();
                testMethod(gammaPosterior);
                Console.ReadLine();
            }
    
            public static void testMethod(Bernoulli[][] gammaPosterior)
            {
                int tp = 0, fp = 0, fn = 0, tn = 0;
                double threshold = 0.5;
                for (int j = 0; j < gammaPosterior.Length - 1; j++)
                {
                    for (int i = 0; i < 3; i++)
                    {
                        if (gammaPosterior[j][i].GetMean() > threshold)
                        {
                            tp++;
                        }
                        else
                        {
                            fn++;
                        }
                    }
                    for (int i = 3; i < gammaPosterior[j].Length; i++)
                    {
                        if (gammaPosterior[j][i].GetMean() > threshold)
                        {
                            fp++;
                            Console.WriteLine(j + "," + i);
                        }
                        else
                        {
                            tn++;
                        }
                    }
                }
                Console.WriteLine("Result: ");
                Console.WriteLine("TP: " + tp + "\nTN: " + tn + "\nFP: " + fp + "\nFN: " + fn);
            }
    
            public static void writeArr(Gaussian[] array, string name/*, StreamWriter resultFile*/)
            {
                // Console.WriteLine();
                for (int i = 0; i < array.Length; i++)
                {
                    Console.WriteLine(name + "(" + i + ")=" + array[i]);
                }
            }
    
            public static void writeArr(Bernoulli[] array, string name/*, StreamWriter resultFile*/)
            {
                // Console.WriteLine();
                for (int i = 0; i < array.Length; i++)
                {
                    Console.WriteLine(name + "(" + i + ")=" + array[i]);
                }
            }
    
            static double[][][][] inputData;
            static double[] baseData;
    
            public static void GenerateData()
            {
    
                Gaussian slabDist = Gaussian.FromMeanAndPrecision(1, 0.01);
                Gaussian slabDist2 = Gaussian.FromMeanAndPrecision(1, 1);
                Gaussian slabDist3 = Gaussian.FromMeanAndPrecision(1, 1);
                Gaussian data = Gaussian.FromMeanAndPrecision(0, 1);
                Gaussian noise = Gaussian.FromMeanAndVariance(0, 0.05);
    
                inputData = new double[numtime][][][];
                //outputData = new double[numtime][];
                baseData = new double[numtime];
    
                Rand.Restart(12347);
                inputData[0] = new double[dataSize][][];
                for (int d = 0; d < dataSize; d++)
                {
                    inputData[0][d] = new double[numReplicate][];
                    for (int r = 0; r < numReplicate; r++)
                    {
                        //Console.WriteLine("t, d, f = " + 0 + " , " + nData + " , ");
                        inputData[0][d][r] = new double[featureSize + 1];//featureSize+1
                    }
                    // for replicate 0, initialize the data
                    for (int f = 0; f < featureSize; f++)
                    {
                        inputData[0][d][0][f] = data.Sample();
                    }
                    inputData[0][d][0][featureSize] = 1;
                    //inputData[0][d][featureSize] = 1;
                }
    
                for (int t = 1; t < numtime; t++)
                {
                    inputData[t] = new double[dataSize][][];
                    for (int d = 0; d < dataSize; d++)
                    {
                        inputData[t][d] = new double[numReplicate][];
                        for (int r = 0; r < numReplicate; r++)
                        {
                            inputData[t][d][r] = new double[featureSize + 1];//featureSize+1
                        }
                        for (int j = 0; j < featureSize; j++)
                        {
                            inputData[t][d][0][j] = 2 * inputData[t - 1][d][0][0] + 3 * inputData[t - 1][d][0][1] + 2 * inputData[t - 1][d][0][2]; ///*+ baseData[i] * 0.2 /*this is for simulation of self-effect/                        
                        }
                        inputData[t][d][0][featureSize] = 1;
    
                    }
                }
    
                // generate replicate data for each experiment by adding noise to replicate (r=0)
                for (int r = 1; r < numReplicate; r++)
                {
                    for (int t = 0; t < numtime; t++)
                    {
                        for (int d = 0; d < dataSize; d++)
                        {
                            for (int f = 0; f < featureSize; f++)
                            {
                                inputData[t][d][r][f] = inputData[t][d][0][f] + noise.Sample();
                            }
                            inputData[t][d][r][featureSize] = 1;
                        }
                    }
                }
    
                // add noise to replicate r0 itself
                for (int t = 0; t < numtime; t++)
                {
                    for (int d = 0; d < dataSize; d++)
                    {
                        for (int f = 0; f < featureSize; f++)
                        {
                            inputData[t][d][0][f] = inputData[t][d][0][f] + noise.Sample();
                        }
                    }
                }
                Console.WriteLine("Hi :-)");
                //baseData[0] = outputData[dataSize - 1];
            }
    
        }
    }

    Monday, December 8, 2014 2:37 PM

Answers

  • Looks like a convergence issue.  When I run it, w doesn't converge.
    • Marked as answer by Capli19 Tuesday, December 16, 2014 6:33 PM
    Friday, December 12, 2014 5:14 PM
    Owner

All replies

  • I'm not seeing any problem.  Are you using version 2.6?
    Monday, December 8, 2014 5:05 PM
    Owner
  • No I'm using the previous version. I will check it with 2.6 now. Thank you.
    Monday, December 8, 2014 6:03 PM
  • Thanks it's working now. Still another question: Is it the right way to add replicate and time dimensions to my data?

    When I increase replicates and time points, the results are worst. That's while it should be better in theory.

    For example when I compare this case with three replicates:

            static int dataSize = 20;
            static int featureSize = 10;
            static int numReplicate = 3;
            static int numtime = 2;

    I infer a network of zero edges (all probabilities for gamma are very close to zero).

    But if I use the following setting with 20 data rows and only one replicate I get the correct network:

            static int dataSize = 20;
            static int featureSize = 10;
            static int numReplicate = 1;
            static int numtime = 2;

    Is there technically something wrong in my implementation? Is there something wrong with adding to the dimensions of my descriptor variables?
    when I have more than one replicates, it doesn't matter how much I increase my dataSize, the result won't improve and it's always zero in all elements of gamma.

    Thanks in advance for any help.



    • Edited by Capli19 Tuesday, December 9, 2014 4:25 PM
    Tuesday, December 9, 2014 4:09 PM
  • One more thing:

    When I have three dimensions in my input data (data x replicate x time), I get different result than when I just paste all of replicate data along with data dimension alltogether.

    I mean when I have two dimensions (dataAndReplicatesTogether x time).

    Is it possible to get different result when I have different implementations (like a 2d or 3d input data)?

    Wednesday, December 10, 2014 2:21 PM
  • Looks like a convergence issue.  When I run it, w doesn't converge.
    • Marked as answer by Capli19 Tuesday, December 16, 2014 6:33 PM
    Friday, December 12, 2014 5:14 PM
    Owner