多元正态分布代码实现

public sealed class MultivariateNormalDistribution : MultivariateContinuousDistribution
    {
        private Vector mean;

        private SymmetricMatrix covarianceMatrix;

        private Matrix factor;

        private double scale;

        private void DoPrecomputations()
        {
            if (this.factor != null)
            {
                return;
            }
            if (this.factor == null)
            {
                bool flag;
                flag = false;
                try
                {
                    this.factor = new CholeskyDecomposition(this.covarianceMatrix, overwrite: false).LowerTriangularFactor;
                    flag = this.factor.EstimateConditionNumber() < 10000.0;
                }
                catch (MatrixNotPositiveDefiniteException)
                {
                }
                catch (MatrixSingularException)
                {
                }
                if (!flag)
                {
                    SymmetricEigenvalueDecomposition symmetricEigenvalueDecomposition;
                    symmetricEigenvalueDecomposition = new SymmetricEigenvalueDecomposition(this.covarianceMatrix);
                    this.factor = symmetricEigenvalueDecomposition.Eigenvectors.ScaleColumns(Vector.Sqrt(symmetricEigenvalueDecomposition.Eigenvalues));
                }
            }
            this.scale = 0.0 - (0.5 * (double)this.mean.Length * 1.8378770664093456 + Math.Log(Math.Abs(this.factor.GetDeterminant())));
        }

        protected override void FillRandomVariateCore(System.Random random, Vector sample)
        {
            if (random == null)
            {
                ThrowException.ArgumentNull("random");
            }
            if (sample == null)
            {
                ThrowException.ArgumentNull("sample");
            }
            this.DoPrecomputations();
            NormalDistribution.Standard.GetRandomVariates(random, sample);
            sample.Multiply(this.factor, MatrixOperationSide.Left);
        }

        public static MultivariateNormalDistribution CreateStandard(int order)
        {
            return new MultivariateNormalDistribution(order);
        }

        public MultivariateNormalDistribution(int order)
            : this(new ConstantVector(order, 0.0))
        {
        }

        public MultivariateNormalDistribution(Vector mean)
            : base(Vector.SafeLength(mean))
        {
            if (mean == null)
            {
                ThrowException.ArgumentNull("mean");
            }
            if (mean.Length == 0)
            {
                ThrowException.ArgumentOutOfRange("mean");
            }
            this.mean = mean;
            this.covarianceMatrix = Matrix.CreateSymmetric(mean.Length);
            this.covarianceMatrix.GetDiagonal().SetValue(1.0);
        }

        public MultivariateNormalDistribution(Vector mean, SymmetricMatrix covarianceMatrix)
            : base(Vector.SafeLength(mean))
        {
            if (mean == null)
            {
                ThrowException.ArgumentNull("mean");
            }
            if (covarianceMatrix == null)
            {
                ThrowException.ArgumentNull("covarianceMatrix");
            }
            if (mean.Length == 0)
            {
                ThrowException.ArgumentOutOfRange("mean");
            }
            if (mean.Length != covarianceMatrix.RowCount)
            {
                throw new DimensionMismatchException(DimensionType.Row, "covarianceMatrix", DimensionType.Length, "mean");
            }
            this.mean = mean;
            this.covarianceMatrix = covarianceMatrix;
        }

        public MultivariateNormalDistribution(Matrix data)
            : base((!(data == null)) ? data.ColumnCount : 0)
        {
            if (data == null)
            {
                ThrowException.ArgumentNull("variable");
            }
            this.mean = data.ApplyToColumns(Stats.Mean);
            if (data.RowCount <= 1)
            {
                this.covarianceMatrix = Matrix.CreateSymmetric(data.ColumnCount);
                this.covarianceMatrix.GetDiagonal().SetValue(1.0);
            }
            else
            {
                this.covarianceMatrix = Stats.CovarianceMatrix(data);
            }
        }

        public MultivariateNormalDistribution(NumericalVariable[] variables)
            : base((variables != null) ? variables.Length : 0)
        {
            if (variables == null)
            {
                ThrowException.ArgumentNull("variable");
            }
            for (int i = 0; i < variables.Length; i++)
            {
                if (variables[i] == null)
                {
                    throw new ArgumentException();
                }
            }
            if (variables.Length == 0)
            {
                ThrowException.ArgumentOutOfRange("variables");
            }
            _ = variables[0].Length;
            this.mean = new GeneralVector(Array.ConvertAll(variables, Stats.Mean), reuseComponentArray: true);
            this.covarianceMatrix = Stats.CovarianceMatrix(variables);
        }

        public override double LogProbabilityDensityFunction(Vector x)
        {
            if (x == null)
            {
                ThrowException.ArgumentNull("x");
            }
            this.DoPrecomputations();
            GeneralVector generalVector;
            generalVector = (x - this.mean).AsGeneralVector();
            this.factor.Solve(generalVector, overwrite: true);
            return this.scale - 0.5 * generalVector.NormSquared();
        }

        public override Vector GetMeans()
        {
            return this.mean;
        }

        public override SymmetricMatrix GetVarianceCovarianceMatrix()
        {
            return this.covarianceMatrix;
        }
    }

 如果对您有帮忙,非常感谢您支持一下创造者的付出!

 感谢支持技术分享,请扫码点赞支持:

技术合作交流qq:2401315930

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

兴诚

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值