实例化新的概率分布

假设在我们的应用程序中，我们需要创建一个具有特定均值和标准差的正态（高斯）分布。为此，我们首先将框架的命名空间导入到我们的源文件中

using Accord.Statistics;
using Accord.Statistics.Distributions.Univariate;

之后，让我们创建一个均值μ=4，标准差σ=4.2的正态分布。

// Create a new Normal distribution with mu: 4 and s: 4.2
var normal = new NormalDistribution(mean: 4, stdDev: 4.2);

现在我们可以根据需要查询分布的任何属性

double mean   = normal.Mean;              // 4.0   (mu)
double median = normal.Median;            // 4.0   (for the Gaussian, same as the mean)
double mode   = normal.Mode;              // 4.0   (for the Gaussian, same as the mean)
double stdDev = normal.StandardDeviation; // 4.2   (sigma) 
double var    = normal.Variance;          // 17.64 (sigma squared)

我们还可以查询其相关函数

// The cumulative distribution function, or the probability that a real-valued 
// random variable will be found to have a value less than or equal to some x:
double cdf = normal.DistributionFunction(x: 1.4);        // 0.26794249453351904

// The probability density function, or the relative likelihood for a real-valued 
// random variable will be found to take on a given specific value of x:
double pdf = normal.ProbabilityDensityFunction(x: 1.4);  // 0.078423391448155175

// The log of the probability density function, useful for applications where
// precision is critical
double lpdf = normal.LogProbabilityDensityFunction(x: 1.4); // -2.5456330358182586

// The complementary distribution function, or the tail function, that gives the
// probability that a real-valued random variable will be found to have a value 
// greater than some x. This function is also known as the Survival function.
double ccdf = normal.ComplementaryDistributionFunction(x: 1.4); // 0.732057505466481

// The inverse distribution function, or the Quantile function, that is able to
// revert probability values back to the real value that produces that probability
double icdf = normal.InverseDistributionFunction(p: cdf); // 1.4

// The Hazard function, or the failure rate, the event rate at time t conditional 
// on survival until time t or later. Note that this function may only make sense
// when using time-defined distributions, such as the Poisson.
double hf = normal.HazardFunction(x: 1.4);            // 0.10712736480747137

// The cumulative hazard function, that gives how much the hazard 
// function accumulated over time until a given time instant x.
double chf = normal.CumulativeHazardFunction(x: 1.4); // 0.31189620872601354

// Every distribution has a friendly string representation
string str = normal.ToString(); // N(x; mu = 4, s2 = 17.64)

请注意，这些函数和度量是通用IDistribution接口的一部分。因此，无论我们创建什么分布，用法始终相同。例如，让我们将正态分布切换为离散的泊松分布

// Create a new Poisson distribution with lambda = 0.7
PoissonDistribution poisson = new PoissonDistribution(0.7);

double mean = poisson.Mean;                // 0.7    (lambda) 
double median = poisson.Median;            // 1.0
double mode = poisson.Mode;                // 0.7    (lambda)  
double stdDev = poisson.StandardDeviation; // 0.836  [sqrt((lambda))]
double var = poisson.Variance;             // 0.7    (lambda) 

// The cumulative distribution function, or the probability that a real-valued 
// random variable will be found to have a value less than or equal to some x:
double cdf = poisson.DistributionFunction(k: 1);        // 0.84419501644539618

// The probability density function, or the relative likelihood for a real-valued 
// random variable will be found to take on a given specific value of x:
double pdf = poisson.ProbabilityMassFunction(k: 1);  // 0.34760971265398666

// The log of the probability density function, useful for applications where
// precision is critical
double lpdf = poisson.LogProbabilityMassFunction(k: 1); // -1.0566749439387324

// The complementary distribution function, or the tail function, that gives the
// probability that a real-valued random variable will be found to have a value 
// greater than some x. This function is also known as the Survival function.
double ccdf = poisson.ComplementaryDistributionFunction(k: 1); // 0.15580498355460382

// The inverse distribution function, or the Quantile function, that is able to
// revert probability values back to the real value that produces that probability
int icdf = poisson.InverseDistributionFunction(p: cdf); // 1

// The Hazard function, or the failure rate, the event rate at time t conditional 
// on survival until time t or later. Note that this function may only make sense
// when using time-defined distributions, such as the Poisson.
double hf = poisson.HazardFunction(x: 1);            // 2.2310564445595058

// The cumulative hazard function, that gives how much the hazard 
// function accumulated over time until a given time instant x.
double chf = poisson.CumulativeHazardFunction(x: 1); // 1.8591501591854034

// Every distribution has a friendly string representation
string str = poisson.ToString(); // Poisson(x; lambda = 0.7)

我们也可以用同样的方法创建多维分布。在下一个例子中，我们来看看如何创建多维高斯分布并计算它的一些度量。

// 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[] mode = dist.Mode;         // { 4, 2 }
double[,] cov = dist.Covariance;   // { { 0.3, 0.1 }, { 0.1, 0.7 } }
double[] var = dist.Variance;      // { 0.3, 0.7 } (diagonal from cov)
int dimensions = dist.Dimension;   // 2

