超几何分布代码实现

public class HypergeometricDistribution : DiscreteDistribution
    {
        private int tagged;

        private int untagged;

        private int samples;

        private double W;

        private double P;

        private double P1;

        private double P2;

        private double P3;

        private double XL;

        private double XR;

        private double LAMDL;

        private double LAMDR;

        private double A;

        public int TaggedPopulation => this.tagged;

        public int UntaggedPopulation => this.untagged;

        public int NumberOfSamples => this.samples;

        public override double Mean => (double)this.tagged * (double)this.samples / (double)(this.tagged + this.untagged);

        public override double Variance
        {
            get
            {
                double num;
                num = this.tagged + this.untagged;
                return (double)this.samples * (double)this.tagged * (double)this.untagged * (num - (double)this.samples) / (num * num * (num - 1.0));
            }
        }

        public override double Skewness => (double)(this.untagged - this.tagged) * (double)(this.untagged + this.tagged - 2 * this.samples) / (double)(this.tagged + this.untagged - 2) * Math.Sqrt((double)(this.tagged + this.untagged - 1) / ((double)this.tagged * (double)this.untagged * (double)this.samples * (double)(this.tagged + this.untagged - this.samples)));

        public override double Kurtosis
        {
            get
            {
                double num;
                num = this.tagged + this.untagged;
                return -3.0 + (-1.0 + num) * num * num * (num * (1.0 + num) - 6.0 * (num - (double)this.samples) * (double)this.samples + 3.0 * (double)this.tagged * (double)this.untagged * (num * num * (double)(-2 + this.samples) + 6.0 * (num - (double)this.samples) * (double)this.samples - num * (double)this.samples * (double)this.samples) / (num * num)) / ((double)this.tagged * (double)this.untagged * (-3.0 + num) * (-2.0 + num) * (num - (double)this.samples) * (double)this.samples);
            }
        }

        public static int GetRandomVariate(System.Random random, int taggedPopulation, int untaggedPopulation, int samples)
        {
            if (random == null)
            {
                ThrowException.ArgumentNull("random");
            }
            return new HypergeometricDistribution(taggedPopulation, untaggedPopulation, samples).GetRandomVariate(random);
        }

        public HypergeometricDistribution(int taggedPopulation, int untaggedPopulation, int samples)
        {
            if (taggedPopulation < 0)
            {
                ThrowException.ArgumentOutOfRange("taggedPopulation");
            }
            if (untaggedPopulation < 0)
            {
                ThrowException.ArgumentOutOfRange("untaggedPopulation");
            }
            if (samples < 0)
            {
                ThrowException.ArgumentOutOfRange("samples");
            }
            if (samples > taggedPopulation + untaggedPopulation)
            {
                ThrowException.ArgumentOutOfRange("samples");
            }
            this.tagged = taggedPopulation;
            this.untagged = untaggedPopulation;
            this.samples = samples;
        }

        public override double Probability(int n)
        {
            if (n > this.tagged)
            {
                return 0.0;
            }
            if (n > this.samples)
            {
                return 0.0;
            }
            if (this.samples > this.untagged && n + this.untagged < this.samples)
            {
                return 0.0;
            }
            return Math.Exp(Combinatorics.LogCombinations(this.tagged, n) + Combinatorics.LogCombinations(this.untagged, this.samples - n) - Combinatorics.LogCombinations(this.tagged + this.untagged, this.samples));
        }

        public override double DistributionFunction(int n)
        {
            if (n < 0)
            {
                return 0.0;
            }
            if (n >= this.samples)
            {
                return 1.0;
            }
            if ((double)n < ((double)this.samples + 1.0) * ((double)this.tagged + 1.0) / ((double)(this.tagged + this.untagged) + 2.0))
            {
                double num;
                num = Math.Exp(Combinatorics.LogCombinations(this.tagged, n) + Combinatorics.LogCombinations(this.untagged, this.samples - n) - Combinatorics.LogCombinations(this.tagged + this.untagged, this.samples));
                double num2;
                num2 = num;
                int num3;
                num3 = n;
                while (num / num2 > 2.2250738585072014E-16)
                {
                    num = num * ((double)num3 * (double)(num3 + this.untagged - this.samples)) / ((double)(this.tagged - num3 + 1) * (double)(this.samples + 1 - num3));
                    num2 += num;
                    num3--;
                }
                return num2;
            }
            n++;
            double num4;
            num4 = Math.Exp(Combinatorics.LogCombinations(this.tagged, n) + Combinatorics.LogCombinations(this.untagged, this.samples - n) - Combinatorics.LogCombinations(this.tagged + this.untagged, this.samples));
            double num5;
            num5 = num4;
            int num6;
            num6 = n + 1;
            while (num4 / num5 > 2.2250738585072014E-16)
            {
                num4 = num4 / ((double)num6 * (double)(num6 + this.untagged - this.samples)) * ((double)(this.tagged - num6 + 1) * (double)(this.samples + 1 - num6));
                num5 += num4;
                num6++;
            }
            return 1.0 - num5;
        }

        public override int GetRandomVariate(System.Random random)
        {
            if (random == null)
            {
                ThrowException.ArgumentNull("random");
            }
            int num;
            num = this.tagged;
            int num2;
            num2 = this.untagged;
            int num3;
            num3 = this.samples;
            if (num < 0 || num2 < 0 || num3 < 0 || num3 > num + num2)
            {
                return -1;
            }
            bool flag;
            flag = true;
            bool flag2;
            flag2 = true;
            bool flag3;
            flag3 = true;
            double num4;
            num4 = num + num2;
            int num5;
            int num6;
            if (num <= num2)
            {
                num5 = num;
                num6 = num2;
            }
            else
            {
                num5 = num2;
                num6 = num;
            }
            int num7;
            num7 = ((!((double)(num3 + num3) >= num4)) ? num3 : ((int)num4 - num3));
            int num8;
            num8 = (int)(((double)num7 + 1.0) * ((double)num5 + 1.0) / (num4 + 2.0));
            int num9;
            num9 = Math.Max(0, num7 - num6);
            int num10;
            num10 = Math.Min(num5, num7);
            if (num9 == num10)
            {
                return num10;
            }
            int num11;
            if (num8 - num9 < 10)
            {
                if (flag2 || flag3)
                {
                    if (num7 < num6)
                    {
                        this.W = Math.Exp(57.56462733 + Combinatorics.LogFactorial(num6) + Combinatorics.LogFactorial(num5 + num6 - num7) - Combinatorics.LogFactorial(num6 - num7) - Combinatorics.LogFactorial(num5 + num6));
                    }
                    else
                    {
                        this.W = Math.Exp(57.56462733 + Combinatorics.LogFactorial(num5) + Combinatorics.LogFactorial(num7) - Combinatorics.LogFactorial(num7 - num6) - Combinatorics.LogFactorial(num5 + num6));
                    }
                }
                while (true)
                {
                    IL_0159:
                    this.P = this.W;
                    num11 = num9;
                    double num12;
                    num12 = random.NextDouble() * 1E+25;
                    while (num12 > this.P)
                    {
                        num12 -= this.P;
                        this.P = this.P * (double)(num5 - num11) * (double)(num7 - num11);
                        num11++;
                        this.P = this.P / (double)num11 / (double)(num6 - num7 + num11);
                        if (num11 > num10)
                        {
                            goto IL_0159;
                        }
                    }
                    break;
                }
            }
            else
            {
                if (flag2 || flag3)
                {
                    double num13;
                    num13 = (double)(int)(1.5 * Math.Sqrt((num4 - (double)num7) * (double)num7 * (double)num5 * (double)num6 / (num4 - 1.0) / num4 / num4)) + 0.5;
                    this.XL = (double)num8 - num13 + 0.5;
                    this.XR = (double)num8 + num13 + 0.5;
                    this.A = Combinatorics.LogFactorial(num8) + Combinatorics.LogFactorial(num5 - num8) + Combinatorics.LogFactorial(num7 - num8) + Combinatorics.LogFactorial(num6 - num7 + num8);
                    double num14;
                    num14 = Math.Exp(this.A - Combinatorics.LogFactorial((int)this.XL) - Combinatorics.LogFactorial((int)((double)num5 - this.XL)) - Combinatorics.LogFactorial((int)((double)num7 - this.XL)) - Combinatorics.LogFactorial((int)((double)(num6 - num7) + this.XL)));
                    double num15;
                    num15 = Math.Exp(this.A - Combinatorics.LogFactorial((int)(this.XR - 1.0)) - Combinatorics.LogFactorial((int)((double)num5 - this.XR + 1.0)) - Combinatorics.LogFactorial((int)((double)num7 - this.XR + 1.0)) - Combinatorics.LogFactorial((int)((double)(num6 - num7) + this.XR - 1.0)));
                    this.LAMDL = 0.0 - Math.Log(this.XL * ((double)(num6 - num7) + this.XL) / ((double)num5 - this.XL + 1.0) / ((double)num7 - this.XL + 1.0));
                    this.LAMDR = 0.0 - Math.Log(((double)num5 - this.XR + 1.0) * ((double)num7 - this.XR + 1.0) / this.XR / ((double)(num6 - num7) + this.XR));
                    this.P1 = num13 + num13;
                    this.P2 = this.P1 + num14 / this.LAMDL;
                    this.P3 = this.P2 + num15 / this.LAMDR;
                }
                do
                {
                    double num16;
                    num16 = random.NextDouble() * this.P3;
                    double num17;
                    num17 = random.NextDouble();
                    if (num16 < this.P1)
                    {
                        num11 = (int)(this.XL + num16);
                        goto IL_04bd;
                    }
                    if (num16 <= this.P2)
                    {
                        num11 = (int)(this.XL + Math.Log(num17) / this.LAMDL);
                        if (num11 >= num9)
                        {
                            num17 = num17 * (num16 - this.P1) * this.LAMDL;
                            goto IL_04bd;
                        }
                    }
                    else
                    {
                        num11 = (int)(this.XR - Math.Log(num17) / this.LAMDR);
                        if (num11 <= num10)
                        {
                            num17 = num17 * (num16 - this.P2) * this.LAMDR;
                            goto IL_04bd;
                        }
                    }
                    continue;
                    IL_04bd:
                    if (num8 < 100 || num11 <= 50)
                    {
                        double num18;
                        num18 = 1.0;
                        if (num8 < num11)
                        {
                            for (double num19 = num8 + 1; num19 <= (double)num11; num19 += 1.0)
                            {
                                num18 = num18 * ((double)num5 - num19 + 1.0) * ((double)num7 - num19 + 1.0) / ((double)(num6 - num7) + num19) / num19;
                            }
                        }
                        else if (num8 > num11)
                        {
                            for (double num20 = num11 + 1; num20 <= (double)num8; num20 += 1.0)
                            {
                                num18 = num18 * num20 * ((double)(num6 - num7) + num20) / ((double)num5 - num20 + 1.0) / ((double)num7 - num20 + 1.0);
                            }
                        }
                        if (num17 <= num18)
                        {
                            flag = false;
                        }
                        continue;
                    }
                    double num21;
                    num21 = num11;
                    double num22;
                    num22 = num21 + 1.0;
                    double num23;
                    num23 = num21 - (double)num8;
                    double num24;
                    num24 = (double)num5 - num21 + 1.0;
                    double num25;
                    num25 = (double)num7 - num21 + 1.0;
                    double num26;
                    num26 = (double)(num6 - num7) + num22;
                    double num27;
                    num27 = (0.0 - num23) / num22;
                    double num28;
                    num28 = num23 / num24;
                    double num29;
                    num29 = num23 / num25;
                    double num30;
                    num30 = (0.0 - num23) / num26;
                    double num31;
                    num31 = num24 * num25 / (num22 * num26) - 1.0;
                    double num32;
                    num32 = 1.0;
                    if (num31 < 0.0)
                    {
                        num32 = 1.0 + num31;
                    }
                    double num33;
                    num33 = num31 * (1.0 + num31 * (-0.5 + num31 / 3.0));
                    double num34;
                    num34 = num33 - 0.25 * (num31 * num31) * (num31 * num31) / num32;
                    double num35;
                    num35 = (double)num8 + 0.5;
                    double num36;
                    num36 = (double)(num5 - num8) + 0.5;
                    double num37;
                    num37 = (double)(num7 - num8) + 0.5;
                    double num38;
                    num38 = (double)(num6 - num7) + num35;
                    double num39;
                    num39 = num21 * num33 - (double)num8 * num34 + 0.0034 + num35 * num27 * (1.0 + num27 * (-0.5 + num27 / 3.0)) + num36 * num28 * (1.0 + num28 * (-0.5 + num28 / 3.0)) + num37 * num29 * (1.0 + num29 * (-0.5 + num29 / 3.0)) + num38 * num30 * (1.0 + num30 * (-0.5 + num30 / 3.0));
                    double num40;
                    num40 = Math.Log(num17);
                    if (num40 > num39)
                    {
                        flag = true;
                        continue;
                    }
                    double num41;
                    num41 = num35 * (num27 * num27) * (num27 * num27);
                    if (num27 < 0.0)
                    {
                        num41 /= 1.0 + num27;
                    }
                    double num42;
                    num42 = num36 * (num28 * num28) * (num28 * num28);
                    if (num28 < 0.0)
                    {
                        num42 /= 1.0 + num28;
                    }
                    double num43;
                    num43 = num37 * (num29 * num29) * (num29 * num29);
                    if (num29 < 0.0)
                    {
                        num43 /= 1.0 + num29;
                    }
                    if (num40 < num39 - 0.25 * (num41 + num42 + num43 + num38 * (num30 * num30) * (num30 * num30)) + (num21 + (double)num8) * (num34 - num33) - 0.0078)
                    {
                        flag = false;
                        flag = ((!(num40 <= this.A - GammaFunctions.LogGamma(1 + num11) - GammaFunctions.LogGamma(1 + num5 - num11) - GammaFunctions.LogGamma(1 + num7 - num11) - GammaFunctions.LogGamma(1 + num6 - num7 + num11))) ? true : false);
                    }
                }
                while (flag);
            }
            if ((double)(num3 + num3) >= num4)
            {
                num11 = ((num <= num2) ? (num - num11) : (num3 - num2 + num11));
            }
            else if (num > num2)
            {
                num11 = num3 - num11;
            }
            return num11;
        }

        public Histogram GetExpectedHistogram(double numberOfSamples)
        {
            return base.GetExpectedHistogram(0, Math.Min(this.tagged, this.samples) + 1, numberOfSamples);
        }
    }

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

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

技术合作交流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、付费专栏及课程。

余额充值