Using inverseCDF 
Returns a random variate from the $\textrm{Poisson}(\lambda)$ distribution.

Expects $\lambda > 0$.
@fullName rpois(lambda)()
@memberof poisson
516738
function rpois(lambda) {
      var mode;

      mode = { val: Math.floor(lambda) };
      mode.prob = Math.exp(-lambda + mode.val * Math.log(lambda) -
         lgamma(mode.val + 1));

      return inverseCDF(
         function getMode() { return mode; },
         function updateLeft() {
            if (this.val === 0) { return false; }
            this.prob *= this.val / lambda;
            this.val -= 1;
            return true;
         },
         function updateRight() {
            this.val += 1;
            this.prob *= lambda / this.val;
            return true;
         }
      );
   }
1208864
def rhypergeometric(n, m, N, size=None):
    
    if n == 0:
        return np.zeros(size, dtype=int)
    elif n == N:
        out = np.empty(size, dtype=int)
        out.fill(m)
        return out
    return np.random.hypergeometric(n, N - n, m, size)
288193
function getRandomInt(min, max) {
    const minInt = Math.ceil(min);
    const maxInt = Math.floor(max);
    return Math.floor(Math.random() * (maxInt - minInt)) + minInt;
}
509654
public function gaussian($max, $lambda = 0.2) {
    $tau = $max;
    $u = 1.0 * $this->getNext(0, self::MAX) / self::MAX;
    $x = - log(1 - (1 - exp(- $lambda * $tau)) * $u) / $lambda;
    $ran = (int) $x;
    if ( $ran >= $max ) $ran=$max-1; // If $tao and $lambda are wack
    return $ran;
  }
432385
@Override
	public double getStandardDeviation(RandomVariable probabilities) {
		if(isDeterministic()) {
			return 0.0;
		}
		if(size() == 0) {
			return Double.NaN;
		}

		return Math.sqrt(getVariance(probabilities));
	}
726500
public SDVariable randomExponential(double lambda, SDVariable shape) {
        return new RandomExponential(sameDiff(), shape, lambda).outputVariable();
    }
504033
@Override
    public INDArray rand(int[] shape, float min, float max, org.nd4j.linalg.api.rng.Random rng) {
        //ensure shapes that wind up being scalar end up with the write shape
        if (shape.length == 1 && shape[0] == 0) {
            shape = new int[] {1, 1};
        }
        return Nd4j.getDistributions().createUniform(min, max).sample(shape);
    }
1114579
public static double multinomialGaussianPdf(double[] mean, double[][] covariance, double[] x)  {
        MultivariateNormalDistribution gaussian = new 
            MultivariateNormalDistribution(mean, covariance);
        return gaussian.density(x);
    }
365773
def R_r(self, r):
        
        phi = np.random.uniform(0, 2*np.pi)
        theta = np.random.uniform(0, np.pi)
        x = r * np.sin(theta) * np.cos(phi)
        y = r * np.sin(theta) * np.sin(phi)
        R = np.sqrt(x**2 + y**2)
        return R, x, y
504761
@Override
	public double getQuantile(double quantile, RandomVariable probabilities) {
		return ((RandomVariableAAD) getRandomVariableInterface()).getRandomVariableInterface().getQuantile(quantile, probabilities);
	}
