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

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 PMOwner
All replies
-
I'm not seeing any problem. Are you using version 2.6?Monday, December 8, 2014 5:05 PMOwner
-
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 PMOwner