prob_dens_func.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file prob_dens_func.h
24  \brief File for probability density functions
25 */
26 #ifndef O2SCL_PROB_DENS_FUNC_H
27 #define O2SCL_PROB_DENS_FUNC_H
28 
29 #include <gsl/gsl_rng.h>
30 #include <gsl/gsl_randist.h>
31 #include <gsl/gsl_cdf.h>
32 
33 // For solving quadratics in bivariate gaussians
34 #include <gsl/gsl_poly.h>
35 
36 #include <random>
37 
38 #include <boost/numeric/ublas/vector.hpp>
39 #include <boost/numeric/ublas/operation.hpp>
40 
41 #include <o2scl/hist.h>
42 #include <o2scl/rng_gsl.h>
43 #include <o2scl/search_vec.h>
44 #include <o2scl/cholesky.h>
45 #include <o2scl/lu.h>
46 #include <o2scl/vec_stats.h>
47 
48 #ifndef DOXYGEN_NO_O2NS
49 namespace o2scl {
50 #endif
51 
52  /** \brief A one-dimensional probability density function
53 
54  This class is experimental.
55 
56  \future Give functions for mean, median, mode, variance, etc?
57 
58  \comment
59  For now, there aren't any pure virtual functions,
60  since this causes problems in creating an
61  std::vector<prob_dens_func> object below (especially
62  with intel compilers)
63  \endcomment
64  */
66 
67  public:
68 
69  /// Sample from the specified density
70  virtual double operator()() const {
71  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
72  return 0.0;
73  }
74 
75  /// The normalized density
76  virtual double pdf(double x) const {
77  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
78  return 0.0;
79  }
80 
81  /// The log of the normalized density
82  virtual double log_pdf(double x) const {
83  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
84  return 0.0;
85  }
86 
87  /// The cumulative distribution function (from the lower tail)
88  virtual double cdf(double x) const {
89  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
90  return 0.0;
91  }
92 
93  /// The inverse cumulative distribution function
94  virtual double invert_cdf(double cdf) const {
95  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
96  return 0.0;
97  }
98 
99  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
100  virtual double entropy() const {
101  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
102  return 0.0;
103  }
104 
105  };
106 
107  /** \brief A one-dimensional Gaussian probability density
108 
109  The distribution
110  \f[
111  P(x)=\frac{1}{\sigma \sqrt{2 \pi}}
112  e^{-\frac{\left(x-x_0\right)^2}{2\sigma^2}}
113  \f]
114 
115  This class is experimental.
116  */
118 
119  protected:
120 
121  /** \brief Central value
122  */
123  double cent_;
124 
125  /** \brief Width parameter
126 
127  A value of -1 indicates it is yet unspecified.
128  */
129  double sigma_;
130 
131  /// Base GSL random number generator
133 
134  public:
135 
136  /** \brief Create a standard normal distribution
137  */
139  cent_=0.0;
140  sigma_=1.0;
141  }
142 
143  /** \brief Create a Gaussian distribution with width \c sigma
144 
145  The value of \c sigma must be larger than zero.
146  */
147  prob_dens_gaussian(double cent, double sigma) {
148  if (sigma<0.0) {
149  O2SCL_ERR2("Tried to create a Gaussian dist. with sigma",
150  "<0 in prob_dens_gaussian::prob_dens_gaussian().",
151  exc_einval);
152  }
153  cent_=cent;
154  sigma_=sigma;
155  }
156 
157  virtual ~prob_dens_gaussian() {
158  }
159 
160  /// Copy constructor
162  cent_=pdg.cent_;
163  sigma_=pdg.sigma_;
164  r=pdg.r;
165  }
166 
167  /// Copy constructor with operator=
169  // Check for self-assignment
170  if (this!=&pdg) {
171  cent_=pdg.cent_;
172  sigma_=pdg.sigma_;
173  r=pdg.r;
174  }
175  return *this;
176  }
177 
178  /// Set the seed
179  void set_seed(unsigned long int s) {
180  r.set_seed(s);
181  return;
182  }
183 
184  /// Set the center
185  void set_center(double cent) {
186  cent_=cent;
187  return;
188  }
189 
190  /// Set the Gaussian width (must be positive)
191  void set_sigma(double sigma) {
192  if (sigma<0.0) {
193  O2SCL_ERR2("Tried to set negative sigma ",
194  "in prob_dens_gaussian::prob_dens_gaussian().",
195  exc_einval);
196  }
197  sigma_=sigma;
198  return;
199  }
200 
201  /// Get the center
202  double mean() {
203  return cent_;
204  }
205 
206  /// Get the Gaussian width
207  double stddev() {
208  if (sigma_<0.0) {
209  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
210  "get_sigma().",exc_einval);
211  }
212  return sigma_;
213  }
214 
215  /// Sample from the specified density
216  virtual double operator()() const {
217  if (sigma_<0.0) {
218  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
219  "operator().",exc_einval);
220  }
221  return cent_+gsl_ran_gaussian(&r,sigma_);
222  }
223 
224  /// The normalized density
225  virtual double pdf(double x) const {
226  if (sigma_<0.0) {
227  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
228  "pdf().",exc_einval);
229  }
230  return gsl_ran_gaussian_pdf(x-cent_,sigma_);
231  }
232 
233  /// The log of the normalized density
234  virtual double log_pdf(double x) const {
235  if (sigma_<0.0) {
236  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
237  "pdf().",exc_einval);
238  }
239  return log(gsl_ran_gaussian_pdf(x-cent_,sigma_));
240  }
241 
242  /// The cumulative distribution function (from the lower tail)
243  virtual double cdf(double x) const {
244  if (sigma_<0.0) {
245  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
246  "cdf().",exc_einval);
247  }
248  return gsl_cdf_gaussian_P(x-cent_,sigma_);
249  }
250 
251  /// The inverse cumulative distribution function
252  virtual double invert_cdf(double in_cdf) const {
253  if (sigma_<0.0) {
254  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
255  "invert_cdf().",exc_einval);
256  }
257  if (in_cdf<0.0 || in_cdf>1.0) {
258  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
259  "prob_dens_gaussian::invert_cdf().",exc_einval);
260  }
261  return gsl_cdf_gaussian_Pinv(in_cdf,sigma_)+cent_;
262  }
263 
264  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
265  virtual double entropy() const {
266  if (sigma_<0.0) {
267  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
268  "invert_cdf().",exc_einval);
269  }
270  return log(2.0*o2scl_const::pi*exp(1.0)*sigma_*sigma_);
271  }
272 
273  };
274 
275  /** \brief A one-dimensional probability density over
276  a finite range
277 
278  This class is experimental.
279  */
281 
282  public:
283 
284  /// Lower limit of the range
285  virtual double lower_limit() const=0;
286 
287  /// Uower limit of the range
288  virtual double upper_limit() const=0;
289 
290  };
291 
292  /** \brief A uniform one-dimensional probability density
293  over a finite range
294 
295  A flat distribution given by \f$ P(x)=1/(b-a) \f$ for \f$ a<x<b
296  \f$, where \f$ a \f$ is the lower limit and \f$ b \f$ is the
297  upper limit.
298 
299  This class is experimental.
300  */
302 
303  protected:
304 
305  /// Lower limit
306  double ll;
307 
308  /// Upper limit
309  double ul;
310 
311  /// The GSL random number generator
313 
314  public:
315 
316  /** \brief Create a blank uniform distribution
317  */
319  ll=1.0;
320  ul=0.0;
321  }
322 
323  /** \brief Create a uniform distribution from \f$ a<x<b \f$
324  */
325  prob_dens_uniform(double a, double b) {
326  // Ensure a<b
327  if (a>b) {
328  double tmp=a;
329  a=b;
330  b=tmp;
331  }
332  ll=a;
333  ul=b;
334  }
335 
336  virtual ~prob_dens_uniform() {
337  }
338 
339  /// Copy constructor
341  ll=pdg.ll;
342  ul=pdg.ul;
343  }
344 
345  /// Copy constructor with operator=
347  // Check for self-assignment
348  if (this==&pdg) return *this;
349  ll=pdg.ll;
350  ul=pdg.ul;
351  return *this;
352  }
353 
354  /// Set the seed
355  void set_seed(unsigned long int s) {
356  r.set_seed(s);
357  return;
358  }
359 
360  /** \brief Set the limits of the uniform distribution
361  */
362  void set_limits(double a, double b) {
363  // Ensure a<b
364  if (a>b) {
365  double tmp=a;
366  a=b;
367  b=tmp;
368  }
369  ll=a;
370  ul=b;
371  return;
372  }
373 
374  /// Lower limit of the range
375  virtual double lower_limit() const {
376  if (ll>ul) {
377  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
378  "lower_limit().",exc_einval);
379  }
380  return ll;
381  }
382 
383  /// Uower limit of the range
384  virtual double upper_limit() const {
385  if (ll>ul) {
386  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
387  "upper_limit().",exc_einval);
388  }
389  return ul;
390  }
391 
392  /// Operator from the specified density
393  virtual double operator()() const {
394  if (ll>ul) {
395  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
396  "operator().",exc_einval);
397  }
398  return gsl_ran_flat(&r,ll,ul);
399  }
400 
401  /// The normalized density
402  virtual double pdf(double x) const {
403  if (ll>ul) {
404  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
405  "pdf().",exc_einval);
406  }
407  if (x<ll || x>ul) return 0.0;
408  return gsl_ran_flat_pdf(x,ll,ul);
409  }
410 
411  /// The log of the normalized density
412  virtual double log_pdf(double x) const {
413  if (ll>ul) {
414  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
415  "pdf().",exc_einval);
416  }
417  if (x<ll || x>ul) return 0.0;
418  return log(gsl_ran_flat_pdf(x,ll,ul));
419  }
420 
421  /// The cumulative distribution function (from the lower tail)
422  virtual double cdf(double x) const {
423  if (ll>ul) {
424  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
425  "cdf().",exc_einval);
426  }
427  if (x<ll) return 0.0;
428  if (x>ul) return 1.0;
429  return gsl_cdf_flat_P(x,ll,ul);
430  }
431 
432  /// The inverse cumulative distribution function
433  virtual double invert_cdf(double in_cdf) const {
434  if (ll>ul) {
435  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
436  "invert_cdf().",exc_einval);
437  }
438  if (in_cdf<0.0 || in_cdf>1.0) {
439  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
440  "prob_dens_uniform::invert_cdf().",exc_einval);
441  }
442  return gsl_cdf_flat_Pinv(in_cdf,ll,ul);
443  }
444 
445  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
446  virtual double entropy() const {
447  if (ll>ul) {
448  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
449  "entropy().",exc_einval);
450  }
451  return log(ul-ll);
452  }
453 
454  };
455 
456  /** \brief A one-dimensional probability density over the
457  positive real numbers
458 
459  This class is experimental.
460  */
462 
463  };
464 
465  /** \brief Lognormal density function
466 
467  The distribution
468  \f[
469  P(x)=\frac{1}{x \sigma \sqrt{2 \pi}}
470  \exp \left[-\frac{\left(\ln x-\mu\right)^2}{2\sigma^2}\right]
471  \f]
472 
473  This class is experimental.
474  */
476 
477  protected:
478 
479  /** \brief Width parameter
480 
481  A value of -1 indicates it is yet unspecified.
482  */
483  double sigma_;
484 
485  /** \brief Central value
486 
487  A value of -1 indicates it is yet unspecified.
488  */
489  double mu_;
490 
491  /// The GSL random number generator
493 
494  public:
495 
496  /** \brief Create a blank lognormal distribution
497  */
499  sigma_=-1.0;
500  mu_=0.0;
501  }
502 
503  /** \brief Create lognormal distribution with mean parameter \c mu
504  and width parameter \c sigma
505 
506  The value of \c sigma must be larger than zero.
507  */
508  prob_dens_lognormal(double mu, double sigma) {
509  if (sigma<0.0) {
510  O2SCL_ERR2("Tried to create log normal dist. with mu or sigma",
511  "<0 in prob_dens_lognormal::prob_dens_lognormal().",
512  exc_einval);
513  }
514  mu_=mu;
515  sigma_=sigma;
516  }
517 
518  virtual ~prob_dens_lognormal() {
519  }
520 
521  /// Copy constructor
523  mu_=pdg.mu_;
524  sigma_=pdg.sigma_;
525  }
526 
527  /// Copy constructor with operator=
529  // Check for self-assignment
530  if (this==&pdg) return *this;
531  mu_=pdg.mu_;
532  sigma_=pdg.sigma_;
533  return *this;
534  }
535 
536  /** \brief Set the maximum and width of the lognormal distribution
537  */
538  void set_mu_sigma(double mu, double sigma) {
539  if (sigma<0.0) {
540  O2SCL_ERR2("Tried to set mu or sigma negative",
541  "in prob_dens_lognormal::prob_dens_lognormal().",
542  exc_einval);
543  }
544  mu_=mu;
545  sigma_=sigma;
546  }
547 
548  /// Set the seed
549  void set_seed(unsigned long int s) {
550  r.set_seed(s);
551  }
552 
553  /// Sample from the specified density
554  virtual double operator()() const {
555  return gsl_ran_lognormal(&r,mu_,sigma_);
556  }
557 
558  /// The normalized density
559  virtual double pdf(double x) const {
560  if (x<0.0) {
561  return 0.0;
562  }
563  return gsl_ran_lognormal_pdf(x,mu_,sigma_);
564  }
565 
566  /// The log of the normalized density
567  virtual double log_pdf(double x) const {
568  if (x<0.0) {
569  return 0.0;
570  }
571  return log(gsl_ran_lognormal_pdf(x,mu_,sigma_));
572  }
573 
574  /// The cumulative distribution function (from the lower tail)
575  virtual double cdf(double x) const {
576  if (x<0.0) {
577  return 0.0;
578  }
579  return gsl_cdf_lognormal_P(x,mu_,sigma_);
580  }
581 
582  /// The inverse cumulative distribution function
583  virtual double invert_cdf(double in_cdf) const {
584  if (in_cdf<0.0 || in_cdf>1.0) {
585  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
586  "prob_dens_lognormal::invert_cdf().",exc_einval);
587  }
588  return gsl_cdf_lognormal_Pinv(in_cdf,mu_,sigma_);
589  }
590 
591  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
592  virtual double entropy() const {
593  if (sigma_<0.0) {
594  O2SCL_ERR2("Parameters not set in prob_dens_lognormal::",
595  "entropy().",exc_einval);
596  }
597  return 0.5+0.5*log(2.0*o2scl_const::pi*sigma_*sigma_)+mu_;
598  }
599 
600  };
601 
602  /** \brief Probability density function based on a histogram
603 
604  \note This class is experimental.
605  */
607 
608  public:
609 
611 
612  protected:
613 
614  /// Search through the partial sums
616 
617  /// Number of original histogram bins
618  size_t n;
619 
620  /** \brief Normalized partial sum of histogram bins
621 
622  This vector has size \ref n plus one.
623  */
625 
626  /** \brief Vector specifying original histogram bins
627 
628  This vector has size \ref n plus one.
629  */
631 
632  /// Random number generator
633  mutable rng_gsl rng;
634 
635  public:
636 
637  prob_dens_hist();
638 
639  ~prob_dens_hist();
640 
641  /// Initialize with histogram \c h
642  void init(hist &h);
643 
644  /// Get reference to partial sums
645  const ubvector &partial_sums() { return sum; }
646 
647  /// Get reference to bin ranges
648  const ubvector &bin_ranges() { return range; }
649 
650  /// Generate a sample
651  virtual double operator()() const;
652 
653  /// Lower limit of the range
654  virtual double lower_limit() const;
655 
656  /// Uower limit of the range
657  virtual double upper_limit() const;
658 
659  /// The normalized density
660  virtual double pdf(double x) const;
661 
662  /// The log of the normalized density
663  virtual double log_pdf(double x) const;
664 
665  /// Cumulative distribution function (from the lower tail)
666  virtual double cdf(double x) const;
667 
668  /// Inverse cumulative distribution function (from the lower tail)
669  virtual double invert_cdf(double x) const;
670 
671  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
672  virtual double entropy() const {
673  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
674  return 0.0;
675  }
676 
677  };
678 
679  /** \brief A multi-dimensional probability density function
680 
681  This class is experimental.
682  */
683  template<class vec_t=boost::numeric::ublas::vector<double> >
685 
686  public:
687 
688  /// Return the dimensionality
689  virtual size_t dim() const {
690  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
691  return 0;
692  }
693 
694  /// The normalized density
695  virtual double pdf(const vec_t &x) const {
696  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
697  return 0.0;
698  }
699 
700  /// The log of the normalized density
701  virtual double log_pdf(const vec_t &x) const {
702  double val=pdf(x);
703  if (!std::isfinite(val) || val<0.0) {
704  O2SCL_ERR2("PDF not finite or negative in ",
705  "prob_dens_mdim::log_pdf().",o2scl::exc_efailed);
706  }
707  double val2=log(pdf(x));
708  if (!std::isfinite(val2)) {
709  std::cout << val << " " << val2 << std::endl;
710  O2SCL_ERR2("Log of PDF not finite in ",
711  "prob_dens_mdim::log_pdf().",o2scl::exc_efailed);
712  }
713  return val2;
714  }
715 
716  /// Sample the distribution
717  virtual void operator()(vec_t &x) const {
718  O2SCL_ERR("Executing blank parent function.",o2scl::exc_eunimpl);
719  return;
720  }
721 
722  };
723 
724  /** \brief A multidimensional distribution formed by the product
725  of several one-dimensional distributions
726  */
727  template<class vec_t=boost::numeric::ublas::vector<double> >
728  class prob_dens_mdim_factor : public prob_dens_mdim<vec_t> {
729 
730  protected:
731 
732  /// Vector of one-dimensional distributions
733  std::vector<prob_dens_func> list;
734 
735  public:
736 
737  prob_dens_mdim_factor(std::vector<prob_dens_func> &p_list) {
738  list=p_list;
739  }
740 
741  /// Copy constructor
743  list=pdmf.list;
744  }
745 
746  /// Copy constructor with operator=
747  prob_dens_mdim_factor &operator=
748  (const prob_dens_mdim_factor &pdmf) {
749  // Check for self-assignment
750  if (this!=&pdmf) {
751  list=pdmf.list;
752  }
753  return *this;
754  }
755 
756  /// Return the dimensionality
757  virtual size_t dim() const {
758  return list.size();
759  }
760 
761  /// The normalized density
762  virtual double pdf(const vec_t &x) const {
763  double ret=1.0;
764  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
765  return ret;
766  }
767 
768  /// The log of the normalized density
769  virtual double log_pdf(const vec_t &x) const {
770  double ret=1.0;
771  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
772  return log(ret);
773  }
774 
775  /// Sample the distribution
776  virtual void operator()(vec_t &x) const {
777  for(size_t i=0;i<list.size();i++) x[i]=list[i]();
778  return;
779  }
780 
781  };
782 
783  /** \brief A bivariate gaussian probability distribution
784 
785  For a two-dimensional gaussian, given a mean
786  \f$ ( \mu_x, \mu_y ) \f$ and a covariance matrix
787  \f[
788  \Sigma = \left(
789  \begin{array}{cc}
790  \sigma_x^2 & \rho \sigma_x \sigma_y \\
791  \rho \sigma_x \sigma_y & \sigma_y^2 \\
792  \end{array}
793  \right)
794  \f]
795  the PDF is
796  \f[
797  pdf(x,y) = \left(2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}\right)^{-1}
798  \exp \left\{ - \frac{1}{2 (1-\rho^2)}
799  \left[ \frac{(x-\mu_x)^2}{\sigma_x^2} +
800  \frac{(y-\mu_y)^2}{\sigma_y^2} -
801  \frac{2 \rho (x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \right]
802  \right\}
803  \f]
804  (taken from the Wikipedia page on the "Multivariate normal
805  distribution").
806 
807  The function \ref o2scl::prob_dens_mdim_biv_gaussian::contour()
808  gives a point on the contour line for a fixed value of the
809  PDF given an angle \f$ \theta \f$. In particular,
810  it solves
811  \f[
812  c = pdf(r \cos \theta + \mu_x, r \sin \theta + \mu_y )
813  \f]
814  for the radius \f$ r \f$ and then stores the values
815  \f$ r \cos \theta + \mu_x \f$ and \f$ r \sin \theta + \mu_y \f$
816  in the reference parameters named \c x and \c y . Thus
817  this function can be used to map out the full contour
818  by selecting values for \f$ \theta \in [0,2 \pi] \f$.
819 
820  The function \ref
821  o2scl::prob_dens_mdim_biv_gaussian::level_fixed_integral() gives
822  the value of the PDF for which the integral inside the
823  corresponding contour is some fraction of the total integral
824  (which is always 1). Given a fraction \f$ f \f$, the argument of
825  the exponential is related to the inverse of the cumulative
826  distribution function for the chi-squared probability
827  distribution for two degrees of freedom for \f$ 1-f \f$. For a
828  fraction \f$ f \f$, the value \f$ \chi^2 \f$ (i.e. the
829  Mahalanobis distance) is \f$ \chi^2 = -2 \log (1-f) \f$ and then
830  the value of the PDF for the corresponding contour is \f$
831  pdf(x,y) = \left(2 \pi \sigma_x \sigma_y
832  \sqrt{1-\rho^2}\right)^{-1} \exp (-\chi^2/2) \f$ .
833 
834  */
835  template<class vec_t=boost::numeric::ublas::vector<double> >
837 
838  private:
839 
840  /// The x coordinate of the centroid
841  double x0;
842 
843  /// The y coordinate of the centroid
844  double y0;
845 
846  /// The x standard deviation
847  double sig_x;
848 
849  /// The y standard deviation
850  double sig_y;
851 
852  /// The covariance
853  double rho;
854 
855  public:
856 
858  }
859 
860  /// Copy constructor
862  x0=pdmbg.x0;
863  y0=pdmbg.y0;
864  sig_x=pdmbg.sig_x;
865  sig_y=pdmbg.sig_y;
866  rho=pdmbg.rho;
867  }
868 
869  /// Copy constructor with operator=
870  prob_dens_mdim_biv_gaussian &operator=
872  // Check for self-assignment
873  if (this!=&pdmbg) {
874  x0=pdmbg.x0;
875  y0=pdmbg.y0;
876  sig_x=pdmbg.sig_x;
877  sig_y=pdmbg.sig_y;
878  rho=pdmbg.rho;
879  }
880  return *this;
881  }
882 
883  /** \brief Set the properties of the distribution
884 
885  \note If \f$ |\rho|\geq 1 \f$ this function will
886  call the error handler.
887  */
888  void set(double x_cent, double y_cent, double x_std, double y_std,
889  double covar) {
890  if (fabs(covar)>=1.0) {
891  O2SCL_ERR2("Covariance cannot have magnitude equal or larger than ",
892  "1 in prob_dens_mdim_biv_gaussian::set().",
894  }
895  x0=x_cent;
896  y0=y_cent;
897  sig_x=x_std;
898  sig_y=y_std;
899  rho=covar;
900  return;
901  }
902 
903  /** \brief Compute the normalized probability density
904  */
905  virtual double pdf(const vec_t &v) const {
906  double x=v[0], y=v[1];
907  double arg=-((x-x0)*(x-x0)/sig_x/sig_x+
908  (y-y0)*(y-y0)/sig_y/sig_y-
909  2.0*rho*(x-x0)*(y-y0)/sig_x/sig_y)/2.0/(1.0-rho*rho);
910  double ret=exp(arg)/2.0/o2scl_const::pi/sig_x/sig_y/
911  sqrt(1.0-rho*rho);
912  return ret;
913  }
914 
915  /** \brief Return the contour level corresponding to a fixed
916  integral
917  */
918  virtual double level_fixed_integral(double integral) {
919  // This comes from inverting the cumulative distribution function
920  // for the chi-squared distribution for two degrees of of freedom,
921  // i.e. exp(-x/2)
922  double arg=-2.0*log(1.0-integral);
923  // Now compute the pdf for the fixed value of the
924  // squared Mahalanobis distance
925  return exp(-0.5*arg)/2.0/o2scl_const::pi/sig_x/
926  sig_y/sqrt(1.0-rho*rho);
927  }
928 
929  /** \brief Return a point on the contour for a specified level
930  given an angle
931  */
932  virtual void contour(double level, double theta, vec_t &x) {
933  if (level<0.0) {
934  O2SCL_ERR2("Cannot produce contours for negative values in ",
935  "prob_dens_mdim_biv_gaussian::contour().",
937  }
938  double max=0.5/sig_x/sig_y/o2scl_const::pi/sqrt(1.0-rho*rho);
939  if (level>max) {
940  O2SCL_ERR2("Cannot produce contours larger than maximum in ",
941  "prob_dens_mdim_biv_gaussian::contour().",
943  }
944  double arg=-log(level*2.0*o2scl_const::pi*sig_x*sig_y*
945  sqrt(1.0-rho*rho))*2.0*(1.0-rho*rho);
946  double r2=arg/(cos(theta)*cos(theta)/sig_x/sig_x+
947  sin(theta)*sin(theta)/sig_y/sig_y-
948  2.0*rho/sig_x/sig_y*cos(theta)*sin(theta));
949  x[0]=sqrt(r2)*cos(theta)+x0;
950  x[1]=sqrt(r2)*sin(theta)+y0;
951  return;
952  }
953 
954  };
955 
956  /** \brief A multi-dimensional Gaussian probability density function
957  using a Cholesky decomposition
958 
959  Given a (square) covariance matrix, \f$ \Sigma \f$, and a mean
960  vector \f$ \mu \f$ the PDF is
961  \f[
962  P(x) = \det \left( 2 \pi \Sigma \right)^{-1/2}
963  \exp \left[ -\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right]
964  \f]
965 
966  Given the Cholesky decomposition \f$ A A^{T} = \Sigma \f$,
967  and a vector, \f$ z \f$ of samples from the standard Gaussian
968  with 0 mean and unit variance, one can create a sample
969  \f$ x \f$ from \f$ x = \mu + A z \f$ .
970 
971  \note This class inverts the matrix, necessary for computing the
972  pdf, but not for sampling the distribution, so for large
973  matrices the inversion can be a waste of computation if the pdf
974  is not needed.
975 
976  A separate class for the two-dimensional case is \ref
977  prob_dens_mdim_biv_gaussian .
978 
979  \note Const functions are not thread-safe because
980  mutable storage is used.
981 
982  \future Create alternate versions based on other
983  matrix decompositions?
984  */
985  template<class vec_t=boost::numeric::ublas::vector<double>,
986  class mat_t=boost::numeric::ublas::matrix<double> >
987  class prob_dens_mdim_gaussian : public prob_dens_mdim<vec_t> {
988 
989  protected:
990 
991  /// Cholesky decomposition
992  mat_t chol;
993 
994  /// Inverse of the covariance matrix
995  mat_t covar_inv;
996 
997  /// Location of the peak
998  vec_t peak;
999 
1000  /// Normalization factor, \f$ \det ( 2 \pi \Sigma)^{-1/2} \f$
1001  double norm;
1002 
1003  /// Number of dimensions
1004  size_t ndim;
1005 
1006  /// Temporary storage 1
1007  mutable vec_t q;
1008 
1009  /// Temporary storage 2
1010  mutable vec_t vtmp;
1011 
1012  public:
1013 
1014  /** \brief Standard normal
1015  \comment
1016  This has to be public so the user can set the random seed,
1017  or we have to create a new set_seed() function.
1018  \endcomment
1019  */
1021 
1022  /** \brief Get the Cholesky decomposition
1023  */
1024  const mat_t &get_chol() {
1025  return chol;
1026  }
1027 
1028  /** \brief Get the inverse of the covariance matrix
1029  */
1030  const mat_t &get_covar_inv() {
1031  return covar_inv;
1032  }
1033 
1034  /** \brief Get the peak location
1035  */
1036  const vec_t &get_peak() {
1037  return peak;
1038  }
1039 
1040  /** \brief Get the normalization
1041  */
1042  const double &get_norm() {
1043  return norm;
1044  }
1045 
1046  /// The dimensionality
1047  virtual size_t dim() const {
1048  return ndim;
1049  }
1050 
1051  /// Create an empty distribution
1053  ndim=0;
1054  }
1055 
1056  /// Copy constructor
1058  ndim=pdmg_loc.ndim;
1059  chol=pdmg_loc.chol;
1060  covar_inv=pdmg_loc.covar_inv;
1061  norm=pdmg_loc.norm;
1062  q.resize(ndim);
1063  vtmp.resize(ndim);
1064  }
1065 
1066  /// Copy constructor with operator=
1068  // Check for self-assignment
1069  if (this!=&pdmg_loc) {
1070  ndim=pdmg_loc.ndim;
1071  chol=pdmg_loc.chol;
1072  covar_inv=pdmg_loc.covar_inv;
1073  norm=pdmg_loc.norm;
1074  q.resize(ndim);
1075  vtmp.resize(ndim);
1076  }
1077  return *this;
1078  }
1079 
1080  /** \brief Create a distribution from a set of samples from a
1081  multidimensional Gaussian, returning the peak values and
1082  covariance matrix
1083 
1084  The matrix \c pts should have a size of \c n_pts in the first
1085  index and \c p_mdim in the second index
1086  */
1087  template<class mat2_t, class vec2_t,
1088  class mat2_col_t=const_matrix_column_gen<mat2_t> >
1089  int set(size_t p_mdim, size_t n_pts, const mat2_t &pts,
1090  const vec2_t &vals, vec_t &peak_arg, mat_t &covar_arg) {
1091 
1092  // Set peak with average and diagonal elements in covariance
1093  // matrix with variance
1094  for(size_t i=0;i<p_mdim;i++) {
1095  const mat2_col_t col(pts,i);
1096  peak_arg[i]=o2scl::wvector_mean<mat2_col_t>(n_pts,col,vals);
1097  // Square standard deviation
1098  covar_arg(i,i)=o2scl::wvector_stddev<mat2_col_t>(n_pts,col,vals);
1099  covar_arg(i,i)*=covar_arg(i,i);
1100  }
1101  // Setup off-diagonal covariance matrix
1102  for(size_t i=0;i<p_mdim;i++) {
1103  mat2_col_t col_i(pts,i);
1104  for(size_t j=i+1;j<p_mdim;j++) {
1105  const mat2_col_t col_j(pts,j);
1106  double cov=o2scl::wvector_covariance(n_pts,col_i,col_j,vals);
1107  covar_arg(i,j)=cov;
1108  covar_arg(j,i)=cov;
1109  }
1110  }
1111  set(p_mdim,peak_arg,covar_arg);
1112  return 0;
1113  }
1114 
1115  /** \brief Create a distribution from a set of samples from a
1116  multidimensional Gaussian
1117 
1118  The matrix \c pts should have a size of \c n_pts in the first
1119  index and \c p_mdim in the second index
1120  */
1121  template<class mat2_t, class vec2_t,
1122  class mat2_col_t=const_matrix_column_gen<mat2_t> >
1123  int set(size_t p_mdim, size_t n_pts, const mat2_t &pts,
1124  const vec2_t &vals) {
1125 
1126  vec_t peak_arg(p_mdim);
1127  mat_t covar_arg(p_mdim,p_mdim);
1128 
1129  set<mat2_t,vec2_t,mat2_col_t>(p_mdim,n_pts,pts,vals,peak_arg,covar_arg);
1130 
1131  return 0;
1132  }
1133  /** \brief Create a distribution from the covariance matrix
1134  */
1135  prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
1136  set(p_ndim,p_peak,covar);
1137  }
1138 
1139  /** \brief Set the peak and covariance matrix for the distribution
1140 
1141  \note This function is called in constructors and thus
1142  should not be virtual.
1143  */
1144  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
1145  if (p_ndim==0) {
1146  O2SCL_ERR("Zero dimension in prob_dens_mdim_gaussian::set().",
1148  }
1149  ndim=p_ndim;
1150  norm=1.0;
1151  peak.resize(ndim);
1152  for(size_t i=0;i<ndim;i++) peak[i]=p_peak[i];
1153  q.resize(ndim);
1154  vtmp.resize(ndim);
1155 
1156  // Perform the Cholesky decomposition of the covariance matrix
1157  chol=covar;
1159 
1160  // Find the inverse
1161  covar_inv=chol;
1162  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1163 
1164  // Force chol to be lower triangular and compute the square root
1165  // of the determinant
1166  double sqrt_det=1.0;
1167  for(size_t i=0;i<ndim;i++) {
1168  if (!std::isfinite(chol(i,i))) {
1169  O2SCL_ERR2("An entry of the Cholesky decomposition was not finite ",
1170  "in prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1171  }
1172  sqrt_det*=chol(i,i);
1173  for(size_t j=0;j<ndim;j++) {
1174  if (i<j) chol(i,j)=0.0;
1175  }
1176  }
1177 
1178  // Compute normalization
1179  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt_det;
1180  if (!std::isfinite(norm)) {
1181  O2SCL_ERR2("Normalization not finite in ",
1182  "prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1183  }
1184  }
1185 
1186  /** \brief Alternate set function for use when covariance matrix
1187  has already been decomposed and inverted
1188  */
1189  void set_alt(size_t p_ndim, vec_t &p_peak, mat_t &p_chol,
1190  mat_t &p_covar_inv, double p_norm) {
1191  ndim=p_ndim;
1192  peak=p_peak;
1193  chol=p_chol;
1194  covar_inv=p_covar_inv;
1195  norm=p_norm;
1196  q.resize(ndim);
1197  vtmp.resize(ndim);
1198  return;
1199  }
1200 
1201  /** \brief Given a data set and a covariance function, construct
1202  probability distribution based on a Gaussian process which
1203  includes noise
1204 
1205  \note The type <tt>mat_col_t</tt> is a matrix column type
1206  for the internal object matrix type <tt>mat_t</tt>, and
1207  not associated with the data type <tt>vec_vec_t</tt>.
1208  Since the default matrix type is
1209  <tt>boost::numeric::ublas::matrix &lt; double &gt; </tt>
1210  a good matrix column type for this function is
1211  <tt>boost::numeric::ublas::matrix_column &lt;
1212  boost::numeric::ublas::matrix &lt; double &gt; &gt;</tt> .
1213  This matrix column type is needed for the LU
1214  decomposition and inversion.
1215  */
1216  template<class vec_vec_t, class mat_col_t, class func_t>
1217  void set_gproc(size_t n_dim, size_t n_init,
1218  vec_vec_t &x, vec_t &y, func_t &fcovar) {
1219 
1220  // Construct the four covariance matrices
1221 
1222  mat_t KXsX(n_dim,n_init);
1223  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1224  for(size_t icol=0;icol<n_init;icol++) {
1225  KXsX(irow-n_init,icol)=fcovar(x[irow],x[icol]);
1226  }
1227  }
1228 
1229  mat_t KXXs=boost::numeric::ublas::trans(KXsX);
1230 
1231  mat_t KXX(n_init,n_init);
1232  for(size_t irow=0;irow<n_init;irow++) {
1233  for(size_t icol=0;icol<n_init;icol++) {
1234  if (irow>icol) {
1235  KXX(irow,icol)=KXX(icol,irow);
1236  } else {
1237  KXX(irow,icol)=fcovar(x[irow],x[icol]);
1238  }
1239  }
1240  }
1241 
1242  mat_t KXsXs(n_dim,n_dim);
1243  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1244  for(size_t icol=n_init;icol<n_dim+n_init;icol++) {
1245  if (irow>icol) {
1246  KXsXs(irow-n_init,icol-n_init)=KXsXs(icol-n_init,irow-n_init);
1247  } else {
1248  KXsXs(irow-n_init,icol-n_init)=fcovar(x[irow],x[icol]);
1249  }
1250  }
1251  }
1252 
1253  // Construct the inverse of KXX
1254  mat_t inv_KXX(n_init,n_init);
1256  int signum;
1257  o2scl_linalg::LU_decomp(n_init,KXX,p,signum);
1258  if (o2scl_linalg::diagonal_has_zero(n_dim,KXX)) {
1259  O2SCL_ERR2("KXX matrix is singular in ",
1260  "prob_dens_mdim_gaussian::set_gproc().",
1262  }
1263  o2scl_linalg::LU_invert<mat_t,mat_t,mat_col_t>(n_init,KXX,p,inv_KXX);
1264 
1265  // Compute the mean vector
1266  vec_t prod(n_init), mean(n_dim);
1267  boost::numeric::ublas::axpy_prod(inv_KXX,y,prod,true);
1268  boost::numeric::ublas::axpy_prod(KXsX,prod,mean,true);
1269 
1270  // Compute the covariance matrix
1271  mat_t covar(n_dim,n_dim), prod2(n_init,n_dim), prod3(n_dim,n_dim);
1272  boost::numeric::ublas::axpy_prod(inv_KXX,KXXs,prod2,true);
1273  boost::numeric::ublas::axpy_prod(KXsX,prod2,prod3,true);
1274  covar=KXsXs-prod3;
1275 
1276  // Now use set() in the parent class
1277  this->set(n_dim,mean,covar);
1278 
1279  }
1280 
1281  /// The normalized density
1282  virtual double pdf(const vec_t &x) const {
1283  if (ndim==0) {
1284  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1285  "pdf().",o2scl::exc_einval);
1286  }
1287  double ret=norm;
1288  for(size_t i=0;i<ndim;i++) {
1289  q[i]=x[i]-peak[i];
1290  }
1291  vtmp=prod(covar_inv,q);
1292  ret*=exp(-0.5*inner_prod(q,vtmp));
1293  return ret;
1294  }
1295 
1296  /// The log of the normalized density
1297  virtual double log_pdf(const vec_t &x) const {
1298  if (ndim==0) {
1299  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1300  "pdf().",o2scl::exc_einval);
1301  }
1302  double ret=log(norm);
1303  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
1304  vtmp=prod(covar_inv,q);
1305  ret+=-0.5*inner_prod(q,vtmp);
1306  return ret;
1307  }
1308 
1309  /// Sample the distribution
1310  virtual void operator()(vec_t &x) const {
1311  if (ndim==0) {
1312  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1313  "operator().",o2scl::exc_einval);
1314  }
1315  for(size_t i=0;i<ndim;i++) q[i]=pdg();
1316  vtmp=prod(chol,q);
1317  for(size_t i=0;i<ndim;i++) x[i]=peak[i]+vtmp[i];
1318  return;
1319  }
1320 
1321  };
1322 
1323  /** \brief Gaussian distribution bounded by a hypercube
1324 
1325  \note This class naively resamples the Gaussian until
1326  a sample is within bounds. This is a temporary hack and
1327  can be very slow depending on the size of the volume
1328  excluded.
1329 
1330  \warning The PDF is not yet properly normalized
1331  */
1332  template<class vec_t=boost::numeric::ublas::vector<double>,
1333  class mat_t=boost::numeric::ublas::matrix<double> >
1335  public prob_dens_mdim_gaussian<vec_t,mat_t> {
1336 
1337  protected:
1338 
1339  /** \brief Lower limits
1340  */
1341  vec_t low;
1342 
1343  /** \brief Upper limits
1344  */
1345  vec_t high;
1346 
1347  public:
1348 
1349  /** \brief Maximum number of samples
1350  */
1351  size_t samp_max;
1352 
1353  /** \brief Create an empty distribution
1354  */
1356  samp_max=100000;
1357  }
1358 
1359  /** \brief Create a distribution with the specified peak, covariance
1360  matrix, lower limits, and upper limits
1361  */
1362  prob_dens_mdim_bound_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar,
1363  vec_t &p_low, vec_t &p_high) {
1364  set(p_ndim,p_peak,covar,p_low,p_high);
1365  samp_max=100000;
1366  }
1367 
1368  /** \brief Set the peak, covariance matrix, lower limits, and upper
1369  limits
1370 
1371  \note This function is called in constructors and thus
1372  should not be virtual.
1373  */
1374  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar,
1375  vec_t &p_low, vec_t &p_high) {
1376  prob_dens_mdim_gaussian<vec_t,mat_t>::set(p_ndim,p_peak,covar);
1377  low=p_low;
1378  high=p_high;
1379  return;
1380  }
1381 
1382  /** \brief Compute the probability density function (arbitrary
1383  normalization)
1384  */
1385  virtual double pdf(const vec_t &x) const {
1386  for(size_t i=0;i<this->ndim;i++) {
1387  if (x[i]<low[i]) {
1388  O2SCL_ERR("Parameter too small in pdf().",
1390  }
1391  if (x[i]>high[i]) {
1392  O2SCL_ERR("Parameter too large in pdf().",
1394  }
1395  }
1397  }
1398 
1399  /** \brief Compute the natural log of the probability density function
1400  (arbitrary normalization)
1401  */
1402  virtual double log_pdf(const vec_t &x) const {
1403  for(size_t i=0;i<this->ndim;i++) {
1404  if (x[i]<low[i]) {
1405  O2SCL_ERR("Parameter too small in log_pdf().",
1407  }
1408  if (x[i]>high[i]) {
1409  O2SCL_ERR("Parameter too large in log_pdf().",
1411  }
1412  }
1414  }
1415 
1416  /** \brief Sample the distribution
1417  */
1418  virtual void operator()(vec_t &x) const {
1419  bool done=false;
1420  size_t j=0;
1421  while (done==false) {
1422  done=true;
1424  j++;
1425  for(size_t i=0;i<this->ndim;i++) {
1426  if (x[i]<low[i]) {
1427  done=false;
1428  //std::cout << "Too small " << i << " " << x[i] << " "
1429  //<< low[i] << std::endl;
1430  i=this->ndim;
1431  } else if (x[i]>high[i]) {
1432  done=false;
1433  //std::cout << "Too large " << i << " " << x[i] << " "
1434  //<< high[i] << std::endl;
1435  i=this->ndim;
1436  }
1437  }
1438  if (j>samp_max) {
1439  O2SCL_ERR2("Sampling failed in ",
1440  "prob_dens_mdim_bound_gaussian::operator().",
1442  }
1443  }
1444  return;
1445  }
1446 
1447  };
1448 
1449  /** \brief A multi-dimensional conditional probability density function
1450 
1451  Note that conditional probabilities are typically written \f$
1452  P(A|B) \f$, i.e. the probability of \f$ A \f$ given \f$ B \f$.
1453  \o2 arranges the function parameters for the functions \ref
1454  o2scl::prob_cond_mdim::pdf, \ref o2scl::prob_cond_mdim::log_pdf
1455  \ref o2scl::prob_cond_mdim::operator()(), so that \f$ B \f$ is
1456  given first, and \f$ A \f$ is second.
1457 
1458  \ref
1459  o2scl::prob_cond_mdim::log_metrop_hast is a vector from \f$ B
1460  \f$ as denoted above.
1461 
1462  This class is experimental.
1463  */
1464  template<class vec_t=boost::numeric::ublas::vector<double> >
1466 
1467  public:
1468 
1469  /// The dimensionality
1470  virtual size_t dim() const {
1471  return 0;
1472  }
1473 
1474  /** \brief The conditional probability of x_A given x_B,
1475  i.e. \f$ P(A|B) \f$
1476  */
1477  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1478 
1479  /** \brief The log of the conditional probability of x_A given x_B
1480  i.e. \f$ \log [P(A|B)] \f$
1481  */
1482  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1483 
1484  /// Sample the distribution
1485  virtual void operator()(const vec_t &x_B, vec_t &x_A) const=0;
1486 
1487  /** \brief Sample the distribution and return the
1488  log of the Metropolis-Hastings ratio
1489 
1490  The Metropolis-Hastings ratio for a step beginning at \f$ x \f$
1491  and ending at \f$ x^{\prime} \f$ is
1492  obeys
1493  \f[
1494  \frac{P(x^{\prime})g(x|x^{\prime})}{P(x)g(x^{\prime}|x)}
1495  \f]
1496  taking the log, this gives
1497  \f[
1498  \log[P(x^{\prime})] - \log[P(x)] +
1499  \log \left[ \frac{g(x|x^{\prime})}{g(x^{\prime}|x)} \right]
1500  \f]
1501  thus this function computes
1502  \f[
1503  \log \left[ g(x|x^{\prime}) \right]
1504  - \log \left[ g(x^{\prime}|x) \right]
1505  \f]
1506  and thus, to keep a similar notation to
1507  \ref prob_cond_mdim::pdf() where \f$ g(x^{\prime}|x) \f$
1508  is obtained from
1509  \code
1510  pdf(x,x_prime)
1511  \endcode
1512  this function computes
1513  \code
1514  h(x,x_prime) = log_pdf(x_prime,x)-log_pdf(x,x_prime);
1515  \endcode
1516 
1517  To check this, in the limit that \f$ g(x|x^{\prime})
1518  \rightarrow P(x) \f$ this function returns
1519  \f[
1520  \log \left[ \frac{P(x)}{P(x^{\prime})} \right]
1521  \f]
1522 
1523  */
1524  virtual double log_metrop_hast(const vec_t &x, vec_t &x_prime) const {
1525  operator()(x,x_prime);
1526  double val1=log_pdf(x_prime,x);
1527  double val2=log_pdf(x,x_prime);
1528  if (!std::isfinite(val1)) {
1529  std::cout << "val1: " << val1 << std::endl;
1530  O2SCL_ERR("Log pdf not finite in log_metrop_hast 1.",
1532  }
1533  if (!std::isfinite(val2)) {
1534  std::cout << "val2: " << val2 << std::endl;
1535  O2SCL_ERR("Log pdf not finite in log_metrop_hast 2.",
1537  }
1538  return val1-val2;
1539  }
1540 
1541  };
1542 
1543  /** \brief Conditional probability for a random walk inside a
1544  hypercube
1545 
1546  \comment
1547  I had previously used std::uniform_real_distribution
1548  instead of rng_gsl, but this caused problems with
1549  intel compilers.
1550  \endcomment
1551 
1552  This conditional probability is most useful in providing a
1553  Metropolis-Hastings distribution with a fixed step size which
1554  properly handles a boundary. The Metropolis-Hastings step is
1555  accepted if \f$ r \in [0,1] \f$ obeys
1556  \f[
1557  r < \frac{P(x^{\prime})g(x|x^{\prime})}
1558  {P(x)g(x^{\prime}|x)}
1559  \f]
1560  The function \f$ g(x^{\prime}|x) = g(x^{\prime},x)/g(x) \f$, and
1561  because of the fixed step size these probabilities are just
1562  proportional to the inverse allowed volume, i.e. \f$ V(x)^{-1}
1563  V^{-1}(x^{\prime}) / V(x)^{-1} = V^{-1}(x^{\prime}) \f$ . If \f$
1564  x^{\prime} \f$ is near a boundary then \f$ V(x^{\prime}) \f$ is
1565  decreased and the conditional probability increases accordingly.
1566  If the distance between \f$ x \f$ and \f$ x^{\prime} \f$ is
1567  unreachable in a step, then the PDF is zero.
1568 
1569  \note This class is experimental.
1570 
1571  \note Const functions are not thread-safe because the
1572  class contains an internal mutable random number generator
1573  object.
1574 
1575  \comment
1576  If we do not include the g ratio, then the edges
1577  will be undersampled because we don't properly include
1578  the rejections
1579 
1580  For \f$ g(x^{\prime}|x) \f$, if x is near the edge, then the
1581  cond. prob. is larger, thus the g ratio is smaller than 1,
1582  encouraging more rejections.
1583  \endcomment
1584  */
1585  template<class vec_t=boost::numeric::ublas::vector<double> >
1587 
1588  protected:
1589 
1590  /** \brief Step sizes
1591  */
1592  std::vector<double> u_step;
1593 
1594  /** \brief Lower limits
1595  */
1596  std::vector<double> u_low;
1597 
1598  /** \brief Upper limits
1599  */
1600  std::vector<double> u_high;
1601 
1602  /** \brief Internal random number generator
1603  */
1604  mutable rng_gsl rg;
1605 
1606  /** \brief Internal set function
1607 
1608  \comment
1609  This can't be virtual because it needs to be called
1610  by the constructor
1611  \endcomment
1612  */
1613  int set_internal(size_t sz, vec_t &step, vec_t &low, vec_t &high) {
1614  u_step.resize(sz);
1615  u_low.resize(sz);
1616  u_high.resize(sz);
1617  for(size_t i=0;i<sz;i++) {
1618  u_step[i]=step[i];
1619 
1620  if (!std::isfinite(low[i]) || !std::isfinite(high[i])) {
1621  O2SCL_ERR2("Limit not finite in prob_cond_mdim_fixed_step::",
1622  "set_internal().",o2scl::exc_einval);
1623  }
1624 
1625  // Force low and high to be properly ordered
1626  if (low[i]>high[i]) {
1627  u_high[i]=low[i];
1628  u_low[i]=high[i];
1629  } else {
1630  u_low[i]=low[i];
1631  u_high[i]=high[i];
1632  }
1633  }
1634  return 0;
1635  }
1636 
1637  public:
1638 
1640  }
1641 
1642  /// Copy constructor
1644  u_step=pcmfs.u_step;
1645  u_low=pcmfs.u_low;
1646  u_high=pcmfs.u_high;
1647  }
1648 
1649  /// Copy constructor with operator=
1651  {
1652  // Check for self-assignment
1653  if (this!=&pcmfs) {
1654  u_step=pcmfs.u_step;
1655  u_low=pcmfs.u_low;
1656  u_high=pcmfs.u_high;
1657  }
1658  return *this;
1659  }
1660 
1661  /** \brief Set the random number generator seed
1662  */
1663  void set_seed(unsigned long int s) {
1664  rg.set_seed(s);
1665  return;
1666  }
1667 
1668  /** \brief Create a conditional probability object
1669  with specified step sizes and limits
1670  */
1671  template<class=vec_t> prob_cond_mdim_fixed_step
1672  (vec_t &step, vec_t &low, vec_t &high) {
1673  if (step.size()!=low.size()) {
1674  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1675  "prob_cond_mdim_fixed_step constructor.",
1677  }
1678  if (step.size()!=high.size()) {
1679  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1680  "prob_cond_mdim_fixed_step constructor.",
1682  }
1683  set_internal(step.size(),step,low,high);
1684  }
1685 
1686  /** \brief Set step sizes and limits
1687  */
1688  virtual int set(vec_t &step, vec_t &low, vec_t &high) {
1689  if (step.size()!=low.size()) {
1690  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1691  "prob_cond_mdim_fixed_step::set().",
1693  }
1694  if (step.size()!=high.size()) {
1695  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1696  "prob_cond_mdim_fixed_step::set().",
1698  }
1699  set_internal(step.size(),step,low,high);
1700  return 0;
1701  }
1702 
1703  /// The dimensionality
1704  virtual size_t dim() const {
1705  return u_step.size();
1706  }
1707 
1708  /** \brief The conditional probability of x_A given x_B,
1709  i.e. \f$ P(A|B) \f$
1710  */
1711  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1712  double vol1=1.0;
1713  for(size_t i=0;i<u_step.size();i++) {
1714 
1715  if (fabs(x_A[i]-x_B[i])>u_step[i]) return 0.0;
1716 
1717  if (x_B[i]-u_step[i]<u_low[i]) {
1718  if (x_B[i]+u_step[i]>u_high[i]) {
1719  // If x+step is too large and x-step is too small
1720  vol1*=u_high[i]-u_low[i];
1721  } else {
1722  // If x-step is too small
1723  vol1*=x_B[i]+u_step[i]-u_low[i];
1724  }
1725  } else {
1726  if (x_B[i]+u_step[i]>u_high[i]) {
1727  // If x+step is too large
1728  vol1*=u_high[i]-(x_B[i]-u_step[i]);
1729  } else {
1730  // The normal case, where the volumes are both inside
1731  // of the boundaries
1732  vol1*=2.0*u_step[i];
1733  }
1734  }
1735  }
1736  return 1.0/vol1;
1737  }
1738 
1739  /** \brief The log of the conditional probability of x_A given x_B
1740  i.e. \f$ \log [P(A|B)] \f$
1741  */
1742  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1743  return log(pdf(x_B,x_A));
1744  }
1745 
1746  /// Sample the distribution
1747  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1748  size_t nv=u_step.size();
1749  for(size_t i=0;i<nv;i++) {
1750  if (x_B[i]<u_low[i] || x_B[i]>u_high[i]) {
1751  std::cout << "Input out of bounds in fixed_step::operator(): "
1752  << i << " " << x_B[i] << " "
1753  << u_low[i] << " " << u_high[i] << std::endl;
1754  O2SCL_ERR("Input out of bounds in fixed_step::operator().",
1756  }
1757  do {
1758  x_A[i]=x_B[i]+u_step[i]*(rg.random()*2.0-1.0);
1759  } while (x_A[i]<u_low[i] || x_A[i]>u_high[i]);
1760  }
1761  return;
1762  }
1763 
1764  };
1765 
1766  /** \brief A multi-dimensional conditional probability density function
1767  independent of the input
1768 
1769  The conditional probability, \f$ P(A|B) = P(A,B)/P(B) \f$. If
1770  the joint probability is factorizable because the events \f$ A
1771  \f$ and \f$ B \f$ are independent, i.e. \f$ P(A,B) = P(A) P(B)
1772  \f$, then \f$ P(A|B) = P(A) \f$ and is independent of \f$ B \f$.
1773  This class handles that particular case.
1774 
1775  This class is experimental.
1776  */
1777  template<class vec_t=boost::numeric::ublas::vector<double> >
1778  class prob_cond_mdim_indep : public prob_cond_mdim<vec_t> {
1779 
1780  protected:
1781 
1782  prob_dens_mdim<vec_t> &base;
1783 
1784  public:
1785 
1786  /** \brief Create a conditional probability distribution
1787  based on the specified probability distribution
1788  */
1790  }
1791 
1792  /// Copy constructor
1794  base=pcmi.base;
1795  }
1796 
1797  /// Copy constructor with operator=
1799  // Check for self-assignment
1800  if (this!=&pcmi) {
1801  base=pcmi.base;
1802  }
1803  return *this;
1804  }
1805 
1806  /// The dimensionality
1807  virtual size_t dim() const {
1808  return base.dim();
1809  }
1810 
1811  /** \brief The conditional probability of x_A given x_B,
1812  i.e. \f$ P(A|B) \f$
1813  */
1814  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1815  return base.pdf(x_A);
1816  }
1817 
1818  /** \brief The log of the conditional probability of x_A given x_B
1819  i.e. \f$ \log [P(A|B)] \f$
1820  */
1821  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1822  return base.log_pdf(x_A);
1823  }
1824 
1825  /// Sample the distribution
1826  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1827  return base(x_A);
1828  }
1829 
1830  };
1831 
1832  /** \brief A multi-dimensional Gaussian conditional probability
1833  density function
1834 
1835  This class is experimental.
1836 
1837  \todo This should be a symmetric conditional probability,
1838  i.e. \f$ P(x|y) = P(y|x) \f$. Test this.
1839 
1840  \note Const functions are not thread-safe because
1841  mutable storage is used.
1842  */
1843  template<class vec_t=boost::numeric::ublas::vector<double>,
1844  class mat_t=boost::numeric::ublas::matrix<double> >
1845  class prob_cond_mdim_gaussian : public prob_cond_mdim<vec_t> {
1846 
1847  protected:
1848 
1849  /// Cholesky decomposition
1850  mat_t chol;
1851 
1852  /// Inverse of the covariance matrix
1853  mat_t covar_inv;
1854 
1855  /// Normalization factor
1856  double norm;
1857 
1858  /// Number of dimensions
1859  size_t ndim;
1860 
1861  /// Temporary storage 1
1862  mutable vec_t q;
1863 
1864  /// Temporary storage 2
1865  mutable vec_t vtmp;
1866 
1867  /// Standard normal
1869 
1870  public:
1871 
1872  /** \brief Create an empty distribution
1873  */
1875  ndim=0;
1876  }
1877 
1878  /// Copy constructor
1880  ndim=pcmg.ndim;
1881  chol=pcmg.chol;
1882  covar_inv=pcmg.covar_inv;
1883  norm=pcmg.norm;
1884  q.resize(ndim);
1885  vtmp.resize(ndim);
1886  }
1887 
1888  /// Copy constructor with operator=
1890  // Check for self-assignment
1891  if (this!=&pcmg) {
1892  ndim=pcmg.ndim;
1893  chol=pcmg.chol;
1894  covar_inv=pcmg.covar_inv;
1895  norm=pcmg.norm;
1896  q.resize(ndim);
1897  vtmp.resize(ndim);
1898  }
1899  return *this;
1900  }
1901 
1902  /** \brief Create a distribution from the covariance matrix
1903  */
1904  prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar) {
1905  set(p_ndim,covar);
1906  }
1907 
1908  /// The dimensionality
1909  virtual size_t dim() const {
1910  return ndim;
1911  }
1912 
1913  /// Set the seed
1914  void set_seed(unsigned long int s) {
1915  pdg.set_seed(s);
1916  return;
1917  }
1918 
1919  /** \brief Set the covariance matrix for the distribution
1920  */
1921  void set(size_t p_ndim, mat_t &covar) {
1922  if (p_ndim==0) {
1923  O2SCL_ERR("Zero dimension in prob_cond_mdim_gaussian::set().",
1925  }
1926  ndim=p_ndim;
1927  norm=1.0;
1928  q.resize(ndim);
1929  vtmp.resize(ndim);
1930 
1931  // Perform the Cholesky decomposition
1932  chol=covar;
1934 
1935  // Find the inverse
1936  covar_inv=chol;
1937  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1938 
1939  // Force chol to be lower triangular and compute the determinant
1940  double sqrt_det=1.0;
1941  for(size_t i=0;i<ndim;i++) {
1942  if (!std::isfinite(chol(i,i))) {
1943  O2SCL_ERR2("An entry of the Cholesky decomposition was not finite ",
1944  "in prob_cond_mdim_gaussian::set().",o2scl::exc_einval);
1945  }
1946  sqrt_det*=chol(i,i);
1947  for(size_t j=0;j<ndim;j++) {
1948  if (i<j) chol(i,j)=0.0;
1949  }
1950  }
1951 
1952  // Compute normalization
1953  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt_det;
1954  if (!std::isfinite(norm)) {
1955  O2SCL_ERR2("Normalization not finite in ",
1956  "prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1957  }
1958  }
1959 
1960  /** \brief The conditional probability of x_A given x_B,
1961  i.e. \f$ P(A|B) \f$
1962  */
1963  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1964  if (ndim==0) {
1965  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1966  "pdf().",o2scl::exc_einval);
1967  }
1968  double ret=norm;
1969  for(size_t i=0;i<ndim;i++) q[i]=x_A[i]-x_B[i];
1970  vtmp=prod(covar_inv,q);
1971  ret*=exp(-0.5*inner_prod(q,vtmp));
1972  return ret;
1973  }
1974 
1975  /** \brief The log of the conditional probability of x_A given x_B
1976  i.e. \f$ \log [P(A|B)] \f$
1977  */
1978  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1979  if (ndim==0) {
1980  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1981  "pdf().",o2scl::exc_einval);
1982  }
1983  double ret=log(norm);
1984  for(size_t i=0;i<ndim;i++) q[i]=x_A[i]-x_B[i];
1985  vtmp=prod(covar_inv,q);
1986  ret+=-0.5*inner_prod(q,vtmp);
1987  /*
1988  std::cout << "pdmg lp: ";
1989  for (size_t i=0;i<ndim;i++) std::cout << x_A[i] << " ";
1990  for (size_t i=0;i<ndim;i++) std::cout << q[i] << " ";
1991  for (size_t i=0;i<ndim;i++) std::cout << vtmp[i] << " ";
1992  std::cout << ret << std::endl;
1993  */
1994  return ret;
1995  }
1996 
1997  /// Sample the distribution
1998  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1999  if (ndim==0) {
2000  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
2001  "operator().",o2scl::exc_einval);
2002  }
2003  for (size_t i=0;i<ndim;i++) q[i]=pdg();
2004  vtmp=prod(chol,q);
2005  for (size_t i=0;i<ndim;i++) x_A[i]=x_B[i]+vtmp[i];
2006  /*
2007  std::cout << "pdmg op: ";
2008  for (size_t i=0;i<ndim;i++) std::cout << x_A[i] << " ";
2009  std::cout << std::endl;
2010  */
2011  return;
2012  }
2013 
2014  };
2015 
2016 #ifdef O2SCL_NEVER_DEFINED
2017  /** \brief A multidimensional normal distribution from
2018  a Gaussian process
2019 
2020  \future The linear algebra only works with ublas and is
2021  not optimized.
2022  */
2023  template<class vec_t=boost::numeric::ublas::vector<double>,
2024  class mat_t=boost::numeric::ublas::matrix<double>,
2025  class mat_col_t=boost::numeric::ublas::matrix_column<mat_t> >
2026  class prob_dens_mdim_gproc :
2027  public o2scl::prob_dens_mdim_gaussian<vec_t> {
2028 
2029  public:
2030 
2031  };
2032 #endif
2033 
2034 #ifndef DOXYGEN_NO_O2NS
2035 }
2036 #endif
2037 
2038 #endif
o2scl::prob_dens_hist::log_pdf
virtual double log_pdf(double x) const
The log of the normalized density.
o2scl::prob_dens_mdim_bound_gaussian::samp_max
size_t samp_max
Maximum number of samples.
Definition: prob_dens_func.h:1351
o2scl::hist
A one-dimensional histogram class.
Definition: hist.h:113
o2scl::prob_dens_positive
A one-dimensional probability density over the positive real numbers.
Definition: prob_dens_func.h:461
o2scl::prob_dens_mdim_gaussian::get_chol
const mat_t & get_chol()
Get the Cholesky decomposition.
Definition: prob_dens_func.h:1024
o2scl::prob_dens_lognormal::entropy
virtual double entropy() const
Entropy of the distribution ( )
Definition: prob_dens_func.h:592
o2scl::prob_cond_mdim::dim
virtual size_t dim() const
The dimensionality.
Definition: prob_dens_func.h:1470
o2scl::prob_cond_mdim_gaussian::covar_inv
mat_t covar_inv
Inverse of the covariance matrix.
Definition: prob_dens_func.h:1853
o2scl::prob_dens_mdim_gaussian::prob_dens_mdim_gaussian
prob_dens_mdim_gaussian(const prob_dens_mdim_gaussian &pdmg_loc)
Copy constructor.
Definition: prob_dens_func.h:1057
o2scl::prob_dens_mdim_bound_gaussian::prob_dens_mdim_bound_gaussian
prob_dens_mdim_bound_gaussian()
Create an empty distribution.
Definition: prob_dens_func.h:1355
o2scl::prob_dens_lognormal::log_pdf
virtual double log_pdf(double x) const
The log of the normalized density.
Definition: prob_dens_func.h:567
o2scl::rng_gsl
Random number generator (GSL)
Definition: rng_gsl.h:55
o2scl::prob_cond_mdim_indep::log_pdf
virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const
The log of the conditional probability of x_A given x_B i.e. .
Definition: prob_dens_func.h:1821
o2scl::prob_cond_mdim_gaussian::dim
virtual size_t dim() const
The dimensionality.
Definition: prob_dens_func.h:1909
o2scl::prob_cond_mdim_fixed_step::set_internal
int set_internal(size_t sz, vec_t &step, vec_t &low, vec_t &high)
Internal set function.
Definition: prob_dens_func.h:1613
o2scl::prob_dens_hist::partial_sums
const ubvector & partial_sums()
Get reference to partial sums.
Definition: prob_dens_func.h:645
o2scl::prob_dens_lognormal::prob_dens_lognormal
prob_dens_lognormal()
Create a blank lognormal distribution.
Definition: prob_dens_func.h:498
o2scl::prob_cond_mdim_fixed_step::dim
virtual size_t dim() const
The dimensionality.
Definition: prob_dens_func.h:1704
o2scl::prob_cond_mdim::log_metrop_hast
virtual double log_metrop_hast(const vec_t &x, vec_t &x_prime) const
Sample the distribution and return the log of the Metropolis-Hastings ratio.
Definition: prob_dens_func.h:1524
o2scl::prob_cond_mdim_gaussian::vtmp
vec_t vtmp
Temporary storage 2.
Definition: prob_dens_func.h:1865
o2scl::prob_dens_gaussian::entropy
virtual double entropy() const
Entropy of the distribution ( )
Definition: prob_dens_func.h:265
o2scl::prob_dens_mdim_factor::operator()
virtual void operator()(vec_t &x) const
Sample the distribution.
Definition: prob_dens_func.h:776
o2scl::prob_dens_gaussian::pdf
virtual double pdf(double x) const
The normalized density.
Definition: prob_dens_func.h:225
o2scl::prob_cond_mdim_fixed_step::rg
rng_gsl rg
Internal random number generator.
Definition: prob_dens_func.h:1604
o2scl::prob_dens_gaussian::prob_dens_gaussian
prob_dens_gaussian()
Create a standard normal distribution.
Definition: prob_dens_func.h:138
o2scl::prob_dens_hist
Probability density function based on a histogram.
Definition: prob_dens_func.h:606
o2scl::permutation
A class for representing permutations.
Definition: permutation.h:70
o2scl::prob_dens_uniform::ll
double ll
Lower limit.
Definition: prob_dens_func.h:306
o2scl::prob_dens_hist::lower_limit
virtual double lower_limit() const
Lower limit of the range.
boost::numeric::ublas::vector< double >
o2scl::prob_cond_mdim_gaussian::set_seed
void set_seed(unsigned long int s)
Set the seed.
Definition: prob_dens_func.h:1914
o2scl::prob_dens_uniform::upper_limit
virtual double upper_limit() const
Uower limit of the range.
Definition: prob_dens_func.h:384
o2scl::prob_dens_hist::pdf
virtual double pdf(double x) const
The normalized density.
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::prob_dens_func::operator()
virtual double operator()() const
Sample from the specified density.
Definition: prob_dens_func.h:70
o2scl::prob_dens_mdim_bound_gaussian::set
void set(size_t p_ndim, vec_t &p_peak, mat_t &covar, vec_t &p_low, vec_t &p_high)
Set the peak, covariance matrix, lower limits, and upper limits.
Definition: prob_dens_func.h:1374
o2scl::prob_cond_mdim_indep::operator()
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
Definition: prob_dens_func.h:1826
o2scl::prob_dens_mdim_gaussian::get_peak
const vec_t & get_peak()
Get the peak location.
Definition: prob_dens_func.h:1036
o2scl::prob_dens_mdim_biv_gaussian::prob_dens_mdim_biv_gaussian
prob_dens_mdim_biv_gaussian(const prob_dens_mdim_biv_gaussian &pdmbg)
Copy constructor.
Definition: prob_dens_func.h:861
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::prob_dens_lognormal
Lognormal density function.
Definition: prob_dens_func.h:475
o2scl::prob_cond_mdim::log_pdf
virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const =0
The log of the conditional probability of x_A given x_B i.e. .
o2scl::prob_dens_mdim_biv_gaussian::y0
double y0
The y coordinate of the centroid.
Definition: prob_dens_func.h:844
o2scl::prob_cond_mdim_fixed_step::u_high
std::vector< double > u_high
Upper limits.
Definition: prob_dens_func.h:1600
o2scl::prob_dens_mdim_gaussian::operator()
virtual void operator()(vec_t &x) const
Sample the distribution.
Definition: prob_dens_func.h:1310
o2scl::prob_dens_mdim_gaussian::get_covar_inv
const mat_t & get_covar_inv()
Get the inverse of the covariance matrix.
Definition: prob_dens_func.h:1030
o2scl_const::pi
const double pi
Definition: constants.h:65
o2scl::prob_dens_mdim_gaussian::dim
virtual size_t dim() const
The dimensionality.
Definition: prob_dens_func.h:1047
o2scl::prob_cond_mdim_indep::dim
virtual size_t dim() const
The dimensionality.
Definition: prob_dens_func.h:1807
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::prob_dens_mdim_biv_gaussian
A bivariate gaussian probability distribution.
Definition: prob_dens_func.h:836
o2scl::prob_dens_mdim_gaussian::vtmp
vec_t vtmp
Temporary storage 2.
Definition: prob_dens_func.h:1010
o2scl::prob_dens_gaussian::prob_dens_gaussian
prob_dens_gaussian(const prob_dens_gaussian &pdg)
Copy constructor.
Definition: prob_dens_func.h:161
o2scl::prob_cond_mdim_indep
A multi-dimensional conditional probability density function independent of the input.
Definition: prob_dens_func.h:1778
o2scl::prob_dens_func::cdf
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
Definition: prob_dens_func.h:88
o2scl::prob_cond_mdim_fixed_step::log_pdf
virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const
The log of the conditional probability of x_A given x_B i.e. .
Definition: prob_dens_func.h:1742
o2scl::wvector_covariance
double wvector_covariance(size_t n, const vec_t &data1, const vec2_t &data2, const vec3_t &weights)
The weighted covariance of two vectors.
Definition: vec_stats.h:1354
o2scl_linalg::LU_decomp
int LU_decomp(const size_t N, mat_t &A, o2scl::permutation &p, int &signum)
Compute the LU decomposition of the matrix A.
Definition: lu_base.h:86
o2scl::prob_dens_lognormal::operator()
virtual double operator()() const
Sample from the specified density.
Definition: prob_dens_func.h:554
o2scl::prob_dens_mdim_gaussian::prob_dens_mdim_gaussian
prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar)
Create a distribution from the covariance matrix.
Definition: prob_dens_func.h:1135
o2scl::prob_dens_mdim_factor::pdf
virtual double pdf(const vec_t &x) const
The normalized density.
Definition: prob_dens_func.h:762
o2scl::prob_dens_gaussian::log_pdf
virtual double log_pdf(double x) const
The log of the normalized density.
Definition: prob_dens_func.h:234
o2scl::prob_cond_mdim
A multi-dimensional conditional probability density function.
Definition: prob_dens_func.h:1465
o2scl::prob_cond_mdim_indep::pdf
virtual double pdf(const vec_t &x_B, const vec_t &x_A) const
The conditional probability of x_A given x_B, i.e. .
Definition: prob_dens_func.h:1814
o2scl::prob_cond_mdim_fixed_step::operator()
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
Definition: prob_dens_func.h:1747
o2scl::prob_dens_lognormal::operator=
prob_dens_lognormal & operator=(const prob_dens_lognormal &pdg)
Copy constructor with operator=.
Definition: prob_dens_func.h:528
o2scl::prob_dens_mdim_factor::dim
virtual size_t dim() const
Return the dimensionality.
Definition: prob_dens_func.h:757
o2scl::prob_dens_mdim_gaussian::set_gproc
void set_gproc(size_t n_dim, size_t n_init, vec_vec_t &x, vec_t &y, func_t &fcovar)
Given a data set and a covariance function, construct probability distribution based on a Gaussian pr...
Definition: prob_dens_func.h:1217
o2scl::exc_eunimpl
@ exc_eunimpl
requested feature not (yet) implemented
Definition: err_hnd.h:99
o2scl::prob_dens_mdim_biv_gaussian::level_fixed_integral
virtual double level_fixed_integral(double integral)
Return the contour level corresponding to a fixed integral.
Definition: prob_dens_func.h:918
o2scl::prob_dens_mdim_factor::log_pdf
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
Definition: prob_dens_func.h:769
o2scl::prob_cond_mdim::pdf
virtual double pdf(const vec_t &x_B, const vec_t &x_A) const =0
The conditional probability of x_A given x_B, i.e. .
o2scl::prob_dens_mdim_gaussian::set
void set(size_t p_ndim, vec_t &p_peak, mat_t &covar)
Set the peak and covariance matrix for the distribution.
Definition: prob_dens_func.h:1144
o2scl::prob_dens_hist::cdf
virtual double cdf(double x) const
Cumulative distribution function (from the lower tail)
o2scl::prob_dens_mdim_gaussian::set
int set(size_t p_mdim, size_t n_pts, const mat2_t &pts, const vec2_t &vals, vec_t &peak_arg, mat_t &covar_arg)
Create a distribution from a set of samples from a multidimensional Gaussian, returning the peak valu...
Definition: prob_dens_func.h:1089
o2scl::prob_dens_mdim_bound_gaussian::log_pdf
virtual double log_pdf(const vec_t &x) const
Compute the natural log of the probability density function (arbitrary normalization)
Definition: prob_dens_func.h:1402
o2scl::prob_dens_gaussian::stddev
double stddev()
Get the Gaussian width.
Definition: prob_dens_func.h:207
o2scl::prob_dens_frange
A one-dimensional probability density over a finite range.
Definition: prob_dens_func.h:280
o2scl::prob_dens_lognormal::cdf
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
Definition: prob_dens_func.h:575
o2scl::prob_cond_mdim_gaussian::chol
mat_t chol
Cholesky decomposition.
Definition: prob_dens_func.h:1850
o2scl::prob_dens_hist::range
ubvector range
Vector specifying original histogram bins.
Definition: prob_dens_func.h:630
o2scl::prob_dens_gaussian::set_center
void set_center(double cent)
Set the center.
Definition: prob_dens_func.h:185
o2scl::prob_dens_lognormal::mu_
double mu_
Central value.
Definition: prob_dens_func.h:489
o2scl::prob_dens_lognormal::sigma_
double sigma_
Width parameter.
Definition: prob_dens_func.h:483
o2scl::prob_cond_mdim_gaussian::norm
double norm
Normalization factor.
Definition: prob_dens_func.h:1856
o2scl::prob_cond_mdim_gaussian::prob_cond_mdim_gaussian
prob_cond_mdim_gaussian()
Create an empty distribution.
Definition: prob_dens_func.h:1874
o2scl::prob_dens_hist::bin_ranges
const ubvector & bin_ranges()
Get reference to bin ranges.
Definition: prob_dens_func.h:648
o2scl::rng_gsl::random
double random()
Return a random number in .
Definition: rng_gsl.h:82
o2scl::rng_gsl::set_seed
void set_seed(unsigned long int s)
Set the seed.
Definition: rng_gsl.h:106
o2scl::prob_dens_uniform::r
rng_gsl r
The GSL random number generator.
Definition: prob_dens_func.h:312
o2scl::prob_dens_uniform::operator=
prob_dens_uniform & operator=(const prob_dens_uniform &pdg)
Copy constructor with operator=.
Definition: prob_dens_func.h:346
o2scl::prob_dens_frange::upper_limit
virtual double upper_limit() const =0
Uower limit of the range.
o2scl::prob_dens_mdim_gaussian::ndim
size_t ndim
Number of dimensions.
Definition: prob_dens_func.h:1004
o2scl::prob_cond_mdim_fixed_step::operator=
prob_cond_mdim_fixed_step & operator=(const prob_cond_mdim_fixed_step &pcmfs)
Copy constructor with operator=.
Definition: prob_dens_func.h:1650
o2scl::prob_dens_uniform::log_pdf
virtual double log_pdf(double x) const
The log of the normalized density.
Definition: prob_dens_func.h:412
o2scl::prob_dens_mdim_factor::list
std::vector< prob_dens_func > list
Vector of one-dimensional distributions.
Definition: prob_dens_func.h:733
o2scl::prob_cond_mdim_gaussian::q
vec_t q
Temporary storage 1.
Definition: prob_dens_func.h:1862
o2scl::prob_dens_mdim_gaussian::pdg
o2scl::prob_dens_gaussian pdg
Standard normal.
Definition: prob_dens_func.h:1020
o2scl::prob_cond_mdim_fixed_step::u_step
std::vector< double > u_step
Step sizes.
Definition: prob_dens_func.h:1592
o2scl::prob_dens_uniform::prob_dens_uniform
prob_dens_uniform(const prob_dens_uniform &pdg)
Copy constructor.
Definition: prob_dens_func.h:340
o2scl::prob_cond_mdim_fixed_step::set
virtual int set(vec_t &step, vec_t &low, vec_t &high)
Set step sizes and limits.
Definition: prob_dens_func.h:1688
o2scl::prob_cond_mdim_gaussian::operator()
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
Definition: prob_dens_func.h:1998
o2scl::prob_dens_hist::invert_cdf
virtual double invert_cdf(double x) const
Inverse cumulative distribution function (from the lower tail)
o2scl::prob_dens_mdim_bound_gaussian::operator()
virtual void operator()(vec_t &x) const
Sample the distribution.
Definition: prob_dens_func.h:1418
o2scl::prob_dens_mdim_factor::prob_dens_mdim_factor
prob_dens_mdim_factor(const prob_dens_mdim_factor &pdmf)
Copy constructor.
Definition: prob_dens_func.h:742
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::prob_dens_mdim_biv_gaussian::pdf
virtual double pdf(const vec_t &v) const
Compute the normalized probability density.
Definition: prob_dens_func.h:905
o2scl::prob_dens_mdim::pdf
virtual double pdf(const vec_t &x) const
The normalized density.
Definition: prob_dens_func.h:695
o2scl::prob_cond_mdim_gaussian::log_pdf
virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const
The log of the conditional probability of x_A given x_B i.e. .
Definition: prob_dens_func.h:1978
o2scl::prob_dens_hist::sum
ubvector sum
Normalized partial sum of histogram bins.
Definition: prob_dens_func.h:624
o2scl::prob_dens_uniform::operator()
virtual double operator()() const
Operator from the specified density.
Definition: prob_dens_func.h:393
o2scl::prob_cond_mdim_gaussian::set
void set(size_t p_ndim, mat_t &covar)
Set the covariance matrix for the distribution.
Definition: prob_dens_func.h:1921
o2scl::prob_dens_mdim_biv_gaussian::set
void set(double x_cent, double y_cent, double x_std, double y_std, double covar)
Set the properties of the distribution.
Definition: prob_dens_func.h:888
o2scl::prob_dens_hist::sv
search_vec< ubvector > sv
Search through the partial sums.
Definition: prob_dens_func.h:615
o2scl::prob_dens_uniform::lower_limit
virtual double lower_limit() const
Lower limit of the range.
Definition: prob_dens_func.h:375
o2scl::prob_cond_mdim_gaussian::prob_cond_mdim_gaussian
prob_cond_mdim_gaussian(const prob_cond_mdim_gaussian &pcmg)
Copy constructor.
Definition: prob_dens_func.h:1879
o2scl_linalg::cholesky_decomp
int cholesky_decomp(const size_t M, mat_t &A, bool err_on_fail=true)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix.
Definition: cholesky_base.h:62
o2scl::prob_dens_hist::operator()
virtual double operator()() const
Generate a sample.
o2scl::prob_dens_gaussian::operator=
prob_dens_gaussian & operator=(const prob_dens_gaussian &pdg)
Copy constructor with operator=.
Definition: prob_dens_func.h:168
o2scl::prob_dens_gaussian::sigma_
double sigma_
Width parameter.
Definition: prob_dens_func.h:129
o2scl::prob_dens_mdim::log_pdf
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
Definition: prob_dens_func.h:701
o2scl::prob_dens_mdim_gaussian::get_norm
const double & get_norm()
Get the normalization.
Definition: prob_dens_func.h:1042
o2scl::prob_dens_gaussian::set_seed
void set_seed(unsigned long int s)
Set the seed.
Definition: prob_dens_func.h:179
o2scl::prob_cond_mdim_indep::prob_cond_mdim_indep
prob_cond_mdim_indep(prob_dens_mdim< vec_t > &out)
Create a conditional probability distribution based on the specified probability distribution.
Definition: prob_dens_func.h:1789
o2scl::prob_dens_mdim_bound_gaussian::high
vec_t high
Upper limits.
Definition: prob_dens_func.h:1345
o2scl::prob_dens_lognormal::set_seed
void set_seed(unsigned long int s)
Set the seed.
Definition: prob_dens_func.h:549
o2scl::const_matrix_column_gen
Generic object which represents a column of a const matrix.
Definition: vector.h:3108
o2scl::prob_cond_mdim::operator()
virtual void operator()(const vec_t &x_B, vec_t &x_A) const =0
Sample the distribution.
o2scl::prob_dens_mdim_gaussian::chol
mat_t chol
Cholesky decomposition.
Definition: prob_dens_func.h:992
o2scl::prob_dens_mdim::operator()
virtual void operator()(vec_t &x) const
Sample the distribution.
Definition: prob_dens_func.h:717
o2scl::prob_dens_uniform::cdf
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
Definition: prob_dens_func.h:422
o2scl::prob_dens_gaussian::operator()
virtual double operator()() const
Sample from the specified density.
Definition: prob_dens_func.h:216
o2scl::prob_dens_uniform::prob_dens_uniform
prob_dens_uniform(double a, double b)
Create a uniform distribution from .
Definition: prob_dens_func.h:325
o2scl::prob_dens_mdim::dim
virtual size_t dim() const
Return the dimensionality.
Definition: prob_dens_func.h:689
o2scl::prob_dens_lognormal::invert_cdf
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
Definition: prob_dens_func.h:583
o2scl::search_vec
Searching class for monotonic data with caching.
Definition: search_vec.h:74
o2scl::prob_dens_uniform::prob_dens_uniform
prob_dens_uniform()
Create a blank uniform distribution.
Definition: prob_dens_func.h:318
o2scl::prob_dens_mdim_bound_gaussian
Gaussian distribution bounded by a hypercube.
Definition: prob_dens_func.h:1334
o2scl::prob_dens_gaussian::cdf
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
Definition: prob_dens_func.h:243
o2scl::prob_dens_mdim_gaussian::operator=
prob_dens_mdim_gaussian & operator=(const prob_dens_mdim_gaussian &pdmg_loc)
Copy constructor with operator=.
Definition: prob_dens_func.h:1067
o2scl::prob_cond_mdim_fixed_step::prob_cond_mdim_fixed_step
prob_cond_mdim_fixed_step(const prob_cond_mdim_fixed_step &pcmfs)
Copy constructor.
Definition: prob_dens_func.h:1643
o2scl::prob_dens_uniform::pdf
virtual double pdf(double x) const
The normalized density.
Definition: prob_dens_func.h:402
o2scl::prob_dens_mdim_gaussian::covar_inv
mat_t covar_inv
Inverse of the covariance matrix.
Definition: prob_dens_func.h:995
o2scl::prob_cond_mdim_fixed_step::set_seed
void set_seed(unsigned long int s)
Set the random number generator seed.
Definition: prob_dens_func.h:1663
o2scl::prob_dens_gaussian::set_sigma
void set_sigma(double sigma)
Set the Gaussian width (must be positive)
Definition: prob_dens_func.h:191
o2scl::prob_cond_mdim_gaussian::pdf
virtual double pdf(const vec_t &x_B, const vec_t &x_A) const
The conditional probability of x_A given x_B, i.e. .
Definition: prob_dens_func.h:1963
o2scl::prob_cond_mdim_gaussian
A multi-dimensional Gaussian conditional probability density function.
Definition: prob_dens_func.h:1845
o2scl::prob_dens_mdim_bound_gaussian::pdf
virtual double pdf(const vec_t &x) const
Compute the probability density function (arbitrary normalization)
Definition: prob_dens_func.h:1385
o2scl::prob_dens_mdim
A multi-dimensional probability density function.
Definition: prob_dens_func.h:684
o2scl::prob_dens_func::entropy
virtual double entropy() const
Entropy of the distribution ( )
Definition: prob_dens_func.h:100
o2scl::prob_dens_uniform::invert_cdf
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
Definition: prob_dens_func.h:433
o2scl::prob_dens_gaussian::r
o2scl::rng_gsl r
Base GSL random number generator.
Definition: prob_dens_func.h:132
o2scl::prob_dens_mdim_gaussian::prob_dens_mdim_gaussian
prob_dens_mdim_gaussian()
Create an empty distribution.
Definition: prob_dens_func.h:1052
o2scl::prob_dens_mdim_gaussian
A multi-dimensional Gaussian probability density function using a Cholesky decomposition.
Definition: prob_dens_func.h:987
o2scl::prob_dens_hist::init
void init(hist &h)
Initialize with histogram h.
o2scl::prob_dens_hist::n
size_t n
Number of original histogram bins.
Definition: prob_dens_func.h:618
o2scl::prob_dens_uniform::entropy
virtual double entropy() const
Entropy of the distribution ( )
Definition: prob_dens_func.h:446
o2scl::prob_cond_mdim_gaussian::ndim
size_t ndim
Number of dimensions.
Definition: prob_dens_func.h:1859
o2scl::prob_cond_mdim_fixed_step
Conditional probability for a random walk inside a hypercube.
Definition: prob_dens_func.h:1586
o2scl::prob_dens_mdim_biv_gaussian::contour
virtual void contour(double level, double theta, vec_t &x)
Return a point on the contour for a specified level given an angle.
Definition: prob_dens_func.h:932
o2scl::prob_dens_func::log_pdf
virtual double log_pdf(double x) const
The log of the normalized density.
Definition: prob_dens_func.h:82
o2scl::prob_dens_mdim_gaussian::q
vec_t q
Temporary storage 1.
Definition: prob_dens_func.h:1007
o2scl::prob_cond_mdim_gaussian::prob_cond_mdim_gaussian
prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar)
Create a distribution from the covariance matrix.
Definition: prob_dens_func.h:1904
o2scl::prob_dens_uniform::set_seed
void set_seed(unsigned long int s)
Set the seed.
Definition: prob_dens_func.h:355
o2scl::prob_dens_gaussian::mean
double mean()
Get the center.
Definition: prob_dens_func.h:202
o2scl::prob_dens_mdim_biv_gaussian::sig_y
double sig_y
The y standard deviation.
Definition: prob_dens_func.h:850
o2scl::prob_cond_mdim_fixed_step::u_low
std::vector< double > u_low
Lower limits.
Definition: prob_dens_func.h:1596
o2scl::prob_dens_uniform::ul
double ul
Upper limit.
Definition: prob_dens_func.h:309
o2scl::prob_dens_mdim_gaussian::set
int set(size_t p_mdim, size_t n_pts, const mat2_t &pts, const vec2_t &vals)
Create a distribution from a set of samples from a multidimensional Gaussian.
Definition: prob_dens_func.h:1123
o2scl::prob_dens_func::pdf
virtual double pdf(double x) const
The normalized density.
Definition: prob_dens_func.h:76
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::prob_dens_lognormal::set_mu_sigma
void set_mu_sigma(double mu, double sigma)
Set the maximum and width of the lognormal distribution.
Definition: prob_dens_func.h:538
o2scl::prob_dens_frange::lower_limit
virtual double lower_limit() const =0
Lower limit of the range.
o2scl::prob_dens_gaussian
A one-dimensional Gaussian probability density.
Definition: prob_dens_func.h:117
o2scl::prob_dens_func::invert_cdf
virtual double invert_cdf(double cdf) const
The inverse cumulative distribution function.
Definition: prob_dens_func.h:94
o2scl::prob_cond_mdim_indep::operator=
prob_cond_mdim_indep & operator=(const prob_cond_mdim_indep &pcmi)
Copy constructor with operator=.
Definition: prob_dens_func.h:1798
o2scl::prob_dens_lognormal::r
rng_gsl r
The GSL random number generator.
Definition: prob_dens_func.h:492
o2scl::prob_cond_mdim_gaussian::pdg
o2scl::prob_dens_gaussian pdg
Standard normal.
Definition: prob_dens_func.h:1868
o2scl::prob_dens_func
A one-dimensional probability density function.
Definition: prob_dens_func.h:65
o2scl::prob_dens_gaussian::cent_
double cent_
Central value.
Definition: prob_dens_func.h:123
o2scl::prob_dens_mdim_gaussian::pdf
virtual double pdf(const vec_t &x) const
The normalized density.
Definition: prob_dens_func.h:1282
o2scl::prob_dens_uniform
A uniform one-dimensional probability density over a finite range.
Definition: prob_dens_func.h:301
o2scl::prob_dens_mdim_bound_gaussian::prob_dens_mdim_bound_gaussian
prob_dens_mdim_bound_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar, vec_t &p_low, vec_t &p_high)
Create a distribution with the specified peak, covariance matrix, lower limits, and upper limits.
Definition: prob_dens_func.h:1362
o2scl::prob_dens_mdim_bound_gaussian::low
vec_t low
Lower limits.
Definition: prob_dens_func.h:1341
o2scl::prob_dens_lognormal::prob_dens_lognormal
prob_dens_lognormal(double mu, double sigma)
Create lognormal distribution with mean parameter mu and width parameter sigma.
Definition: prob_dens_func.h:508
o2scl::prob_dens_gaussian::prob_dens_gaussian
prob_dens_gaussian(double cent, double sigma)
Create a Gaussian distribution with width sigma.
Definition: prob_dens_func.h:147
o2scl::prob_dens_gaussian::invert_cdf
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
Definition: prob_dens_func.h:252
o2scl::prob_dens_lognormal::pdf
virtual double pdf(double x) const
The normalized density.
Definition: prob_dens_func.h:559
o2scl::prob_dens_mdim_biv_gaussian::sig_x
double sig_x
The x standard deviation.
Definition: prob_dens_func.h:847
o2scl::prob_dens_mdim_gaussian::peak
vec_t peak
Location of the peak.
Definition: prob_dens_func.h:998
o2scl::prob_cond_mdim_gaussian::operator=
prob_cond_mdim_gaussian & operator=(const prob_cond_mdim_gaussian &pcmg)
Copy constructor with operator=.
Definition: prob_dens_func.h:1889
o2scl::prob_cond_mdim_fixed_step::pdf
virtual double pdf(const vec_t &x_B, const vec_t &x_A) const
The conditional probability of x_A given x_B, i.e. .
Definition: prob_dens_func.h:1711
o2scl::prob_dens_mdim_gaussian::log_pdf
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
Definition: prob_dens_func.h:1297
o2scl_linalg::diagonal_has_zero
int diagonal_has_zero(const size_t N, mat_t &A)
Return 1 if at least one of the elements in the diagonal is zero.
Definition: lu_base.h:54
o2scl::prob_dens_hist::rng
rng_gsl rng
Random number generator.
Definition: prob_dens_func.h:633
o2scl::prob_cond_mdim_indep::prob_cond_mdim_indep
prob_cond_mdim_indep(const prob_cond_mdim_indep &pcmi)
Copy constructor.
Definition: prob_dens_func.h:1793
o2scl::prob_dens_hist::upper_limit
virtual double upper_limit() const
Uower limit of the range.
o2scl::prob_dens_mdim_biv_gaussian::x0
double x0
The x coordinate of the centroid.
Definition: prob_dens_func.h:841
o2scl::prob_dens_hist::entropy
virtual double entropy() const
Entropy of the distribution ( )
Definition: prob_dens_func.h:672
o2scl::prob_dens_mdim_gaussian::norm
double norm
Normalization factor, .
Definition: prob_dens_func.h:1001
o2scl::prob_dens_mdim_factor
A multidimensional distribution formed by the product of several one-dimensional distributions.
Definition: prob_dens_func.h:728
o2scl::prob_dens_uniform::set_limits
void set_limits(double a, double b)
Set the limits of the uniform distribution.
Definition: prob_dens_func.h:362
o2scl::prob_dens_lognormal::prob_dens_lognormal
prob_dens_lognormal(const prob_dens_lognormal &pdg)
Copy constructor.
Definition: prob_dens_func.h:522
o2scl::prob_dens_mdim_gaussian::set_alt
void set_alt(size_t p_ndim, vec_t &p_peak, mat_t &p_chol, mat_t &p_covar_inv, double p_norm)
Alternate set function for use when covariance matrix has already been decomposed and inverted.
Definition: prob_dens_func.h:1189
o2scl::prob_dens_mdim_biv_gaussian::rho
double rho
The covariance.
Definition: prob_dens_func.h:853

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).