C# Multivariate.MultivariateNormalDistribution类代码示例

        public void MultivariateNormalGenerateTest()
        {
            // mean vector
            double[] mu = { 2.0, 6.0 };

            // covariance
            double[,] cov = 
            {
                { 2, 1 },
                { 1, 5 } 
            };

            // Create a multivariate Normal distribution
            var normal = new MultivariateNormalDistribution(mu, cov);

            // Generate 1000000 samples from it
            double[][] samples = normal.Generate(1000000);

            // Try to estimate a new Normal distribution from
            // generated samples to check if they indeed match
            var actual = MultivariateNormalDistribution.Estimate(samples);

            Assert.IsTrue(mu.IsEqual(actual.Mean, 0.1));
            Assert.IsTrue(cov.IsEqual(actual.Covariance, 0.1));
        }

        public static HiddenMarkovModel<MultivariateNormalDistribution> deserialize(string filename)
        {
            using (StreamReader r = new StreamReader(filename))
            {
                String data = r.ReadToEnd();
                String[] dataArr = data.Split('#');

                // probDelimitedStr dataArr[0]
                double[] prob = createSingleDimDoubleArray(dataArr[0], '|');

                // transDelimitedStr dataArr[1]
                double[,] trans = createDoubleDimDoubleArray(dataArr[1], '|', '*');

                // emissionsDelimitedStr dataArr[2]
                String[] emissions = dataArr[2].Split('$');

                MultivariateNormalDistribution[] e2 = new MultivariateNormalDistribution[emissions.Length];
                for (int i = 0; i < emissions.Length; i++)
                {
                    String[] meansNCovariance = emissions[i].Split('&');
                    String meansStr = meansNCovariance[0];
                    String covarianceStr = meansNCovariance[1];

                    double[] means = createSingleDimDoubleArray(meansStr, '|');
                    double[,] covariance = createDoubleDimDoubleArray(covarianceStr, '|', '*');
                    MultivariateNormalDistribution dist = new MultivariateNormalDistribution(means, covariance);
                    e2[i] = dist;
                }
                HiddenMarkovModel<MultivariateNormalDistribution> hmm = new HiddenMarkovModel<MultivariateNormalDistribution>(trans, e2, prob);
                return hmm;
            }
        }

        /// <summary>
        ///   Creates the initial scatter plot graphs containing some random
        ///   data. This data is generated by sampling Gaussian distributions.
        /// </summary>
        /// 
        private void btnGenerateRandom_Click(object sender, EventArgs e)
        {
            k = (int)numClusters.Value;

            // Generate data with n Gaussian distributions
            double[][][] data = new double[k][][];

            for (int i = 0; i < k; i++)
            {
                // Create random centroid to place the Gaussian distribution
                double[] mean = Matrix.Random(2, -6.0, +6.0);

                // Create random covariance matrix for the distribution
                double[,] covariance = Accord.Statistics.Tools.RandomCovariance(2, -5, 5);

                // Create the Gaussian distribution
                var gaussian = new MultivariateNormalDistribution(mean, covariance);

                int samples = Accord.Math.Tools.Random.Next(150, 250);
                data[i] = gaussian.Generate(samples);
            }

            // Join the generated data
            mixture = Matrix.Stack(data);

            // Update the scatter plot
            CreateScatterplot(graph, mixture, k);

            // Forget previous initialization
            kmeans = null;
        }

        public static HiddenMarkovClassifier<MultivariateNormalDistribution> CreateModel1()
        {
            // Create a Continuous density Hidden Markov Model Sequence Classifier
            // to detect a multivariate sequence and the same sequence backwards.
            double[][][] sequences = new double[][][]
            {
                new double[][] 
                { 
                    // This is the first  sequence with label = 0
                    new double[] { 0 },
                    new double[] { 1 },
                    new double[] { 2 },
                    new double[] { 3 },
                    new double[] { 4 },
                }, 

                new double[][]
                {
                     // This is the second sequence with label = 1
                    new double[] { 4 },
                    new double[] { 3 },
                    new double[] { 2 },
                    new double[] { 1 },
                    new double[] { 0 },
                }
            };

            // Labels for the sequences
            int[] labels = { 0, 1 };

            // Creates a sequence classifier containing 2 hidden Markov Models
            //  with 2 states and an underlying Normal distribution as density.
            MultivariateNormalDistribution density = new MultivariateNormalDistribution(1);
            var classifier = new HiddenMarkovClassifier<MultivariateNormalDistribution>(2, new Ergodic(2), density);

            // Configure the learning algorithms to train the sequence classifier
            var teacher = new HiddenMarkovClassifierLearning<MultivariateNormalDistribution>(classifier,

                // Train each model until the log-likelihood changes less than 0.001
                modelIndex => new BaumWelchLearning<MultivariateNormalDistribution>(classifier.Models[modelIndex])
                {
                    Tolerance = 0.0001,
                    Iterations = 0
                }
            );

            // Train the sequence classifier using the algorithm
            double logLikelihood = teacher.Run(sequences, labels);


            return classifier;
        }

        public void ConstructorTest1()
        {
            MultivariateNormalDistribution[] components = new MultivariateNormalDistribution[2];
            components[0] = new MultivariateNormalDistribution(new double[] { 2 }, new double[,] { { 1 } });
            components[1] = new MultivariateNormalDistribution(new double[] { 5 }, new double[,] { { 1 } });

            var mixture = new MultivariateMixture<MultivariateNormalDistribution>(components);

            double[] expected = { 0.5, 0.5 };

            Assert.IsTrue(expected.IsEqual(mixture.Coefficients));
            Assert.AreEqual(components, mixture.Components);
        }

        public void ConstructorTest4()
        {
            // Create a multivariate Gaussian distribution 
            var dist = new MultivariateNormalDistribution(

                // mean vector mu
                mean: new double[] 
                {
                    4, 2 
                },

                // covariance matrix sigma
                covariance: new double[,] 
                {
                    { 0.3, 0.1 },
                    { 0.1, 0.7 }
                }
            );

            // Common measures
            double[] mean = dist.Mean;     // { 4, 2 }
            double[] median = dist.Median; // { 4, 2 }
            double[] var = dist.Variance;  // { 0.3, 0.7 } (diagonal from cov)
            double[,] cov = dist.Covariance; // { { 0.3, 0.1 }, { 0.1, 0.7 } }

            // Probability mass functions
            double pdf1 = dist.ProbabilityDensityFunction(new double[] { 2, 5 }); // 0.000000018917884164743237
            double pdf2 = dist.ProbabilityDensityFunction(new double[] { 4, 2 }); // 0.35588127170858852
            double pdf3 = dist.ProbabilityDensityFunction(new double[] { 3, 7 }); // 0.000000000036520107734505265
            double lpdf = dist.LogProbabilityDensityFunction(new double[] { 3, 7 }); // -24.033158110192296

            // Cumulative distribution function (for up to two dimensions)
            double cdf = dist.DistributionFunction(new double[] { 3, 5 }); // 0.033944035782101548


            Assert.AreEqual(4, mean[0]);
            Assert.AreEqual(2, mean[1]);
            Assert.AreEqual(4, median[0]);
            Assert.AreEqual(2, median[1]);
            Assert.AreEqual(0.3, var[0]);
            Assert.AreEqual(0.7, var[1]);
            Assert.AreEqual(0.3, cov[0, 0]);
            Assert.AreEqual(0.1, cov[0, 1]);
            Assert.AreEqual(0.1, cov[1, 0]);
            Assert.AreEqual(0.7, cov[1, 1]);
            Assert.AreEqual(0.000000018917884164743237, pdf1);
            Assert.AreEqual(0.35588127170858852, pdf2);
            Assert.AreEqual(0.000000000036520107734505265, pdf3);
            Assert.AreEqual(-24.033158110192296, lpdf);
            Assert.AreEqual(0.033944035782101548, cdf);
        }

        public void ProbabilityDensityFunctionTest()
        {
            MultivariateNormalDistribution[] components = new MultivariateNormalDistribution[2];
            components[0] = new MultivariateNormalDistribution(new double[] { 2 }, new double[,] { { 1 } });
            components[1] = new MultivariateNormalDistribution(new double[] { 5 }, new double[,] { { 1 } });

            double[] coefficients = { 0.3, 0.7 };
            var mixture = new MultivariateMixture<MultivariateNormalDistribution>(coefficients, components);

            double[] x = { 1.2 };

            double expected =
                0.3 * components[0].ProbabilityDensityFunction(x) +
                0.7 * components[1].ProbabilityDensityFunction(x);

            double actual = mixture.ProbabilityDensityFunction(x);

            Assert.AreEqual(expected, actual);
        }

        /*public static void SaveSequenceList(SequenceList seqList, string path)
        {
            Stream writeStream = new FileStream(path, FileMode.Create, FileAccess.Write, FileShare.None);
            seqList.Save(writeStream);
            writeStream.Close();
        }

        public static SequenceList LoadSequenceList(string path)
        {
            Stream readStream = new FileStream(path, FileMode.Open, FileAccess.Read, FileShare.Read);
            SequenceList seqList = SequenceList.Load(readStream);
            readStream.Close();
            return seqList;
        }*/
        public static HiddenMarkovModel<MultivariateNormalDistribution> CreateModelFromFrames(List<List<Frame>> frames)
        {
            SequenceList sequences = Utils.FramesToSequenceList(frames);

            HiddenMarkovModel<MultivariateNormalDistribution> hmm;
            MultivariateNormalDistribution mnd = new MultivariateNormalDistribution(sequences.GetDimensions());
            hmm = new HiddenMarkovModel<MultivariateNormalDistribution>(new Forward(5), mnd);

            var teacher = new BaumWelchLearning<MultivariateNormalDistribution>(hmm);
            teacher.Tolerance = 0.0001;
            teacher.Iterations = 0;
            teacher.FittingOptions = new NormalOptions()
            {
                Diagonal = true,      // only diagonal covariance matrices
                Regularization = 1e-5 // avoid non-positive definite errors
            };

            teacher.Run(sequences.GetArray());

            return hmm;
        }

        public void ConstructorTest1()
        {

            NormalDistribution normal = new NormalDistribution(4.2, 1.2);
            MultivariateNormalDistribution target = new MultivariateNormalDistribution(new[] { 4.2 }, new[,] { { 1.2 * 1.2 } });

            double[] mean = target.Mean;
            double[] median = target.Median;
            double[] var = target.Variance;
            double[,] cov = target.Covariance;

            double apdf1 = target.ProbabilityDensityFunction(new double[] { 2 });
            double apdf2 = target.ProbabilityDensityFunction(new double[] { 4 });
            double apdf3 = target.ProbabilityDensityFunction(new double[] { 3 });
            double alpdf = target.LogProbabilityDensityFunction(new double[] { 3 });
            double acdf = target.DistributionFunction(new double[] { 3 });
            double accdf = target.ComplementaryDistributionFunction(new double[] { 3 });

            double epdf1 = target.ProbabilityDensityFunction(new double[] { 2 });
            double epdf2 = target.ProbabilityDensityFunction(new double[] { 4 });
            double epdf3 = target.ProbabilityDensityFunction(new double[] { 3 });
            double elpdf = target.LogProbabilityDensityFunction(new double[] { 3 });
            double ecdf = target.DistributionFunction(new double[] { 3 });
            double eccdf = target.ComplementaryDistributionFunction(new double[] { 3 });


            Assert.AreEqual(normal.Mean, target.Mean[0]);
            Assert.AreEqual(normal.Median, target.Median[0]);
            Assert.AreEqual(normal.Variance, target.Variance[0]);
            Assert.AreEqual(normal.Variance, target.Covariance[0, 0]);

            Assert.AreEqual(epdf1, apdf1);
            Assert.AreEqual(epdf2, apdf2);
            Assert.AreEqual(epdf3, apdf3);
            Assert.AreEqual(elpdf, alpdf);
            Assert.AreEqual(ecdf, acdf);
            Assert.AreEqual(eccdf, accdf);
            Assert.AreEqual(1.0 - ecdf, eccdf);
        }

        public void MixtureWeightsFitTest()
        {
            // Randomly initialize some mixture components
            MultivariateNormalDistribution[] components = new MultivariateNormalDistribution[2];
            components[0] = new MultivariateNormalDistribution(new double[] { 2 }, new double[,] { { 1 } });
            components[1] = new MultivariateNormalDistribution(new double[] { 5 }, new double[,] { { 1 } });

            // Create an initial mixture
            var mixture = new MultivariateMixture<MultivariateNormalDistribution>(components);

            // Now, suppose we have a weighted data
            // set. Those will be the input points:

            double[][] points = new double[] { 0, 3, 1, 7, 3, 5, 1, 2, -1, 2, 7, 6, 8, 6 } // (14 points)
                .ToArray();

            // And those are their respective unormalized weights:
            double[] weights = { 1, 1, 1, 2, 2, 1, 1, 1, 2, 1, 2, 3, 1, 1 }; // (14 weights)

            // Let's normalize the weights so they sum up to one:
            weights = weights.Divide(weights.Sum());

            // Now we can fit our model to the data:
            mixture.Fit(points, weights);   // done!

            // Our model will be:
            double mean1 = mixture.Components[0].Mean[0]; // 1.41126
            double mean2 = mixture.Components[1].Mean[0]; // 6.53301

            // If we need the GaussianMixtureModel functionality, we can
            // use the estimated mixture to initialize a new model:
            GaussianMixtureModel gmm = new GaussianMixtureModel(mixture);

            Assert.AreEqual(mean1, gmm.Gaussians[0].Mean[0]);
            Assert.AreEqual(mean2, gmm.Gaussians[1].Mean[0]);

            Assert.AreEqual(1.4112610766836404, mean1, 1e-10);
            Assert.AreEqual(6.5330177004151082, mean2, 1e-10);

            Assert.AreEqual(mixture.Coefficients[0], gmm.Gaussians[0].Proportion);
            Assert.AreEqual(mixture.Coefficients[1], gmm.Gaussians[1].Proportion);
        }

        public void MixtureWeightsFitTest2()
        {
            MemoryStream stream = new MemoryStream(Resources.CircleWithWeights);

            ExcelReader reader = new ExcelReader(stream, xlsx: false);

            DataTable table = reader.GetWorksheet("Sheet1");

            double[,] matrix = table.ToMatrix();

            double[][] points = matrix.Submatrix(null, 0, 1).ToArray();
            double[] weights = matrix.GetColumn(2);

            // Randomly initialize some mixture components
            MultivariateNormalDistribution[] components = new MultivariateNormalDistribution[2];
            components[0] = new MultivariateNormalDistribution(new double[] { 0, 1 }, Matrix.Identity(2));
            components[1] = new MultivariateNormalDistribution(new double[] { 1, 0 }, Matrix.Identity(2));

            // Create an initial mixture
            var mixture = new MultivariateMixture<MultivariateNormalDistribution>(components);

            mixture.Fit(points, weights);

            // Our model will be:
            double mean00 = mixture.Components[0].Mean[0];
            double mean01 = mixture.Components[0].Mean[1];
            double mean10 = mixture.Components[1].Mean[0];
            double mean11 = mixture.Components[1].Mean[1];

            Assert.AreEqual(-0.11704994950834195, mean00, 1e-10);
            Assert.AreEqual(0.11603470123007256, mean01, 1e-10);
            Assert.AreEqual(0.11814483652855159, mean10, 1e-10);
            Assert.AreEqual(-0.12029275652994373, mean11, 1e-10);
        }

        public void ProbabilityDensityFunctionTest()
        {
            double[] mean = { 1, -1 };
            double[,] covariance = 
            {
                { 0.9, 0.4 },
                { 0.4, 0.3 },
            };

            var target = new MultivariateNormalDistribution(mean, covariance);

            double[] x = { 1.2, -0.8 };
            double expected = 0.446209421363460;
            double actual = target.ProbabilityDensityFunction(x);

            Assert.AreEqual(expected, actual, 0.00000001);
        }

        public void FitTest2()
        {
            double[] coefficients = { 0.50, 0.50 };

            MultivariateNormalDistribution[] components = new MultivariateNormalDistribution[2];
            components[0] = new MultivariateNormalDistribution(new double[] { 2 }, new double[,] { { 1 } });
            components[1] = new MultivariateNormalDistribution(new double[] { 5 }, new double[,] { { 1 } });

            var target = new MultivariateMixture<MultivariateNormalDistribution>(coefficients, components);

            double[][] values = { new double[] { 2512512312 },
                                  new double[] { 1 }, 
                                  new double[] { 1 },
                                  new double[] { 0 },
                                  new double[] { 1 },
                                  new double[] { 6 },
                                  new double[] { 6 },
                                  new double[] { 5 },
                                  new double[] { 7 },
                                  new double[] { 5 } };

            double[] weights = { 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
            weights = weights.Divide(weights.Sum());

            double[][] part1 = values.Submatrix(1, 4);
            double[][] part2 = values.Submatrix(5, 9);


            target.Fit(values, weights);

            var mean1 = Accord.Statistics.Tools.Mean(part1);
            var var1 = Accord.Statistics.Tools.Variance(part1);
            Assert.AreEqual(mean1[0], target.Components[0].Mean[0], 1e-5);
            Assert.AreEqual(var1[0], target.Components[0].Variance[0], 1e-5);

            var mean2 = Accord.Statistics.Tools.Mean(part2);
            var var2 = Accord.Statistics.Tools.Variance(part2);
            Assert.AreEqual(mean2[0], target.Components[1].Mean[0], 1e-5);
            Assert.AreEqual(var2[0], target.Components[1].Variance[0], 1e-5);


            var expectedMean = Accord.Statistics.Tools.WeightedMean(values, weights);
            var expectedVar = Accord.Statistics.Tools.WeightedCovariance(values, weights);

            var actualMean = target.Mean;
            var actualVar = target.Covariance;

            Assert.AreEqual(expectedMean[0], actualMean[0], 0.0000001);
            Assert.AreEqual(expectedVar[0, 0], actualVar[0, 0], 0.68);
        }

        public void FittingOptionsTest()
        {
            // Create a degenerate problem
            double[][] sequences = new double[][] 
            {
                new double[] { 1,1,1,1,1,0,1,1,1,1 },
                new double[] { 1,1,1,1,0,1,1,1,1,1 },
                new double[] { 1,1,1,1,1,1,1,1,1,1 },
                new double[] { 1,1,1,1,1,1         },
                new double[] { 1,1,1,1,1,1,1       },
                new double[] { 1,1,1,1,1,1,1,1,1,1 },
                new double[] { 1,1,1,1,1,1,1,1,1,1 },
            };

            // Creates a continuous hidden Markov Model with two states organized in a ergodic
            //  topology and an underlying multivariate Normal distribution as density.
            var density = new MultivariateNormalDistribution(1);

            var model = new HiddenMarkovModel<MultivariateNormalDistribution>(new Ergodic(2), density);

            // Configure the learning algorithms to train the sequence classifier
            var teacher = new BaumWelchLearning<MultivariateNormalDistribution>(model)
            {
                Tolerance = 0.0001,
                Iterations = 0,

                // Configure options for fitting the normal distribution
                FittingOptions = new NormalOptions() { Regularization = 0.0001, }
            };

            // Fit the model. No exceptions will be thrown
            double logLikelihood = teacher.Run(sequences);
            double likelihood = Math.Exp(logLikelihood);

            Assert.AreEqual(47.434837528491286, logLikelihood, 1e-15);
            Assert.IsFalse(double.IsNaN(logLikelihood));

            Assert.AreEqual(0.0001, (teacher.FittingOptions as NormalOptions).Regularization);



            // Try without a regularization constant to get an exception
            bool thrown;

            thrown = false;
            density = new MultivariateNormalDistribution(1);
            model = new HiddenMarkovModel<MultivariateNormalDistribution>(new Ergodic(2), density);
            teacher = new BaumWelchLearning<MultivariateNormalDistribution>(model) { Tolerance = 0.0001, Iterations = 0, };
            Assert.IsNull(teacher.FittingOptions);
            try { teacher.Run(sequences); }
            catch { thrown = true; }
            Assert.IsTrue(thrown);

            thrown = false;
            density = new Accord.Statistics.Distributions.Multivariate.MultivariateNormalDistribution(1);
            model = new HiddenMarkovModel<MultivariateNormalDistribution>(new Ergodic(2), density);
            teacher = new BaumWelchLearning<MultivariateNormalDistribution>(model)
            {
                Tolerance = 0.0001,
                Iterations = 0,
                FittingOptions = new NormalOptions() { Regularization = 0 }
            };
            Assert.IsNotNull(teacher.FittingOptions);
            try { teacher.Run(sequences); }
            catch { thrown = true; }
            Assert.IsTrue(thrown);
        }

        public void FitTest2()
        {
            double[][] observations = 
            {
                new double[] { 1, 2 },
                new double[] { 1, 2 },
                new double[] { 1, 2 },
                new double[] { 1, 2 }
            };


            var target = new MultivariateNormalDistribution(2);

            bool thrown = false;
            try { target.Fit(observations); }
            catch (NonPositiveDefiniteMatrixException) { thrown = true; }

            Assert.IsTrue(thrown);

            NormalOptions options = new NormalOptions() { Regularization = double.Epsilon };

            // No exception thrown
            target.Fit(observations, options);
        }