public double Gamma(double alpha, double beta)
        {
            int d = (int) alpha;
            double delta = alpha - d;

            double dlr = 0;

            if (delta > 0)
            {
                double r = Math.E / (Math.E + delta);
                double n1 = 0, n2 = 1;

                do
                {

                    double alpha1 = rand.NextDouble();
                    double alpha2 = rand.NextDouble();



                    if (alpha1 < r)
                    {
                        n1 = Math.Pow(alpha1 / r, 1.0 / delta);
                        n2 = alpha2 * Math.Pow(n1, delta - 1);
                    }
                    else
                    {
                        n1 = 1 - Math.Log((alpha1 - r) / (1 - r));
                        n2 = alpha2 * Math.Exp(-n1);
                    }
                } while (n2 > Math.Pow(n1, delta-1) * Math.Exp(-n1));

                dlr = n1;
            }

            double dr = 0;
            for (int i = 0; i < d; i++)
            {
                dr += Exponential(1);
            }



            return (dlr + dr) / beta;
        }