C#, numerical calculation – calculation method and source program of Poisson Deviance (Poissondev)

1 Text format

using System;

namespace Legalsoft.Truffer
{
public class Poissondev : Ran
{
private double lambda { get; set; }
private double sqlam { get; set; }
private double loglam { get; set; }
private double lamexp { get; set; }
private double lambold { get; set; }
private double[] logfact { get; set; }

public Poissondev(double llambda, ulong ix) : base(ix)
{
this.lambda = llambda;
this.logfact = new double[1024];
for (int i = 0; i < 1024; i + + )
{
logfact[i] = -1.0;
}
this.lambold = -1.0;
}

public intdev()
{
double v2 = 0.0;
int k;
if (lambda < 5.0)
{
if (lambda != lambold)
{
lamexp = Math.Exp(-lambda);
}
k = -1;
double t = 1.0;
do
{
+ + k;
t *= doub();
} while (t > lamexp);
}
else
{
if (lambda != lambold)
{
sqlam = Math.Sqrt(lambda);
loglam = Math.Log(lambda);
}
for (; ; )
{
double u = 0.64 * doub();
double v = -0.68 + 1.28 * doub();
if (lambda > 13.5)
{
v2 = Globals.SQR(v);
if (v >= 0.0)
{
if (v2 > 6.5 * u * (0.64 – u) * (u + 0.2))
{
continue;
}
}
else
{
if (v2 > 9.6 * u * (0.66 – u) * (u + 0.07))
{
continue;
}
}
}
k = (int)Math.Floor(sqlam * (v / u) + lambda + 0.5);
if (k < 0)
{
continue;
}
double u2 = Globals.SQR(u);
if (lambda > 13.5)
{
if (v >= 0.0)
{
if (v2 < 15.2 * u2 * (0.61 - u) * (0.8 - u))
{
break;
}
}
else
{
if (v2 < 6.76 * u2 * (0.62 - u) * (1.4 - u))
{
break;
}
}
}
double lfac;
if (k < 1024)
{
if (logfact[k] < 0.0)
{
logfact[k] = Globals.gammln(k + 1.0);
}
lfac = logfact[k];
}
else
{
lfac = Globals.gammln(k + 1.0);
}
double p = sqlam * Math.Exp(-lambda + k * loglam – lfac);
if (u2 < p)
{
break;
}
}
}
lambold = lambda;
return k;
}

public int dev(double llambda)
{
lambda = llambda;
return dev();
}
}
}

2 Code format

using System;

namespace Legalsoft.Truffer
{
    public class Poissondev : Ran
    {
        private double lambda { get; set; }
        private double sqlam { get; set; }
        private double loglam { get; set; }
        private double lamexp { get; set; }
        private double lambold { get; set; }
        private double[] logfact { get; set; }

        public Poissondev(double llambda, ulong ix) : base(ix)
        {
            this.lambda = llambda;
            this.logfact = new double[1024];
            for (int i = 0; i < 1024; i + + )
            {
                logfact[i] = -1.0;
            }
            this.lambold = -1.0;
        }

        public intdev()
        {
            double v2 = 0.0;
            int k;
            if (lambda < 5.0)
            {
                if (lambda != lambold)
                {
                    lamexp = Math.Exp(-lambda);
                }
                k = -1;
                double t = 1.0;
                do
                {
                     + + k;
                    t *= doub();
                } while (t > lamexp);
            }
            else
            {
                if (lambda != lambold)
                {
                    sqlam = Math.Sqrt(lambda);
                    loglam = Math.Log(lambda);
                }
                for (; ; )
                {
                    double u = 0.64 * doub();
                    double v = -0.68 + 1.28 * doub();
                    if (lambda > 13.5)
                    {
                        v2 = Globals.SQR(v);
                        if (v >= 0.0)
                        {
                            if (v2 > 6.5 * u * (0.64 - u) * (u + 0.2))
                            {
                                continue;
                            }
                        }
                        else
                        {
                            if (v2 > 9.6 * u * (0.66 - u) * (u + 0.07))
                            {
                                continue;
                            }
                        }
                    }
                    k = (int)Math.Floor(sqlam * (v / u) + lambda + 0.5);
                    if (k < 0)
                    {
                        continue;
                    }
                    double u2 = Globals.SQR(u);
                    if (lambda > 13.5)
                    {
                        if (v >= 0.0)
                        {
                            if (v2 < 15.2 * u2 * (0.61 - u) * (0.8 - u))
                            {
                                break;
                            }
                        }
                        else
                        {
                            if (v2 < 6.76 * u2 * (0.62 - u) * (1.4 - u))
                            {
                                break;
                            }
                        }
                    }
                    double lfac;
                    if (k < 1024)
                    {
                        if (logfact[k] < 0.0)
                        {
                            logfact[k] = Globals.gammln(k + 1.0);
                        }
                        lfac = logfact[k];
                    }
                    else
                    {
                        lfac = Globals.gammln(k + 1.0);
                    }
                    double p = sqlam * Math.Exp(-lambda + k * loglam - lfac);
                    if (u2 < p)
                    {
                        break;
                    }
                }
            }
            lambold = lambda;
            return k;
        }

        public int dev(double llambda)
        {
            lambda = llambda;
            return dev();
        }
    }
}

The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Algorithm skill tree Home page Overview 54745 people are learning the system