smash::random::BesselSampler Class Reference

#include <random.h>

The intention of this class is to efficiently sample \( (N_1, N_2) \) from the Bessel distribution \( p(N_1,N_2) \sim \mathrm{Poi}(\nu_1) \mathrm{Poi}(\nu_2) \delta(N_1 - N_2 = N)\), where \(\mathrm{Poi}(\nu)\) denotes Poisson distribution with mean \(\nu\).

In other words, this class samples two Poisson numbers with a given mean and a fixed difference. The intended use is to sample the number of baryons and antibaryons, given their means and net baryon number.

The distribution of \( min(N_1,N_2) \) is a so-called Bessel distribution. Denoting \( a = \sqrt{\nu_1 \nu_2}\), \( p(N_{smaller} = k) = \frac{(a/2)^{2k+N}}{I_N(a) k! (N+k)!} \). We sample this distribution using the method suggested by Yuan and Kalbfleisch Yuan2000: if \( m = \frac{1}{2} (\sqrt{a^2 + N^2} - N) > 6\), then the distribution is approximated well by a Gaussian, else probabilities are computed explicitely and a table sampling is applied.

Definition at line 374 of file random.h.

Collaboration diagram for smash::random::BesselSampler:

Public Member Functions
	BesselSampler (const double poisson_mean1, const double poisson_mean2, const int fixed_difference)
	Construct a BesselSampler. More...
std::pair< int, int >	sample ()
	Sample two numbers from given Poissonians with a fixed difference. More...

Static Private Member Functions
static double	r_ (int n, double a)
	Compute the ratio of two Bessel functions r(n,a) = bessel_I(n+1,a)/bessel_I(n,a) using the continued fraction representation (see Yuan2000). More...

Private Attributes
random::discrete_dist< double >	dist_
	Vector to store tabulated values of probabilities for small m case (m <6). More...
double	m_
	Mode of the Bessel function, see Yuan2000 for details. More...
const double	a_
	Second parameter of Bessel distribution, see Yuan2000 for details. More...
const int	N_
	First parameter of Bessel distribution (= \( \nu \) in Yuan2000). More...
const bool	N_is_positive_
	Boolean variable to verify that N > 0. More...
double	mu_
	Mean of the Bessel distribution. More...
double	sigma_
	Standard deviation of the Bessel distribution. More...

Static Private Attributes
static constexpr double	m_switch_method_ = 6.0
	Switching mode to normal approximation. More...
static constexpr double	negligible_probability_ = 1.e-12
	Probabilities smaller than negligibly_probability are neglected. More...

Constructor & Destructor Documentation

BesselSampler()

smash::random::BesselSampler::BesselSampler	(	const double	poisson_mean1,
		const double	poisson_mean2,
		const int	fixed_difference
	)

Construct a BesselSampler.

Parameters

[in]	poisson_mean1	Mean of the first number's Poisson distribution.
[in]	poisson_mean2	Mean of the second number's Poisson distribution.
[in]	fixed_difference	Difference between the sampled numbers.

Returns

Constructed sampler.

Definition at line 17 of file random.cc.

    : a_(2.0 * std::sqrt(poisson_mean1 * poisson_mean2)),
      N_(std::abs(fixed_difference)),
      N_is_positive_(fixed_difference >= 0) {
  const auto &log = logger<LogArea::GrandcanThermalizer>();
  assert(poisson_mean1 >= 0.0);
  assert(poisson_mean2 >= 0.0);
  log.debug("Bessel sampler", ": Poisson mean N1 = ", poisson_mean1,
            ", Poisson mean N2 = ", poisson_mean2, ", N1 - N2 fixed to ",
            fixed_difference);
  m_ = 0.5 * (std::sqrt(a_ * a_ + N_ * N_) - N_);
  if (m_ >= m_switch_method_) {
    mu_ = 0.5 * a_ * r_(N_, a_);
    const double mean_sqr = mu_ * (1.0 + 0.5 * a_ * r_(N_ + 1, a_));
    sigma_ = std::sqrt(mean_sqr - mu_ * mu_);
    log.debug("m = ", m_, " -> using gaussian sampling with mean = ", mu_,
              ", sigma = ", sigma_);
  } else {
    log.debug("m = ", m_, " -> using direct sampling method");
    std::vector<double> probabilities;
    double wi = 1.0, sum = 0.0;
    int i = 0;
    do {
      sum += wi;
      probabilities.push_back(wi);
      wi *= 0.25 * a_ * a_ / (i + 1) / (N_ + i + 1);
      i++;
    } while (wi > negligible_probability_);
    i = 0;
    for (double p : probabilities) {
      p /= sum;
      log.debug("Probability (", i, ") = ", p);
      i++;
    }
    dist_.reset_weights(probabilities);
  }
}

