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