prob_dens_func.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2018, 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  */
624  ubvector sum;
625 
626  /** \brief Vector specifying original histogram bins
627 
628  This vector has size \ref n plus one.
629  */
630  ubvector range;
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.ndim;
1059  chol=pdmg.chol;
1060  covar_inv=pdmg.covar_inv;
1061  norm=pdmg.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) {
1070  ndim=pdmg.ndim;
1071  chol=pdmg.chol;
1072  covar_inv=pdmg.covar_inv;
1073  norm=pdmg.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
1082  */
1083  template<class mat2_t, class vec2_t,
1084  class mat2_col_t=const_matrix_column_gen<mat2_t> >
1085  prob_dens_mdim_gaussian(size_t p_mdim, size_t n_pts, const mat2_t &pts,
1086  const vec2_t &vals) {
1087 
1088  vec_t peak2(p_mdim);
1089  mat_t covar(p_mdim,p_mdim);
1090 
1091  // Set peak with average and diagonal elements in covariance
1092  // matrix with variance
1093  for(size_t i=0;i<p_mdim;i++) {
1094  const mat2_col_t col(pts,i);
1095  peak2[i]=o2scl::wvector_mean<mat2_col_t>(n_pts,col,vals);
1096  // Square standard deviation
1097  covar(i,i)=o2scl::wvector_stddev<mat2_col_t>(n_pts,col,vals);
1098  covar(i,i)*=covar(i,i);
1099  }
1100  // Setup off-diagonal covariance matrix
1101  for(size_t i=0;i<p_mdim;i++) {
1102  mat2_col_t col_i(pts,i);
1103  for(size_t j=i+1;j<p_mdim;j++) {
1104  const mat2_col_t col_j(pts,j);
1105  double cov=o2scl::wvector_covariance(n_pts,col_i,col_j,vals);
1106  covar(i,j)=cov;
1107  covar(j,i)=cov;
1108  }
1109  }
1110  set(p_mdim,peak2,covar);
1111  }
1112 
1113  /** \brief Create a distribution from the covariance matrix
1114  */
1115  prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
1116  set(p_ndim,p_peak,covar);
1117  }
1118 
1119  /** \brief Set the peak and covariance matrix for the distribution
1120 
1121  \note This function is called in constructors and thus
1122  should not be virtual.
1123  */
1124  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
1125  if (p_ndim==0) {
1126  O2SCL_ERR("Zero dimension in prob_dens_mdim_gaussian::set().",
1128  }
1129  ndim=p_ndim;
1130  norm=1.0;
1131  peak.resize(ndim);
1132  for(size_t i=0;i<ndim;i++) peak[i]=p_peak[i];
1133  q.resize(ndim);
1134  vtmp.resize(ndim);
1135 
1136  // Perform the Cholesky decomposition of the covariance matrix
1137  chol=covar;
1138  o2scl_linalg::cholesky_decomp(ndim,chol);
1139 
1140  // Find the inverse
1141  covar_inv=chol;
1142  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1143 
1144  // Force chol to be lower triangular and compute the square root
1145  // of the determinant
1146  double sqrt_det=1.0;
1147  for(size_t i=0;i<ndim;i++) {
1148  if (!std::isfinite(chol(i,i))) {
1149  O2SCL_ERR2("An entry of the Cholesky decomposition was not finite ",
1150  "in prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1151  }
1152  sqrt_det*=chol(i,i);
1153  for(size_t j=0;j<ndim;j++) {
1154  if (i<j) chol(i,j)=0.0;
1155  }
1156  }
1157 
1158  // Compute normalization
1159  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt_det;
1160  if (!std::isfinite(norm)) {
1161  O2SCL_ERR2("Normalization not finite in ",
1162  "prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1163  }
1164  }
1165 
1166  /** \brief Alternate set function for use when covariance matrix
1167  has already been decomposed and inverted
1168  */
1169  void set_alt(size_t p_ndim, vec_t &p_peak, mat_t &p_chol,
1170  mat_t &p_covar_inv, double p_norm) {
1171  ndim=p_ndim;
1172  peak=p_peak;
1173  chol=p_chol;
1174  covar_inv=p_covar_inv;
1175  norm=p_norm;
1176  q.resize(ndim);
1177  vtmp.resize(ndim);
1178  return;
1179  }
1180 
1181  /** \brief Given a data set and a covariance function, construct
1182  probability distribution based on a Gaussian process which
1183  includes noise
1184 
1185  \note The type <tt>mat_col_t</tt> is a matrix column type
1186  for the internal object matrix type <tt>mat_t</tt>, and
1187  not associated with the data type <tt>vec_vec_t</tt>.
1188  Since the default matrix type is
1189  <tt>boost::numeric::ublas::matrix &lt; double &gt; </tt>
1190  a good matrix column type for this function is
1191  <tt>boost::numeric::ublas::matrix_column &lt;
1192  boost::numeric::ublas::matrix &lt; double &gt; &gt;</tt> .
1193  This matrix column type is needed for the LU
1194  decomposition and inversion.
1195  */
1196  template<class vec_vec_t, class mat_col_t, class func_t>
1197  void set_gproc(size_t n_dim, size_t n_init,
1198  vec_vec_t &x, vec_t &y, func_t &fcovar) {
1199 
1200  // Construct the four covariance matrices
1201 
1202  mat_t KXsX(n_dim,n_init);
1203  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1204  for(size_t icol=0;icol<n_init;icol++) {
1205  KXsX(irow-n_init,icol)=fcovar(x[irow],x[icol]);
1206  }
1207  }
1208 
1209  mat_t KXXs=boost::numeric::ublas::trans(KXsX);
1210 
1211  mat_t KXX(n_init,n_init);
1212  for(size_t irow=0;irow<n_init;irow++) {
1213  for(size_t icol=0;icol<n_init;icol++) {
1214  if (irow>icol) {
1215  KXX(irow,icol)=KXX(icol,irow);
1216  } else {
1217  KXX(irow,icol)=fcovar(x[irow],x[icol]);
1218  }
1219  }
1220  }
1221 
1222  mat_t KXsXs(n_dim,n_dim);
1223  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1224  for(size_t icol=n_init;icol<n_dim+n_init;icol++) {
1225  if (irow>icol) {
1226  KXsXs(irow-n_init,icol-n_init)=KXsXs(icol-n_init,irow-n_init);
1227  } else {
1228  KXsXs(irow-n_init,icol-n_init)=fcovar(x[irow],x[icol]);
1229  }
1230  }
1231  }
1232 
1233  // Construct the inverse of KXX
1234  mat_t inv_KXX(n_init,n_init);
1236  int signum;
1237  o2scl_linalg::LU_decomp(n_init,KXX,p,signum);
1238  if (o2scl_linalg::diagonal_has_zero(n_dim,KXX)) {
1239  O2SCL_ERR2("KXX matrix is singular in ",
1240  "prob_dens_mdim_gaussian::set_gproc().",
1242  }
1243  o2scl_linalg::LU_invert<mat_t,mat_t,mat_col_t>(n_init,KXX,p,inv_KXX);
1244 
1245  // Compute the mean vector
1246  vec_t prod(n_init), mean(n_dim);
1247  boost::numeric::ublas::axpy_prod(inv_KXX,y,prod,true);
1248  boost::numeric::ublas::axpy_prod(KXsX,prod,mean,true);
1249 
1250  // Compute the covariance matrix
1251  mat_t covar(n_dim,n_dim), prod2(n_init,n_dim), prod3(n_dim,n_dim);
1252  boost::numeric::ublas::axpy_prod(inv_KXX,KXXs,prod2,true);
1253  boost::numeric::ublas::axpy_prod(KXsX,prod2,prod3,true);
1254  covar=KXsXs-prod3;
1255 
1256  // Now use set() in the parent class
1257  this->set(n_dim,mean,covar);
1258 
1259  }
1260 
1261  /// The normalized density
1262  virtual double pdf(const vec_t &x) const {
1263  if (ndim==0) {
1264  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1265  "pdf().",o2scl::exc_einval);
1266  }
1267  double ret=norm;
1268  for(size_t i=0;i<ndim;i++) {
1269  q[i]=x[i]-peak[i];
1270  }
1271  vtmp=prod(covar_inv,q);
1272  ret*=exp(-0.5*inner_prod(q,vtmp));
1273  return ret;
1274  }
1275 
1276  /// The log of the normalized density
1277  virtual double log_pdf(const vec_t &x) const {
1278  if (ndim==0) {
1279  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1280  "pdf().",o2scl::exc_einval);
1281  }
1282  double ret=log(norm);
1283  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
1284  vtmp=prod(covar_inv,q);
1285  ret+=-0.5*inner_prod(q,vtmp);
1286  return ret;
1287  }
1288 
1289  /// Sample the distribution
1290  virtual void operator()(vec_t &x) const {
1291  if (ndim==0) {
1292  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1293  "operator().",o2scl::exc_einval);
1294  }
1295  for(size_t i=0;i<ndim;i++) q[i]=pdg();
1296  vtmp=prod(chol,q);
1297  for(size_t i=0;i<ndim;i++) x[i]=peak[i]+vtmp[i];
1298  return;
1299  }
1300 
1301  };
1302 
1303  /** \brief Gaussian distribution bounded by a hypercube
1304 
1305  \note This class naively resamples the Gaussian until
1306  a sample is within bounds. This is a temporary hack and
1307  can be very slow depending on the size of the volume
1308  excluded.
1309 
1310  \warning The PDF is not yet properly normalized
1311  */
1312  template<class vec_t=boost::numeric::ublas::vector<double>,
1313  class mat_t=boost::numeric::ublas::matrix<double> >
1315  public prob_dens_mdim_gaussian<vec_t,mat_t> {
1316 
1317  protected:
1318 
1319  /** \brief Lower limits
1320  */
1321  vec_t low;
1322 
1323  /** \brief Upper limits
1324  */
1325  vec_t high;
1326 
1327  public:
1328 
1329  /** \brief Maximum number of samples
1330  */
1331  size_t samp_max;
1332 
1333  /** \brief Create an empty distribution
1334  */
1336  samp_max=100000;
1337  }
1338 
1339  /** \brief Create a distribution with the specified peak, covariance
1340  matrix, lower limits, and upper limits
1341  */
1342  prob_dens_mdim_bound_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar,
1343  vec_t &p_low, vec_t &p_high) {
1344  set(p_ndim,p_peak,covar,p_low,p_high);
1345  samp_max=100000;
1346  }
1347 
1348  /** \brief Set the peak, covariance matrix, lower limits, and upper
1349  limits
1350 
1351  \note This function is called in constructors and thus
1352  should not be virtual.
1353  */
1354  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar,
1355  vec_t &p_low, vec_t &p_high) {
1356  prob_dens_mdim_gaussian<vec_t,mat_t>::set(p_ndim,p_peak,covar);
1357  low=p_low;
1358  high=p_high;
1359  return;
1360  }
1361 
1362  /** \brief Compute the probability density function (arbitrary
1363  normalization)
1364  */
1365  virtual double pdf(const vec_t &x) const {
1366  for(size_t i=0;i<this->ndim;i++) {
1367  if (x[i]<low[i]) {
1368  O2SCL_ERR("Parameter too small in pdf().",
1370  }
1371  if (x[i]>high[i]) {
1372  O2SCL_ERR("Parameter too large in pdf().",
1374  }
1375  }
1377  }
1378 
1379  /** \brief Compute the natural log of the probability density function
1380  (arbitrary normalization)
1381  */
1382  virtual double log_pdf(const vec_t &x) const {
1383  for(size_t i=0;i<this->ndim;i++) {
1384  if (x[i]<low[i]) {
1385  O2SCL_ERR("Parameter too small in log_pdf().",
1387  }
1388  if (x[i]>high[i]) {
1389  O2SCL_ERR("Parameter too large in log_pdf().",
1391  }
1392  }
1394  }
1395 
1396  /** \brief Sample the distribution
1397  */
1398  virtual void operator()(vec_t &x) const {
1399  bool done=false;
1400  size_t j=0;
1401  while (done==false) {
1402  done=true;
1404  j++;
1405  for(size_t i=0;i<this->ndim;i++) {
1406  if (x[i]<low[i]) {
1407  done=false;
1408  //std::cout << "Too small " << i << " " << x[i] << " "
1409  //<< low[i] << std::endl;
1410  i=this->ndim;
1411  } else if (x[i]>high[i]) {
1412  done=false;
1413  //std::cout << "Too large " << i << " " << x[i] << " "
1414  //<< high[i] << std::endl;
1415  i=this->ndim;
1416  }
1417  }
1418  if (j>samp_max) {
1419  O2SCL_ERR2("Sampling failed in ",
1420  "prob_dens_mdim_bound_gaussian::operator().",
1422  }
1423  }
1424  return;
1425  }
1426 
1427  };
1428 
1429  /** \brief A multi-dimensional conditional probability density function
1430 
1431  Note that conditional probabilities are typically written \f$
1432  P(A|B) \f$, i.e. the probability of \f$ A \f$ given \f$ B \f$.
1433  \o2 arranges the function parameters for the functions \ref
1434  o2scl::prob_cond_mdim::pdf, \ref o2scl::prob_cond_mdim::log_pdf
1435  \ref o2scl::prob_cond_mdim::operator()(), so that \f$ B \f$ is
1436  given first, and \f$ A \f$ is second.
1437 
1438  \ref
1439  o2scl::prob_cond_mdim::log_metrop_hast is a vector from \f$ B
1440  \f$ as denoted above.
1441 
1442  This class is experimental.
1443  */
1444  template<class vec_t=boost::numeric::ublas::vector<double> >
1446 
1447  public:
1448 
1449  /// The dimensionality
1450  virtual size_t dim() const {
1451  return 0;
1452  }
1453 
1454  /** \brief The conditional probability of x_A given x_B,
1455  i.e. \f$ P(A|B) \f$
1456  */
1457  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1458 
1459  /** \brief The log of the conditional probability of x_A given x_B
1460  i.e. \f$ \log [P(A|B)] \f$
1461  */
1462  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1463 
1464  /// Sample the distribution
1465  virtual void operator()(const vec_t &x_B, vec_t &x_A) const=0;
1466 
1467  /** \brief Sample the distribution and return the
1468  log of the Metropolis-Hastings ratio
1469 
1470  The Metropolis-Hastings ratio for a step beginning at \f$ x \f$
1471  and ending at \f$ x^{\prime} \f$ is
1472  obeys
1473  \f[
1474  \frac{P(x^{\prime})g(x|x^{\prime})}{P(x)g(x^{\prime}|x)}
1475  \f]
1476  taking the log, this gives
1477  \f[
1478  \log[P(x^{\prime})] - \log[P(x)] +
1479  \log \left[ \frac{g(x|x^{\prime})}{g(x^{\prime}|x)} \right]
1480  \f]
1481  thus this function computes
1482  \f[
1483  \log \left[ g(x|x^{\prime}) \right]
1484  - \log \left[ g(x^{\prime}|x) \right]
1485  \f]
1486  and thus, to keep a similar notation to
1487  \ref prob_cond_mdim::pdf() where \f$ g(x^{\prime}|x) \f$
1488  is obtained from
1489  \code
1490  pdf(x,x_prime)
1491  \endcode
1492  this function computes
1493  \code
1494  h(x,x_prime) = log_pdf(x_prime,x)-log_pdf(x,x_prime);
1495  \endcode
1496 
1497  To check this, in the limit that \f$ g(x|x^{\prime})
1498  \rightarrow P(x) \f$ this function returns
1499  \f[
1500  \log \left[ \frac{P(x)}{P(x^{\prime})} \right]
1501  \f]
1502 
1503  */
1504  virtual double log_metrop_hast(const vec_t &x, vec_t &x_prime) const {
1505  operator()(x,x_prime);
1506  double val1=log_pdf(x_prime,x);
1507  double val2=log_pdf(x,x_prime);
1508  if (!std::isfinite(val1)) {
1509  std::cout << "val1: " << val1 << std::endl;
1510  O2SCL_ERR("Log pdf not finite in log_metrop_hast 1.",
1512  }
1513  if (!std::isfinite(val2)) {
1514  std::cout << "val2: " << val2 << std::endl;
1515  O2SCL_ERR("Log pdf not finite in log_metrop_hast 2.",
1517  }
1518  return val1-val2;
1519  }
1520 
1521  };
1522 
1523  /** \brief Conditional probability for a random walk inside a
1524  hypercube
1525 
1526  \comment
1527  I had previously used std::uniform_real_distribution
1528  instead of rng_gsl, but this caused problems with
1529  intel compilers.
1530  \endcomment
1531 
1532  This conditional probability is most useful in providing a
1533  Metropolis-Hastings distribution with a fixed step size which
1534  properly handles a boundary. The Metropolis-Hastings step is
1535  accepted if \f$ r \in [0,1] \f$ obeys
1536  \f[
1537  r < \frac{P(x^{\prime})g(x|x^{\prime})}
1538  {P(x)g(x^{\prime}|x)}
1539  \f]
1540  The function \f$ g(x^{\prime}|x) = g(x^{\prime},x)/g(x) \f$, and
1541  because of the fixed step size these probabilities are just
1542  proportional to the inverse allowed volume, i.e. \f$ V(x)^{-1}
1543  V^{-1}(x^{\prime}) / V(x)^{-1} = V^{-1}(x^{\prime}) \f$ . If \f$
1544  x^{\prime} \f$ is near a boundary then \f$ V(x^{\prime}) \f$ is
1545  decreased and the conditional probability increases accordingly.
1546  If the distance between \f$ x \f$ and \f$ x^{\prime} \f$ is
1547  unreachable in a step, then the PDF is zero.
1548 
1549  \note This class is experimental.
1550 
1551  \note Const functions are not thread-safe because the
1552  class contains an internal mutable random number generator
1553  object.
1554 
1555  \comment
1556  If we do not include the g ratio, then the edges
1557  will be undersampled because we don't properly include
1558  the rejections
1559 
1560  For \f$ g(x^{\prime}|x) \f$, if x is near the edge, then the
1561  cond. prob. is larger, thus the g ratio is smaller than 1,
1562  encouraging more rejections.
1563  \endcomment
1564  */
1565  template<class vec_t=boost::numeric::ublas::vector<double> >
1567 
1568  protected:
1569 
1570  /** \brief Step sizes
1571  */
1572  std::vector<double> u_step;
1573 
1574  /** \brief Lower limits
1575  */
1576  std::vector<double> u_low;
1577 
1578  /** \brief Upper limits
1579  */
1580  std::vector<double> u_high;
1581 
1582  /** \brief Internal random number generator
1583  */
1584  mutable rng_gsl rg;
1585 
1586  /** \brief Internal set function
1587 
1588  \comment
1589  This can't be virtual because it needs to be called
1590  by the constructor
1591  \endcomment
1592  */
1593  int set_internal(size_t sz, vec_t &step, vec_t &low, vec_t &high) {
1594  u_step.resize(sz);
1595  u_low.resize(sz);
1596  u_high.resize(sz);
1597  for(size_t i=0;i<sz;i++) {
1598  u_step[i]=step[i];
1599 
1600  if (!std::isfinite(low[i]) || !std::isfinite(high[i])) {
1601  O2SCL_ERR2("Limit not finite in prob_cond_mdim_fixed_step::",
1602  "set_internal().",o2scl::exc_einval);
1603  }
1604 
1605  // Force low and high to be properly ordered
1606  if (low[i]>high[i]) {
1607  u_high[i]=low[i];
1608  u_low[i]=high[i];
1609  } else {
1610  u_low[i]=low[i];
1611  u_high[i]=high[i];
1612  }
1613  }
1614  return 0;
1615  }
1616 
1617  public:
1618 
1620  }
1621 
1622  /// Copy constructor
1624  u_step=pcmfs.u_step;
1625  u_low=pcmfs.u_low;
1626  u_high=pcmfs.u_high;
1627  }
1628 
1629  /// Copy constructor with operator=
1631  {
1632  // Check for self-assignment
1633  if (this!=&pcmfs) {
1634  u_step=pcmfs.u_step;
1635  u_low=pcmfs.u_low;
1636  u_high=pcmfs.u_high;
1637  }
1638  return *this;
1639  }
1640 
1641  /** \brief Set the random number generator seed
1642  */
1643  void set_seed(unsigned long int s) {
1644  rg.set_seed(s);
1645  return;
1646  }
1647 
1648  /** \brief Create a conditional probability object
1649  with specified step sizes and limits
1650  */
1651  template<class=vec_t> prob_cond_mdim_fixed_step
1652  (vec_t &step, vec_t &low, vec_t &high) {
1653  if (step.size()!=low.size()) {
1654  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1655  "prob_cond_mdim_fixed_step constructor.",
1657  }
1658  if (step.size()!=high.size()) {
1659  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1660  "prob_cond_mdim_fixed_step constructor.",
1662  }
1663  set_internal(step.size(),step,low,high);
1664  }
1665 
1666  /** \brief Set step sizes and limits
1667  */
1668  virtual int set(vec_t &step, vec_t &low, vec_t &high) {
1669  if (step.size()!=low.size()) {
1670  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1671  "prob_cond_mdim_fixed_step::set().",
1673  }
1674  if (step.size()!=high.size()) {
1675  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1676  "prob_cond_mdim_fixed_step::set().",
1678  }
1679  set_internal(step.size(),step,low,high);
1680  return 0;
1681  }
1682 
1683  /// The dimensionality
1684  virtual size_t dim() const {
1685  return u_step.size();
1686  }
1687 
1688  /** \brief The conditional probability of x_A given x_B,
1689  i.e. \f$ P(A|B) \f$
1690  */
1691  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1692  double vol1=1.0;
1693  for(size_t i=0;i<u_step.size();i++) {
1694 
1695  if (fabs(x_A[i]-x_B[i])>u_step[i]) return 0.0;
1696 
1697  if (x_B[i]-u_step[i]<u_low[i]) {
1698  if (x_B[i]+u_step[i]>u_high[i]) {
1699  // If x+step is too large and x-step is too small
1700  vol1*=u_high[i]-u_low[i];
1701  } else {
1702  // If x-step is too small
1703  vol1*=x_B[i]+u_step[i]-u_low[i];
1704  }
1705  } else {
1706  if (x_B[i]+u_step[i]>u_high[i]) {
1707  // If x+step is too large
1708  vol1*=u_high[i]-(x_B[i]-u_step[i]);
1709  } else {
1710  // The normal case, where the volumes are both inside
1711  // of the boundaries
1712  vol1*=2.0*u_step[i];
1713  }
1714  }
1715  }
1716  return 1.0/vol1;
1717  }
1718 
1719  /** \brief The log of the conditional probability of x_A given x_B
1720  i.e. \f$ \log [P(A|B)] \f$
1721  */
1722  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1723  return log(pdf(x_B,x_A));
1724  }
1725 
1726  /// Sample the distribution
1727  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1728  size_t nv=u_step.size();
1729  for(size_t i=0;i<nv;i++) {
1730  if (x_B[i]<u_low[i] || x_B[i]>u_high[i]) {
1731  std::cout << "Input out of bounds in fixed_step::operator(): "
1732  << i << " " << x_B[i] << " "
1733  << u_low[i] << " " << u_high[i] << std::endl;
1734  O2SCL_ERR("Input out of bounds in fixed_step::operator().",
1736  }
1737  do {
1738  x_A[i]=x_B[i]+u_step[i]*(rg.random()*2.0-1.0);
1739  } while (x_A[i]<u_low[i] || x_A[i]>u_high[i]);
1740  }
1741  return;
1742  }
1743 
1744  };
1745 
1746  /** \brief A multi-dimensional conditional probability density function
1747  independent of the input
1748 
1749  The conditional probability, \f$ P(A|B) = P(A,B)/P(B) \f$. If
1750  the joint probability is factorizable because the events \f$ A
1751  \f$ and \f$ B \f$ are independent, i.e. \f$ P(A,B) = P(A) P(B)
1752  \f$, then \f$ P(A|B) = P(A) \f$ and is independent of \f$ B \f$.
1753  This class handles that particular case.
1754 
1755  This class is experimental.
1756  */
1757  template<class vec_t=boost::numeric::ublas::vector<double> >
1758  class prob_cond_mdim_indep : public prob_cond_mdim<vec_t> {
1759 
1760  protected:
1761 
1762  prob_dens_mdim<vec_t> &base;
1763 
1764  public:
1765 
1766  /** \brief Create a conditional probability distribution
1767  based on the specified probability distribution
1768  */
1770  }
1771 
1772  /// Copy constructor
1774  base=pcmi.base;
1775  }
1776 
1777  /// Copy constructor with operator=
1779  // Check for self-assignment
1780  if (this!=&pcmi) {
1781  base=pcmi.base;
1782  }
1783  return *this;
1784  }
1785 
1786  /// The dimensionality
1787  virtual size_t dim() const {
1788  return base.dim();
1789  }
1790 
1791  /** \brief The conditional probability of x_A given x_B,
1792  i.e. \f$ P(A|B) \f$
1793  */
1794  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1795  return base.pdf(x_A);
1796  }
1797 
1798  /** \brief The log of the conditional probability of x_A given x_B
1799  i.e. \f$ \log [P(A|B)] \f$
1800  */
1801  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1802  return base.log_pdf(x_A);
1803  }
1804 
1805  /// Sample the distribution
1806  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1807  return base(x_A);
1808  }
1809 
1810  };
1811 
1812  /** \brief A multi-dimensional Gaussian conditional probability
1813  density function
1814 
1815  This class is experimental.
1816 
1817  \todo This should be a symmetric conditional probability,
1818  i.e. \f$ P(x|y) = P(y|x) \f$. Test this.
1819 
1820  \note Const functions are not thread-safe because
1821  mutable storage is used.
1822  */
1823  template<class vec_t=boost::numeric::ublas::vector<double>,
1824  class mat_t=boost::numeric::ublas::matrix<double> >
1825  class prob_cond_mdim_gaussian : public prob_cond_mdim<vec_t> {
1826 
1827  protected:
1828 
1829  /// Cholesky decomposition
1830  mat_t chol;
1831 
1832  /// Inverse of the covariance matrix
1833  mat_t covar_inv;
1834 
1835  /// Normalization factor
1836  double norm;
1837 
1838  /// Number of dimensions
1839  size_t ndim;
1840 
1841  /// Temporary storage 1
1842  mutable vec_t q;
1843 
1844  /// Temporary storage 2
1845  mutable vec_t vtmp;
1846 
1847  /// Standard normal
1849 
1850  public:
1851 
1852  /** \brief Create an empty distribution
1853  */
1855  ndim=0;
1856  }
1857 
1858  /// Copy constructor
1860  ndim=pcmg.ndim;
1861  chol=pcmg.chol;
1862  covar_inv=pcmg.covar_inv;
1863  norm=pcmg.norm;
1864  q.resize(ndim);
1865  vtmp.resize(ndim);
1866  }
1867 
1868  /// Copy constructor with operator=
1870  // Check for self-assignment
1871  if (this!=&pcmg) {
1872  ndim=pcmg.ndim;
1873  chol=pcmg.chol;
1874  covar_inv=pcmg.covar_inv;
1875  norm=pcmg.norm;
1876  q.resize(ndim);
1877  vtmp.resize(ndim);
1878  }
1879  return *this;
1880  }
1881 
1882  /** \brief Create a distribution from the covariance matrix
1883  */
1884  prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar) {
1885  set(p_ndim,covar);
1886  }
1887 
1888  /// The dimensionality
1889  virtual size_t dim() const {
1890  return ndim;
1891  }
1892 
1893  /// Set the seed
1894  void set_seed(unsigned long int s) {
1895  pdg.set_seed(s);
1896  return;
1897  }
1898 
1899  /** \brief Set the covariance matrix for the distribution
1900  */
1901  void set(size_t p_ndim, mat_t &covar) {
1902  if (p_ndim==0) {
1903  O2SCL_ERR("Zero dimension in prob_cond_mdim_gaussian::set().",
1905  }
1906  ndim=p_ndim;
1907  norm=1.0;
1908  q.resize(ndim);
1909  vtmp.resize(ndim);
1910 
1911  // Perform the Cholesky decomposition
1912  chol=covar;
1913  o2scl_linalg::cholesky_decomp(ndim,chol);
1914 
1915  // Find the inverse
1916  covar_inv=chol;
1917  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1918 
1919  // Force chol to be lower triangular and compute the determinant
1920  double sqrt_det=1.0;
1921  for(size_t i=0;i<ndim;i++) {
1922  if (!std::isfinite(chol(i,i))) {
1923  O2SCL_ERR2("An entry of the Cholesky decomposition was not finite ",
1924  "in prob_cond_mdim_gaussian::set().",o2scl::exc_einval);
1925  }
1926  sqrt_det*=chol(i,i);
1927  for(size_t j=0;j<ndim;j++) {
1928  if (i<j) chol(i,j)=0.0;
1929  }
1930  }
1931 
1932  // Compute normalization
1933  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt_det;
1934  if (!std::isfinite(norm)) {
1935  O2SCL_ERR2("Normalization not finite in ",
1936  "prob_dens_mdim_gaussian::set().",o2scl::exc_einval);
1937  }
1938  }
1939 
1940  /** \brief The conditional probability of x_A given x_B,
1941  i.e. \f$ P(A|B) \f$
1942  */
1943  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const {
1944  if (ndim==0) {
1945  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1946  "pdf().",o2scl::exc_einval);
1947  }
1948  double ret=norm;
1949  for(size_t i=0;i<ndim;i++) q[i]=x_A[i]-x_B[i];
1950  vtmp=prod(covar_inv,q);
1951  ret*=exp(-0.5*inner_prod(q,vtmp));
1952  return ret;
1953  }
1954 
1955  /** \brief The log of the conditional probability of x_A given x_B
1956  i.e. \f$ \log [P(A|B)] \f$
1957  */
1958  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const {
1959  if (ndim==0) {
1960  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1961  "pdf().",o2scl::exc_einval);
1962  }
1963  double ret=log(norm);
1964  for(size_t i=0;i<ndim;i++) q[i]=x_A[i]-x_B[i];
1965  vtmp=prod(covar_inv,q);
1966  ret+=-0.5*inner_prod(q,vtmp);
1967  /*
1968  std::cout << "pdmg lp: ";
1969  for (size_t i=0;i<ndim;i++) std::cout << x_A[i] << " ";
1970  for (size_t i=0;i<ndim;i++) std::cout << q[i] << " ";
1971  for (size_t i=0;i<ndim;i++) std::cout << vtmp[i] << " ";
1972  std::cout << ret << std::endl;
1973  */
1974  return ret;
1975  }
1976 
1977  /// Sample the distribution
1978  virtual void operator()(const vec_t &x_B, vec_t &x_A) const {
1979  if (ndim==0) {
1980  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1981  "operator().",o2scl::exc_einval);
1982  }
1983  for (size_t i=0;i<ndim;i++) q[i]=pdg();
1984  vtmp=prod(chol,q);
1985  for (size_t i=0;i<ndim;i++) x_A[i]=x_B[i]+vtmp[i];
1986  /*
1987  std::cout << "pdmg op: ";
1988  for (size_t i=0;i<ndim;i++) std::cout << x_A[i] << " ";
1989  std::cout << std::endl;
1990  */
1991  return;
1992  }
1993 
1994  };
1995 
1996 #ifdef O2SCL_NEVER_DEFINED
1997  /** \brief A multidimensional normal distribution from
1998  a Gaussian process
1999 
2000  \future The linear algebra only works with ublas and is
2001  not optimized.
2002  */
2003  template<class vec_t=boost::numeric::ublas::vector<double>,
2004  class mat_t=boost::numeric::ublas::matrix<double>,
2005  class mat_col_t=boost::numeric::ublas::matrix_column<mat_t> >
2006  class prob_dens_mdim_gproc :
2007  public o2scl::prob_dens_mdim_gaussian<vec_t> {
2008 
2009  public:
2010 
2011  };
2012 #endif
2013 
2014 #ifndef DOXYGEN_NO_O2NS
2015 }
2016 #endif
2017 
2018 #endif
vec_t peak
Location of the peak.
double sigma_
Width parameter.
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...
virtual double entropy() const
Entropy of the distribution ( )
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...
prob_cond_mdim_gaussian & operator=(const prob_cond_mdim_gaussian &pcmg)
Copy constructor with operator=.
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. .
rng_gsl rg
Internal random number generator.
Gaussian distribution bounded by a hypercube.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
prob_cond_mdim_indep & operator=(const prob_cond_mdim_indep &pcmi)
Copy constructor with operator=.
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
double mu_
Central value.
prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar)
Create a distribution from the covariance matrix.
void set_limits(double a, double b)
Set the limits of the uniform distribution.
double mean()
Get the center.
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
virtual size_t dim() const
The dimensionality.
int set_internal(size_t sz, vec_t &step, vec_t &low, vec_t &high)
Internal set function.
const double pi
Definition: constants.h:62
virtual double lower_limit() const
Lower limit of the range.
virtual double pdf(const vec_t &v) const
Compute the normalized probability density.
virtual double pdf(const vec_t &x) const
The normalized density.
virtual size_t dim() const
The dimensionality.
Conditional probability for a random walk inside a hypercube.
void set_seed(unsigned long int s)
Set the seed.
Definition: rng_gsl.h:106
prob_dens_gaussian & operator=(const prob_dens_gaussian &pdg)
Copy constructor with operator=.
prob_cond_mdim_gaussian(const prob_cond_mdim_gaussian &pcmg)
Copy constructor.
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. .
o2scl::prob_dens_gaussian pdg
Standard normal.
prob_dens_mdim_gaussian(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.
Generic object which represents a column of a const matrix.
Definition: vector.h:2647
void set_sigma(double sigma)
Set the Gaussian width (must be positive)
invalid argument supplied by user
Definition: err_hnd.h:59
prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar)
Create a distribution from the covariance matrix.
ubvector sum
Normalized partial sum of histogram bins.
A multi-dimensional conditional probability density function independent of the input.
size_t n
Number of original histogram bins.
const vec_t & get_peak()
Get the peak location.
vec_t vtmp
Temporary storage 2.
virtual double entropy() const
Entropy of the distribution ( )
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
const ubvector & bin_ranges()
Get reference to bin ranges.
virtual double level_fixed_integral(double integral)
Return the contour level corresponding to a fixed integral.
prob_dens_lognormal & operator=(const prob_dens_lognormal &pdg)
Copy constructor with operator=.
virtual void operator()(vec_t &x) const
Sample the distribution.
A class for representing permutations.
Definition: permutation.h:70
virtual double log_pdf(const vec_t &x) const
Compute the natural log of the probability density function (arbitrary normalization) ...
double sig_x
The x standard deviation.
prob_dens_lognormal()
Create a blank lognormal distribution.
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::rng_gsl r
Base GSL random number generator.
o2scl::prob_dens_gaussian pdg
Standard normal.
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. .
double random()
Return a random number in .
Definition: rng_gsl.h:82
double x0
The x coordinate of the centroid.
prob_cond_mdim_fixed_step & operator=(const prob_cond_mdim_fixed_step &pcmfs)
Copy constructor with operator=.
A one-dimensional Gaussian probability density.
prob_dens_uniform(double a, double b)
Create a uniform distribution from .
A multi-dimensional conditional probability density function.
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
rng_gsl r
The GSL random number generator.
prob_dens_gaussian()
Create a standard normal distribution.
mat_t chol
Cholesky decomposition.
virtual double operator()() const
Operator from the specified density.
virtual size_t dim() const
The dimensionality.
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
mat_t covar_inv
Inverse of the covariance matrix.
virtual double invert_cdf(double cdf) const
The inverse cumulative distribution function.
generic failure
Definition: err_hnd.h:61
void set_mu_sigma(double mu, double sigma)
Set the maximum and width of the lognormal distribution.
A one-dimensional histogram class.
Definition: hist.h:113
prob_dens_gaussian(const prob_dens_gaussian &pdg)
Copy constructor.
virtual double pdf(const vec_t &x) const
The normalized density.
const double & get_norm()
Get the normalization.
prob_dens_lognormal(const prob_dens_lognormal &pdg)
Copy constructor.
requested feature not (yet) implemented
Definition: err_hnd.h:99
std::vector< double > u_step
Step sizes.
vec_t vtmp
Temporary storage 2.
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
virtual size_t dim() const
The dimensionality.
prob_cond_mdim_gaussian()
Create an empty distribution.
A uniform one-dimensional probability density over a finite range.
prob_dens_gaussian(double cent, double sigma)
Create a Gaussian distribution with width sigma.
A multi-dimensional Gaussian conditional probability density function.
void set_seed(unsigned long int s)
Set the seed.
rng_gsl r
The GSL random number generator.
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
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:1520
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. .
A multidimensional distribution formed by the product of several one-dimensional distributions.
double norm
Normalization factor.
prob_dens_mdim_factor(const prob_dens_mdim_factor &pdmf)
Copy constructor.
prob_dens_mdim_gaussian(const prob_dens_mdim_gaussian &pdmg)
Copy constructor.
A one-dimensional probability density function.
virtual double pdf(double x) const
The normalized density.
virtual double entropy() const
Entropy of the distribution ( )
A one-dimensional probability density over a finite range.
virtual double entropy() const
Entropy of the distribution ( )
virtual double upper_limit() const
Uower limit of the range.
virtual double entropy() const
Entropy of the distribution ( )
virtual void operator()(vec_t &x) const
Sample the distribution.
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
const mat_t & get_chol()
Get the Cholesky decomposition.
A multi-dimensional probability density function.
const mat_t & get_covar_inv()
Get the inverse of the covariance matrix.
void set_seed(unsigned long int s)
Set the random number generator seed.
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
virtual void operator()(const vec_t &x_B, vec_t &x_A) const
Sample the distribution.
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
ubvector range
Vector specifying original histogram bins.
size_t samp_max
Maximum number of samples.
A bivariate gaussian probability distribution.
virtual double log_pdf(double x) const
The log of the normalized density.
double cent_
Central value.
virtual double pdf(const vec_t &x) const
The normalized density.
virtual size_t dim() const
Return the dimensionality.
search_vec< ubvector > sv
Search through the partial sums.
virtual double log_pdf(double x) const
The log of the normalized density.
virtual double operator()() const
Sample from the specified density.
std::vector< double > u_high
Upper limits.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
virtual double operator()() const
Sample from the specified density.
double sigma_
Width parameter.
vec_t q
Temporary storage 1.
virtual double log_pdf(double x) const
The log of the normalized density.
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. .
A one-dimensional probability density over the positive real numbers.
prob_dens_mdim_biv_gaussian(const prob_dens_mdim_biv_gaussian &pdmbg)
Copy constructor.
vec_t q
Temporary storage 1.
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.
double ul
Upper limit.
virtual double log_pdf(const vec_t &x) const
The log of the normalized density.
A multi-dimensional Gaussian probability density function using a Cholesky decomposition.
virtual size_t dim() const
The dimensionality.
prob_cond_mdim_indep(prob_dens_mdim< vec_t > &out)
Create a conditional probability distribution based on the specified probability distribution.
Random number generator (GSL)
Definition: rng_gsl.h:55
virtual void operator()(vec_t &x) const
Sample the distribution.
prob_dens_uniform()
Create a blank uniform distribution.
Searching class for monotonic data with caching.
Definition: search_vec.h:74
prob_dens_lognormal(double mu, double sigma)
Create lognormal distribution with mean parameter mu and width parameter sigma.
double norm
Normalization factor, .
virtual void contour(double level, double theta, vec_t &x)
Return a point on the contour for a specified level given an angle.
rng_gsl rng
Random number generator.
Probability density function based on a histogram.
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...
size_t ndim
Number of dimensions.
std::vector< prob_dens_func > list
Vector of one-dimensional distributions.
virtual double pdf(double x) const
The normalized density.
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. .
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
virtual size_t dim() const
Return the dimensionality.
virtual double operator()() const
Sample from the specified density.
mat_t covar_inv
Inverse of the covariance matrix.
virtual double pdf(const vec_t &x) const
Compute the probability density function (arbitrary normalization)
prob_dens_uniform(const prob_dens_uniform &pdg)
Copy constructor.
void set_seed(unsigned long int s)
Set the seed.
prob_dens_mdim_gaussian & operator=(const prob_dens_mdim_gaussian &pdmg)
Copy constructor with operator=.
prob_cond_mdim_fixed_step(const prob_cond_mdim_fixed_step &pcmfs)
Copy constructor.
prob_cond_mdim_indep(const prob_cond_mdim_indep &pcmi)
Copy constructor.
void set_seed(unsigned long int s)
Set the seed.
std::vector< double > u_low
Lower limits.
virtual void operator()(vec_t &x) const
Sample the distribution.
size_t ndim
Number of dimensions.
virtual double pdf(double x) const
The normalized density.
virtual double pdf(double x) const
The normalized density.
virtual double log_pdf(double x) const
The log of the normalized density.
prob_dens_mdim_gaussian()
Create an empty distribution.
void set_seed(unsigned long int s)
Set the seed.
const ubvector & partial_sums()
Get reference to partial sums.
double sig_y
The y standard deviation.
double ll
Lower limit.
void set(size_t p_ndim, vec_t &p_peak, mat_t &covar)
Set the peak and covariance matrix for the distribution.
Lognormal density function.
prob_dens_mdim_bound_gaussian()
Create an empty distribution.
double stddev()
Get the Gaussian width.
mat_t chol
Cholesky decomposition.
void set_center(double cent)
Set the center.
prob_dens_uniform & operator=(const prob_dens_uniform &pdg)
Copy constructor with operator=.
double y0
The y coordinate of the centroid.

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