225817
@Override
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if (p < 0.0 || p > 1.0) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();
        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();
        if (p == 1.0) {
            return upperBound;
        }

        final double mu = getNumericalMean();
        final double sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;
        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig) || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {
                lowerBound = -1.0;
                while (cumulativeProbability(lowerBound) >= p) {
                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {
                upperBound = 1.0;
                while (cumulativeProbability(upperBound) < p) {
                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {

            public double value(final double x) {
                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {
            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();
            if (x - dx >= getSupportLowerBound()) {
                double px = cumulativeProbability(x);
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x;
                    while (upperBound - lowerBound > dx) {
                        final double midPoint = 0.5 * (lowerBound + upperBound);
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }
                    return upperBound;
                }
            }
        }
        return x;
    }
554801
public double randInverseCDF() {
        final double a0 = 2.50662823884;
        final double a1 = -18.61500062529;
        final double a2 = 41.39119773534;
        final double a3 = -25.44106049637;
        final double b0 = -8.47351093090;
        final double b1 = 23.08336743743;
        final double b2 = -21.06224101826;
        final double b3 = 3.13082909833;
        final double c0 = 0.3374754822726147;
        final double c1 = 0.9761690190917186;
        final double c2 = 0.1607979714918209;
        final double c3 = 0.0276438810333863;
        final double c4 = 0.0038405729373609;
        final double c5 = 0.0003951896511919;
        final double c6 = 0.0000321767881768;
        final double c7 = 0.0000002888167364;
        final double c8 = 0.0000003960315187;

        double y, r, x;

        double u = Math.random();
        while (u == 0.0) {
            u = Math.random();
        }
        
        y = u - 0.5;

        if (Math.abs(y) < 0.42) {
            r = y * y;
            x = y * (((a3 * r + a2) * r + a1) * r + a0)
                    / ((((b3 * r + b2) * r + b1) * r + b0) * r + 1);

        } else {
            r = u;
            if (y > 0) {
                r = 1 - u;
            }
            r = Math.log(-Math.log(r));
            x = c0 + r * (c1 + r * (c2 + r * (c3 + r * (c4 + r * (c5 + r * (c6 + r * (c7 + r * c8)))))));
            if (y < 0) {
                x = -(x);
            }
        }
        
        return mu + sigma * x;
    }
580199
public static double cdf(double x, double mu, double sigma, double lambda) {
    if(x == Double.NEGATIVE_INFINITY) {
      return 0.;
    }
    if(x == Double.POSITIVE_INFINITY) {
      return 1.;
    }
    final double u = lambda * (x - mu);
    final double v = lambda * sigma;
    final double v2 = v * v;
    final double logphi = FastMath.log(NormalDistribution.cdf(u, v2, v));
    return NormalDistribution.cdf(u, 0., v) - FastMath.exp(-u + v2 * .5 + logphi);
  }
1178910
def _uniform_phi(M):
        
        return np.random.uniform(-np.pi, np.pi, M)
128386
public function inverse($p) {
        $p->x -= $this->x0;
        $p->y -= $this->y0;
        $x = $p->x / $this->a;
        $y = $p->y / $this->a;

        if( $this->sphere ) {
            $cosz = 0.0;
            #$rh;
            $sinz = 0.0;

            $rh = sqrt(($x * $x) + ($y * $y));
            $phi = $rh * 0.5;

            if ($phi > 1.0) {
                Proj4php::reportError('laea:Inv:DataError');
                return null;
            }

            $phi = 2. * asin( $phi );

            if ($this->mode == $this->OBLIQ || $this->mode == $this->EQUIT) {
                $sinz = sin($phi);
                $cosz = cos($phi);
            }
            switch( $this->mode ) {
                case $this->EQUIT:
                    $phi = (abs( $rh ) <= Common::EPSLN) ? 0. : asin( $y * $sinz / $rh );
                    $x *= $sinz;
                    $y = $cosz * $rh;
                    break;
                case $this->OBLIQ:
                    $phi = (abs( $rh ) <= Common::EPSLN) ? $this->phi0 : asin( $cosz * $this->sinph0 + $y * $sinz * $this->cosph0 / $rh );
                    $x *= $sinz * $this->cosph0;
                    $y = ($cosz - sin( $phi ) * $this->sinph0) * $rh;
                    break;
                case $this->N_POLE:
                    $y = -$y;
                    $phi = Common::HALF_PI - $phi;
                    break;
                case $this->S_POLE:
                    $phi -= Common::HALF_PI;
                    break;
            }

            $lam = ($y == 0.0 && ($this->mode == $this->EQUIT || $this->mode == $this->OBLIQ)) ? 0. : atan2($x, $y);
        } else {
            /*
            $cCe;
            $sCe;
            $q;
            $rho;
            */
            $ab = 0.0;

            switch($this->mode) {
                case $this->EQUIT:
                case $this->OBLIQ:
                    $x /= $this->dd;
                    $y *= $this->dd;
                    $rho = sqrt( $x * $x + $y * $y );

                    if ($rho < Common::EPSLN) {
                        $p->x = 0.;
                        $p->y = $this->phi0;
                        return $p;
                    }

                    $sCe = 2. * asin( .5 * $rho / $this->rq );
                    $cCe = cos( $sCe );
                    $x *= ($sCe = sin($sCe));

                    if ($this->mode == $this->OBLIQ) {
                        $ab = $cCe * $this->sinb1 + $y * $sCe * $this->cosb1 / $rho;
                        $q = $this->qp * $ab;
                        $y = $rho * $this->cosb1 * $cCe - $y * $this->sinb1 * $sCe;
                    } else {
                        $ab = $y * $sCe / $rho;
                        $q = $this->qp * $ab;
                        $y = $rho * $cCe;
                    }
                    break;
                case $this->N_POLE:
                    $y = -$y;
                case $this->S_POLE:
                    $q = ($x * $x + $y * $y);

                    if ( ! $q) {
                        $p->x = 0.;
                        $p->y = $this->phi0;
                        return $p;
                    }
                    /*
                      q = $this->qp - q;
                     */
                    $ab = 1.0 - $q / $this->qp;

                    if ($this->mode == $this->S_POLE) {
                        $ab = - $ab;
                    }
                    break;
            }
            $lam = atan2($x, $y);
            $phi = $this->authlat(asin($ab), $this->apa);
        }

        /*
          $Rh = sqrt($p->x *$p->x +$p->y * $p->y);
          $temp = Rh / (2.0 * $this->a);

          if (temp > 1) {
          Proj4php::reportError("laea:Inv:DataError");
          return null;
          }

          $z = 2.0 * Common::asinz(temp);
          $sin_z=sin(z);
          $cos_z=cos(z);

          $lon =$this->long0;
          if (abs(Rh) > Common::EPSLN) {
          $lat = Common::asinz($this->sin_lat_o * cos_z +$this-> cos_lat_o * sin_z *$p->y / Rh);
          $temp =abs($this->lat0) - Common::HALF_PI;
          if (abs(temp) > Common::EPSLN) {
          temp = cos_z -$this->sin_lat_o * sin(lat);
          if(temp!=0.0) lon=Common::adjust_lon($this->long0+atan2($p->x*sin_z*$this->cos_lat_o,temp*Rh));
          } else if ($this->lat0 < 0.0) {
          lon = Common::adjust_lon($this->long0 - atan2(-$p->x,$p->y));
          } else {
          lon = Common::adjust_lon($this->long0 + atan2($p->x, -$p->y));
          }
          } else {
          lat = $this->lat0;
          }
         */
        //return(OK);
        $p->x = Common::adjust_lon($this->long0 + $lam);
        $p->y = $phi;
        return $p;
    }
1042221
@Override
    public INDArray randn(long rows, long columns, org.nd4j.linalg.api.rng.Random r) {
        return randn(new long[] {rows, columns}, r);
    }
1368368
private int randomValue(IntVar x) {
        int i = rnd.nextInt(x.getDomainSize());
        DisposableValueIterator ite = x.getValueIterator(true);
        int pos = -1;
        try {
            while (i >= 0) {
                pos = ite.next();
                i--;
            }
        } finally {
            ite.dispose();
        }
        return pos;
    }