Here is the call graph for this function:

Member Function Documentation

sample()

std::pair< int, int > smash::random::BesselSampler::sample

(

)

Sample two numbers from given Poissonians with a fixed difference.

Returns

Pair of first and second sampled number.

Definition at line 57 of file random.cc.

                                              {
  const int N_smaller = (m_ >= m_switch_method_)
                            ? std::round(random::normal(mu_, sigma_))
                            : dist_();
  return N_is_positive_ ? std::make_pair(N_smaller + N_, N_smaller)
                        : std::make_pair(N_smaller, N_smaller + N_);
}

Here is the call graph for this function:

Here is the caller graph for this function:

r_()

double smash::random::BesselSampler::r_ ( int  n, double  a  )

staticprivate

Compute the ratio of two Bessel functions r(n,a) = bessel_I(n+1,a)/bessel_I(n,a) using the continued fraction representation (see Yuan2000).

Parameters

[in]	n	First Bessel parameter.
[in]	a	Second Bessel parameter.

Returns

Ratio bessel_I(n+1,a)/bessel_I(n,a).

Definition at line 65 of file random.cc.

                                              {
  const double a_inv = 1.0 / a;
  double res = 0.0;
  // |x - continued fraction of order n| < 2^(-n+1), see the book
  // "Continued fractions" by Khinchin. For 10^-16 = ~2^-50 precision
  // 50 iterations should be sufficient. However, I found that for some
  // numerical reason at least 100 terms are needed.
  int i = 200;
  for (; i > 0; i--) {
    res = 1.0 / (a_inv * 2 * (n + i) + res);
  }
  // Check the known property of r(n,a) function, see iref{Yuan2000}.
  assert(a / (std::sqrt(a * a + (n + 1) * (n + 1)) + n + 1) <= res);
  assert(res <= a / (std::sqrt(a * a + n * n) + n));
  return res;
}

Here is the caller graph for this function:

Member Data Documentation

dist_

random::discrete_dist<double> smash::random::BesselSampler::dist_

private

Vector to store tabulated values of probabilities for small m case (m <6).

Definition at line 407 of file random.h.

m_

double smash::random::BesselSampler::m_

private

Mode of the Bessel function, see Yuan2000 for details.

Definition at line 410 of file random.h.

a_

const double smash::random::BesselSampler::a_

private

Second parameter of Bessel distribution, see Yuan2000 for details.

Definition at line 413 of file random.h.

N_

const int smash::random::BesselSampler::N_

private

First parameter of Bessel distribution (= \( \nu \) in Yuan2000).

Definition at line 416 of file random.h.

N_is_positive_

const bool smash::random::BesselSampler::N_is_positive_

private

Boolean variable to verify that N > 0.

Definition at line 419 of file random.h.

m_switch_method_

constexpr double smash::random::BesselSampler::m_switch_method_ = 6.0

staticprivate

Switching mode to normal approximation.

Note

Normal approximation of Bessel functions is possible for modes >= 6. See Yuan2000 for details.

Definition at line 426 of file random.h.

negligible_probability_

constexpr double smash::random::BesselSampler::negligible_probability_ = 1.e-12

staticprivate

Probabilities smaller than negligibly_probability are neglected.

Definition at line 429 of file random.h.

mu_

double smash::random::BesselSampler::mu_

private

Mean of the Bessel distribution.

Definition at line 432 of file random.h.

sigma_

double smash::random::BesselSampler::sigma_

private

Standard deviation of the Bessel distribution.

Definition at line 435 of file random.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/random.h
• src/random.cc