// Probability density functions
double pdf1 = dist.ProbabilityDensityFunction(2, 5);    // 0.000000018917884164743237
double pdf2 = dist.ProbabilityDensityFunction(4, 2);    // 0.35588127170858852
double pdf3 = dist.ProbabilityDensityFunction(3, 7);    // 0.000000000036520107734505265
double lpdf = dist.LogProbabilityDensityFunction(3, 7); // -24.033158110192296

// Cumulative distribution function (for up to two dimensions)
double cdf = dist.DistributionFunction(3, 5);               // 0.033944035782101534
double ccdf = dist.ComplementaryDistributionFunction(3, 5); // 0.00016755510356109232

从数据值估计分布

可能出现的情况是，我们不想使用均值和标准差的具体值，而是想直接从一组值中估计正态分布。例如，假设我们有以下值：

// Suppose those are our observed values
double[] observations = { 0.10, 0.40, 2.00, 2.00 };

// Now, in order to fit a distribution to those values, 
// all we have to do is to create a base distribution:
var normal = new NormalDistribution();

// and then pass those values to its .Fit method:
normal.Fit(observations);

// Result is N(x; 1.125, 1.017)
string text = normal.ToString() 

现在，假设我们的一些样本实际上比其他样本更重要。它们中的一些具有更大的*权重*。在这种情况下，我们可以为每个观测值指定一个实数值，表示每个特定观测值相对于其他观测值的权重。然后，我们可以将观测值和权重这两个向量传递给分布拟合过程，如下所示。

// Suppose those are our observed values and their individual weights
double[] observations = { 0.10, 0.40, 2.00, 2.00 };
double[] weights = { 0.25, 0.25, 0.25, 0.25 };

// We create a base distribution:
var normal = new NormalDistribution();

// and then pass the values and weights to its .Fit method:
normal.Fit(observations, weights);

// Result is N(x; 1.125, 1.017)
string text = normal.ToString() 

下面的示例展示了如何将一个多维高斯分布拟合到一组观测值。不幸的是，这些观测值是线性相关的，导致协方差矩阵是非正定的。为了解决这个问题，我们传递了一个拟合选项对象，指定应该以稳健的方式进行估计。

// Suppose we would like to fit a multidimensional Normal distribution to
// those samples. Those samples are linearly dependent between themselves,
// and attempting to fit them directly might result in an exception.
double[][] observations = 
{
    new double[] { 1, 2 },
    new double[] { 2, 4 },
    new double[] { 3, 6 },
    new double[] { 4, 8 }
};

// Create a base distribution with two dimensions
var normal = new MultivariateNormalDistribution(2);

// Fit the distribution to the data
normal.Fit(observations, new NormalOptions()
{
    Robust = true // do a robust estimation
});

从现有分布生成新样本

可以从框架中的*任何*单变量分布生成新样本。有些分布具有比其他分布更快的专门样本生成方法，但所有分布都可以进行采样，无论是通过专门方法还是使用默认的分位数函数方法。

以下示例展示了如何从均值为2，标准差为5的正态分布中生成1000000个样本。

// Create a Normal with mean 2 and sigma 5
var normal = new NormalDistribution(2, 5);

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

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

string result = actual.ToString("N2"); // N(x; mu = 2.01, s2 = 25.03)

以下示例展示了如何从试验次数为4、成功概率为0.2的二项分布中生成1000000个样本。

// Create a Binomial with n = 4 and p = 0.2
var binomial = new BinomialDistribution(4, 0.2);

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

// Try to estimate a new Binomial distribution from
// generated samples to check if they indeed match
var actual = new BinomialDistribution(4);
actual.Fit(samples);

string result = actual.ToString("N2"); // Binomial(x; n = 4.00, p = 0.20)

以下示例展示了如何从均值向量为[2, 6]，协方差矩阵为[2 1; 1 5]的多元正态分布中生成1000000个样本。

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

创建新的概率分布

在Accord.NET框架中，有两种创建新概率分布的方法：使用Lambda函数或使用继承。虽然后者是创建新分布的首选方法，因为它更容易在不同类和项目之间共享这些分布，但当我们需要尽快测试某些东西时，使用Lambda函数可能更容易编写快速而粗糙的代码。

使用Lambda函数

