Version: SMASH-3.1
Public Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
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 [63] : 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 377 of file random.h.

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 [63]). 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 [63] for details. More...
 
const double a_
 Second parameter of Bessel distribution, see [63] for details. More...
 
const int N_
 First parameter of Bessel distribution (= \( \nu \) in [63]). 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_mean1Mean of the first number's Poisson distribution.
[in]poisson_mean2Mean of the second number's Poisson distribution.
[in]fixed_differenceDifference between the sampled numbers.
Returns
Constructed sampler.

Definition at line 31 of file random.cc.

34  : a_(2.0 * std::sqrt(poisson_mean1 * poisson_mean2)),
35  N_(std::abs(fixed_difference)),
36  N_is_positive_(fixed_difference >= 0) {
37  assert(poisson_mean1 >= 0.0);
38  assert(poisson_mean2 >= 0.0);
39  logg[LGrandcanThermalizer].debug("Bessel sampler",
40  ": Poisson mean N1 = ", poisson_mean1,
41  ", Poisson mean N2 = ", poisson_mean2,
42  ", N1 - N2 fixed to ", fixed_difference);
43  m_ = 0.5 * (std::sqrt(a_ * a_ + N_ * N_) - N_);
44  if (m_ >= m_switch_method_) {
45  mu_ = 0.5 * a_ * r_(N_, a_);
46  const double mean_sqr = mu_ * (1.0 + 0.5 * a_ * r_(N_ + 1, a_));
47  sigma_ = std::sqrt(mean_sqr - mu_ * mu_);
48  logg[LGrandcanThermalizer].debug(
49  "m = ", m_, " -> using gaussian sampling with mean = ", mu_,
50  ", sigma = ", sigma_);
51  } else {
52  logg[LGrandcanThermalizer].debug("m = ", m_,
53  " -> using direct sampling method");
54  std::vector<double> probabilities;
55  double wi = 1.0, sum = 0.0;
56  int i = 0;
57  do {
58  sum += wi;
59  probabilities.push_back(wi);
60  wi *= 0.25 * a_ * a_ / (i + 1) / (N_ + i + 1);
61  i++;
62  } while (wi > negligible_probability_);
63  i = 0;
64  for (double p : probabilities) {
65  p /= sum;
66  logg[LGrandcanThermalizer].debug("Probability (", i, ") = ", p);
67  i++;
68  }
69  dist_.reset_weights(probabilities);
70  }
71 }
smash::random::BesselSampler::mu_
double mu_
Mean of the Bessel distribution.
Definition: random.h:435
smash::random::BesselSampler::sigma_
double sigma_
Standard deviation of the Bessel distribution.
Definition: random.h:438
smash::random::BesselSampler::N_is_positive_
const bool N_is_positive_
Boolean variable to verify that N > 0.
Definition: random.h:422
smash::random::BesselSampler::r_
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 ...
Definition: random.cc:81
smash::random::BesselSampler::negligible_probability_
static constexpr double negligible_probability_
Probabilities smaller than negligibly_probability are neglected.
Definition: random.h:432
smash::random::BesselSampler::m_
double m_
Mode of the Bessel function, see for details.
Definition: random.h:413
smash::random::BesselSampler::m_switch_method_
static constexpr double m_switch_method_
Switching mode to normal approximation.
Definition: random.h:429
smash::random::BesselSampler::dist_
random::discrete_dist< double > dist_
Vector to store tabulated values of probabilities for small m case (m <6).
Definition: random.h:410
smash::random::BesselSampler::N_
const int N_
First parameter of Bessel distribution (= in ).
Definition: random.h:419
smash::random::BesselSampler::a_
const double a_
Second parameter of Bessel distribution, see for details.
Definition: random.h:416
smash::random::discrete_dist::reset_weights
void reset_weights(const std::vector< T > &plist)
Reset the discrete distribution from a new probability list.
Definition: random.h:280
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
smash::pdg::p
constexpr int p
Proton.
Definition: pdgcode_constants.h:28
smash::LGrandcanThermalizer
static constexpr int LGrandcanThermalizer
Definition: grandcan_thermalizer.cc:20

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 73 of file random.cc.

73  {
74  const int N_smaller = (m_ >= m_switch_method_)
75  ? std::round(random::normal(mu_, sigma_))
76  : dist_();
77  return N_is_positive_ ? std::make_pair(N_smaller + N_, N_smaller)
78  : std::make_pair(N_smaller, N_smaller + N_);
79 }
smash::random::normal
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:250

◆ 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 [63]).

Parameters
[in]nFirst Bessel parameter.
[in]aSecond Bessel parameter.
Returns
Ratio bessel_I(n+1,a)/bessel_I(n,a).

Definition at line 81 of file random.cc.

81  {
82  const double a_inv = 1.0 / a;
83  double res = 0.0;
84  // |x - continued fraction of order n| < 2^(-n+1), see the book
85  // "Continued fractions" by Khinchin. For 10^-16 = ~2^-50 precision
86  // 50 iterations should be sufficient. However, I found that for some
87  // numerical reason at least 100 terms are needed.
88  int i = 200;
89  for (; i > 0; i--) {
90  res = 1.0 / (a_inv * 2 * (n + i) + res);
91  }
92  // Check the known property of r(n,a) function, see \cite Yuan2000.
93  assert(a / (std::sqrt(a * a + (n + 1) * (n + 1)) + n + 1) <= res);
94  assert(res <= a / (std::sqrt(a * a + n * n) + n));
95  return res;
96 }
smash::pdg::n
constexpr int n
Neutron.
Definition: pdgcode_constants.h:30

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 410 of file random.h.

◆ m_

double smash::random::BesselSampler::m_
private

Mode of the Bessel function, see [63] for details.

Definition at line 413 of file random.h.

◆ a_

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

Second parameter of Bessel distribution, see [63] for details.

Definition at line 416 of file random.h.

◆ N_

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

First parameter of Bessel distribution (= \( \nu \) in [63]).

Definition at line 419 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 422 of file random.h.

◆ m_switch_method_

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

Switching mode to normal approximation.

Note
Normal approximation of Bessel functions is possible for modes >= 6. See [63] for details.

Definition at line 429 of file random.h.

◆ negligible_probability_

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

Probabilities smaller than negligibly_probability are neglected.

Definition at line 432 of file random.h.

◆ mu_

double smash::random::BesselSampler::mu_
private

Mean of the Bessel distribution.

Definition at line 435 of file random.h.

◆ sigma_

double smash::random::BesselSampler::sigma_
private

Standard deviation of the Bessel distribution.

Definition at line 438 of file random.h.

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