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];
}
}
}