要使用lambda函数创建新分布，我们可以使用GeneralContinuousDistribution类。这个类可以基于定义所需分布的CDF、PDF或两者兼有的lambda函数来创建新分布。

  // Let's suppose we have a formula that defines a probability distribution
  // but we dont know much else about it. We don't know the form of its cumulative
  // distribution function, for example. We would then like to know more about
  // it, such as the underlying distribution's moments, characteristics, and 
  // properties.

  // Let's suppose the formula we have is this one:
  double mu = 5;
  double sigma = 4.2;

  Func<double, double> df = x =>
      1.0 / (sigma * Math.Sqrt(2 * Math.PI)) * Math.Exp(-Math.Pow(x - mu, 2) / (2 * sigma * sigma));

  // And for the moment, let's also pretend we don't know it is actually the
  // p.d.f. of a Gaussian distribution with mean 5 and std. deviation of 4.2.

  // So, let's create a distribution based _solely_ on the formula we have:
  var distribution = GeneralContinuousDistribution.FromDensityFunction(df);

  // Now, we can check everything that we can know about it:
  double mean = distribution.Mean;     // 5
  double median = distribution.Median; // 5
  double var = distribution.Variance;  // 17.64
  double mode = distribution.Mode;     // 5

  double cdf = distribution.DistributionFunction(x: 1.4); // 0.19568296915377595
  double pdf = distribution.ProbabilityDensityFunction(x: 1.4); // 0.065784567984404935
  double lpdf = distribution.LogProbabilityDensityFunction(x: 1.4); // -2.7213699972695058

  double ccdf = distribution.ComplementaryDistributionFunction(x: 1.4); // 0.80431703084622408
  double icdf = distribution.InverseDistributionFunction(p: cdf); // 1.3999999997024655

  double hf = distribution.HazardFunction(x: 1.4); // 0.081789351041333558
  double chf = distribution.CumulativeHazardFunction(x: 1.4); // 0.21776177055276186

使用继承

要使用继承创建新分布，首先需要知道您想创建哪种类型的分布。您需要创建一个新类，然后选择以下类之一进行继承：UnivariateContinuousDistribution、UnivariateDiscreteDistribution、MultivariateContinuousDistribution或MultivariateDiscreteDistribution。一旦您知道要选择哪个基类，请实现该类所需的所有虚方法。

例如，假设我们要实现维基百科中给出的Rademacher分布。

namespace MyApplicationNamespace
{
    class RademacherDistribution
    {
    }
}

根据维基百科的描述，它写道：“在概率论和统计学中，拉德马赫分布（以汉斯·拉德马赫命名）是一种离散的概率分布，其中随机变量X有50%的几率是+1或-1。^[1]”

由于这是一个离散概率分布，现在我们知道我们需要继承自UnivariateDiscreteDistribution。所以让我们这样做，右键单击类名，并实现其所有虚方法

    class RademacherDistribution : UnivariateDiscreteDistribution
    {
        public override double Mean
        {
            get { throw new NotImplementedException(); }
        }

        public override double Variance
        {
            get { throw new NotImplementedException(); }
        }

        public override double Entropy
        {
            get { throw new NotImplementedException(); }
        }

        public override AForge.IntRange Support
        {
            get { throw new NotImplementedException(); }
        }

        public override double DistributionFunction(int k)
        {
            throw new NotImplementedException();
        }

        public override double ProbabilityMassFunction(int k)
        {
            throw new NotImplementedException();
        }

        public override string ToString(string format, IFormatProvider formatProvider)
        {
            throw new NotImplementedException();
        }

        public override object Clone()
        {
            throw new NotImplementedException();
        }
    }

现在，我们所要做的就是填写那些缺失的字段

    class RademacherDistribution : UnivariateDiscreteDistribution
    {
        public override double Mean
        {
            get { return 0; }
        }

        public override double Variance
        {
            get { return 1; }
        }

        public override double Entropy
        {
            get { return System.Math.Log(2); }
        }

        public override AForge.IntRange Support
        {
            get { return new AForge.IntRange(-1, +1); }
        }

        public override double DistributionFunction(int k)
        {
            if (k < -1) return 0;
            if (k >= 1) return 1;
            return 0.5;
        }

        public override double ProbabilityMassFunction(int k)
        {
            if (k == -1) return 0.5;
            if (k == +1) return 0.5;
            return 0.0;
        }

        public override string ToString(string format, IFormatProvider formatProvider)
        {
            return "Rademacher(x)";
        }

        public override object Clone()
        {
            return new RademacherDistribution();
        }
    }

分布创建完成后，我们甚至可以查询我们未实现的功能，例如其逆分布函数、风险函数等。

var dist = new RademacherDistribution();

double median = dist.Median; // 0
double mode = dist.Mode;     // 0

double lpdf = dist.LogProbabilityMassFunction(k: -1);       // -0.69314718055994529
double ccdf = dist.ComplementaryDistributionFunction(k: 1); //  0.0

int icdf1 = dist.InverseDistributionFunction(p: 0); // -1
int icdf2 = dist.InverseDistributionFunction(p: 1); // +1

double hf = dist.HazardFunction(x: 0);            // 0.0
double chf = dist.CumulativeHazardFunction(x: 0); // 0.69314718055994529

string str = dist.ToString(); // Rademacher(x)

猜测数据的合理分布

// Create a new distribution, such as a Poisson
var poisson = new PoissonDistribution(lambda: 0.42);

// Draw enough samples from it
double[] samples = poisson.Generate(1000000).ToDouble();

// Let's pretend we don't know from which distribution
// those sample come from, and create an analysis object
// to check it for us:
var analysis = new DistributionAnalysis(samples);

// Compute the analysis
analysis.Compute();

// Get the most likely distribution (first)
var mostLikely = analysis.GoodnessOfFit[0];

// The result should be Poisson(x; lambda = 0.420961)
var result = mostLikely.Distribution.ToString();