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 
47 #ifndef DOXYGEN_NO_O2NS
48 namespace o2scl {
49 #endif
50 
51  /** \brief A one-dimensional probability density function
52 
53  This class is experimental.
54 
55  \future Give functions for mean, median, mode, variance, etc?
56 
57  \comment
58  For now, there aren't any pure virtual functions,
59  since this causes problems in creating an
60  std::vector<prob_dens_func> object below (especially
61  with intel compilers)
62  \endcomment
63  */
65 
66  public:
67 
68  /// Sample from the specified density
69  virtual double operator()() const {
70  return 0.0;
71  }
72 
73  /// The normalized density
74  virtual double pdf(double x) const {
75  return 0.0;
76  }
77 
78  /// The log of the normalized density
79  virtual double log_pdf(double x) const {
80  return 0.0;
81  }
82 
83  /// The cumulative distribution function (from the lower tail)
84  virtual double cdf(double x) const {
85  return 0.0;
86  }
87 
88  /// The inverse cumulative distribution function
89  virtual double invert_cdf(double cdf) const {
90  return 0.0;
91  }
92 
93  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
94  virtual double entropy() const {
95  return 0.0;
96  }
97 
98  };
99 
100  /** \brief A one-dimensional Gaussian probability density
101 
102  The distribution
103  \f[
104  P(x)=\frac{1}{\sigma \sqrt{2 \pi}}
105  e^{-\frac{\left(x-x_0\right)^2}{2\sigma^2}}
106  \f]
107 
108  This class is experimental.
109  */
111 
112  protected:
113 
114  /** \brief Central value
115  */
116  double cent_;
117 
118  /** \brief Width parameter
119 
120  A value of -1 indicates it is yet unspecified.
121  */
122  double sigma_;
123 
124  /// Base GSL random number generator
126 
127  public:
128 
129  /** \brief Create a standard normal distribution
130  */
132  cent_=0.0;
133  sigma_=1.0;
134  }
135 
136  /** \brief Create a Gaussian distribution with width \c sigma
137 
138  The value of \c sigma must be larger than zero.
139  */
140  prob_dens_gaussian(double cent, double sigma) {
141  if (sigma<0.0) {
142  O2SCL_ERR2("Tried to create a Gaussian dist. with sigma",
143  "<0 in prob_dens_gaussian::prob_dens_gaussian().",
144  exc_einval);
145  }
146  cent_=cent;
147  sigma_=sigma;
148  }
149 
150  virtual ~prob_dens_gaussian() {
151  }
152 
153  /// Copy constructor
155  cent_=pdg.cent_;
156  sigma_=pdg.sigma_;
157  r=pdg.r;
158  }
159 
160  /// Copy constructor with operator=
162  // Check for self-assignment
163  if (this!=&pdg) {
164  cent_=pdg.cent_;
165  sigma_=pdg.sigma_;
166  r=pdg.r;
167  }
168  return *this;
169  }
170 
171  /// Set the seed
172  void set_seed(unsigned long int s) {
173  r.set_seed(s);
174  return;
175  }
176 
177  /// Set the center
178  void set_center(double cent) {
179  cent_=cent;
180  return;
181  }
182 
183  /// Set the Gaussian width (must be positive)
184  void set_sigma(double sigma) {
185  if (sigma<0.0) {
186  O2SCL_ERR2("Tried to set negative sigma ",
187  "in prob_dens_gaussian::prob_dens_gaussian().",
188  exc_einval);
189  }
190  sigma_=sigma;
191  return;
192  }
193 
194  /// Get the center
195  double mean() {
196  return cent_;
197  }
198 
199  /// Get the Gaussian width
200  double stddev() {
201  if (sigma_<0.0) {
202  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
203  "get_sigma().",exc_einval);
204  }
205  return sigma_;
206  }
207 
208  /// Sample from the specified density
209  virtual double operator()() const {
210  if (sigma_<0.0) {
211  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
212  "operator().",exc_einval);
213  }
214  return cent_+gsl_ran_gaussian(&r,sigma_);
215  }
216 
217  /// The normalized density
218  virtual double pdf(double x) const {
219  if (sigma_<0.0) {
220  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
221  "pdf().",exc_einval);
222  }
223  return gsl_ran_gaussian_pdf(x-cent_,sigma_);
224  }
225 
226  /// The log of the normalized density
227  virtual double log_pdf(double x) const {
228  if (sigma_<0.0) {
229  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
230  "pdf().",exc_einval);
231  }
232  return log(gsl_ran_gaussian_pdf(x-cent_,sigma_));
233  }
234 
235  /// The cumulative distribution function (from the lower tail)
236  virtual double cdf(double x) const {
237  if (sigma_<0.0) {
238  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
239  "cdf().",exc_einval);
240  }
241  return gsl_cdf_gaussian_P(x-cent_,sigma_);
242  }
243 
244  /// The inverse cumulative distribution function
245  virtual double invert_cdf(double in_cdf) const {
246  if (sigma_<0.0) {
247  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
248  "invert_cdf().",exc_einval);
249  }
250  if (in_cdf<0.0 || in_cdf>1.0) {
251  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
252  "prob_dens_gaussian::invert_cdf().",exc_einval);
253  }
254  return gsl_cdf_gaussian_Pinv(in_cdf,sigma_)+cent_;
255  }
256 
257  /// The inverse cumulative distribution function
258  virtual double entropy() const {
259  if (sigma_<0.0) {
260  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
261  "invert_cdf().",exc_einval);
262  }
263  return log(2.0*o2scl_const::pi*exp(1.0)*sigma_*sigma_);
264  }
265 
266  };
267 
268  /** \brief A one-dimensional probability density over
269  a finite range
270 
271  This class is experimental.
272  */
274 
275  public:
276 
277  /// Lower limit of the range
278  virtual double lower_limit() const=0;
279 
280  /// Uower limit of the range
281  virtual double upper_limit() const=0;
282 
283  };
284 
285  /** \brief A uniform one-dimensional probability density
286  over a finite range
287 
288  A flat distribution given by \f$ P(x)=1/(b-a) \f$ for \f$ a<x<b
289  \f$, where \f$ a \f$ is the lower limit and \f$ b \f$ is the
290  upper limit.
291 
292  This class is experimental.
293  */
295 
296  protected:
297 
298  /// Lower limit
299  double ll;
300 
301  /// Upper limit
302  double ul;
303 
304  /// The GSL random number generator
306 
307  public:
308 
309  /** \brief Create a blank uniform distribution
310  */
312  ll=1.0;
313  ul=0.0;
314  }
315 
316  /** \brief Create a uniform distribution from \f$ a<x<b \f$
317  */
318  prob_dens_uniform(double a, double b) {
319  // Ensure a<b
320  if (a>b) {
321  double tmp=a;
322  a=b;
323  b=tmp;
324  }
325  ll=a;
326  ul=b;
327  }
328 
329  virtual ~prob_dens_uniform() {
330  }
331 
332  /// Copy constructor
334  ll=pdg.ll;
335  ul=pdg.ul;
336  }
337 
338  /// Copy constructor with operator=
340  // Check for self-assignment
341  if (this==&pdg) return *this;
342  ll=pdg.ll;
343  ul=pdg.ul;
344  return *this;
345  }
346 
347  /// Set the seed
348  void set_seed(unsigned long int s) {
349  r.set_seed(s);
350  return;
351  }
352 
353  /** \brief Set the limits of the uniform distribution
354  */
355  void set_limits(double a, double b) {
356  // Ensure a<b
357  if (a>b) {
358  double tmp=a;
359  a=b;
360  b=tmp;
361  }
362  ll=a;
363  ul=b;
364  return;
365  }
366 
367  /// Lower limit of the range
368  virtual double lower_limit() const {
369  if (ll>ul) {
370  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
371  "lower_limit().",exc_einval);
372  }
373  return ll;
374  }
375 
376  /// Uower limit of the range
377  virtual double upper_limit() const {
378  if (ll>ul) {
379  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
380  "upper_limit().",exc_einval);
381  }
382  return ul;
383  }
384 
385  /// Operator from the specified density
386  virtual double operator()() const {
387  if (ll>ul) {
388  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
389  "operator().",exc_einval);
390  }
391  return gsl_ran_flat(&r,ll,ul);
392  }
393 
394  /// The normalized density
395  virtual double pdf(double x) const {
396  if (ll>ul) {
397  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
398  "pdf().",exc_einval);
399  }
400  if (x<ll || x>ul) return 0.0;
401  return gsl_ran_flat_pdf(x,ll,ul);
402  }
403 
404  /// The log of the normalized density
405  virtual double log_pdf(double x) const {
406  if (ll>ul) {
407  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
408  "pdf().",exc_einval);
409  }
410  if (x<ll || x>ul) return 0.0;
411  return log(gsl_ran_flat_pdf(x,ll,ul));
412  }
413 
414  /// The cumulative distribution function (from the lower tail)
415  virtual double cdf(double x) const {
416  if (ll>ul) {
417  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
418  "cdf().",exc_einval);
419  }
420  if (x<ll) return 0.0;
421  if (x>ul) return 1.0;
422  return gsl_cdf_flat_P(x,ll,ul);
423  }
424 
425  /// The inverse cumulative distribution function
426  virtual double invert_cdf(double in_cdf) const {
427  if (ll>ul) {
428  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
429  "invert_cdf().",exc_einval);
430  }
431  if (in_cdf<0.0 || in_cdf>1.0) {
432  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
433  "prob_dens_uniform::invert_cdf().",exc_einval);
434  }
435  return gsl_cdf_flat_Pinv(in_cdf,ll,ul);
436  }
437 
438  /// The inverse cumulative distribution function
439  virtual double entropy() const {
440  if (ll>ul) {
441  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
442  "entropy().",exc_einval);
443  }
444  return log(ul-ll);
445  }
446 
447  };
448 
449  /** \brief A one-dimensional probability density over the
450  positive real numbers
451 
452  This class is experimental.
453  */
455 
456  };
457 
458  /** \brief Lognormal density function
459 
460  The distribution
461  \f[
462  P(x)=\frac{1}{x \sigma \sqrt{2 \pi}}
463  \exp \left[-\frac{\left(\ln x-\mu\right)^2}{2\sigma^2}\right]
464  \f]
465 
466  This class is experimental.
467  */
469 
470  protected:
471 
472  /** \brief Width parameter
473 
474  A value of -1 indicates it is yet unspecified.
475  */
476  double sigma_;
477 
478  /** \brief Central value
479 
480  A value of -1 indicates it is yet unspecified.
481  */
482  double mu_;
483 
484  /// The GSL random number generator
486 
487  public:
488 
489  /** \brief Create a blank lognormal distribution
490  */
492  sigma_=-1.0;
493  mu_=0.0;
494  }
495 
496  /** \brief Create lognormal distribution with mean parameter \c mu
497  and width parameter \c sigma
498 
499  The value of \c sigma must be larger than zero.
500  */
501  prob_dens_lognormal(double mu, double sigma) {
502  if (sigma<0.0) {
503  O2SCL_ERR2("Tried to create log normal dist. with mu or sigma",
504  "<0 in prob_dens_lognormal::prob_dens_lognormal().",
505  exc_einval);
506  }
507  mu_=mu;
508  sigma_=sigma;
509  }
510 
511  virtual ~prob_dens_lognormal() {
512  }
513 
514  /// Copy constructor
516  mu_=pdg.mu_;
517  sigma_=pdg.sigma_;
518  }
519 
520  /// Copy constructor with operator=
522  // Check for self-assignment
523  if (this==&pdg) return *this;
524  mu_=pdg.mu_;
525  sigma_=pdg.sigma_;
526  return *this;
527  }
528 
529  /** \brief Set the maximum and width of the lognormal distribution
530  */
531  void set_mu_sigma(double mu, double sigma) {
532  if (sigma<0.0) {
533  O2SCL_ERR2("Tried to set mu or sigma negative",
534  "in prob_dens_lognormal::prob_dens_lognormal().",
535  exc_einval);
536  }
537  mu_=mu;
538  sigma_=sigma;
539  }
540 
541  /// Set the seed
542  void set_seed(unsigned long int s) {
543  r.set_seed(s);
544  }
545 
546  /// Sample from the specified density
547  virtual double operator()() const {
548  return gsl_ran_lognormal(&r,mu_,sigma_);
549  }
550 
551  /// The normalized density
552  virtual double pdf(double x) const {
553  if (x<0.0) {
554  return 0.0;
555  }
556  return gsl_ran_lognormal_pdf(x,mu_,sigma_);
557  }
558 
559  /// The log of the normalized density
560  virtual double log_pdf(double x) const {
561  if (x<0.0) {
562  return 0.0;
563  }
564  return log(gsl_ran_lognormal_pdf(x,mu_,sigma_));
565  }
566 
567  /// The cumulative distribution function (from the lower tail)
568  virtual double cdf(double x) const {
569  if (x<0.0) {
570  return 0.0;
571  }
572  return gsl_cdf_lognormal_P(x,mu_,sigma_);
573  }
574 
575  /// The inverse cumulative distribution function
576  virtual double invert_cdf(double in_cdf) const {
577  if (in_cdf<0.0 || in_cdf>1.0) {
578  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
579  "prob_dens_lognormal::invert_cdf().",exc_einval);
580  }
581  return gsl_cdf_lognormal_Pinv(in_cdf,mu_,sigma_);
582  }
583 
584  /// The inverse cumulative distribution function
585  virtual double entropy() const {
586  if (sigma_<0.0) {
587  O2SCL_ERR2("Parameters not set in prob_dens_lognormal::",
588  "entropy().",exc_einval);
589  }
590  return 0.5+0.5*log(2.0*o2scl_const::pi*sigma_*sigma_)+mu_;
591  }
592 
593  };
594 
595  /** \brief Probability density function based on a histogram
596 
597  This class is experimental.
598  */
600 
601  public:
602 
604 
605  protected:
606 
607  /// Search through the partial sums
609 
610  /// Number of original histogram bins
611  size_t n;
612 
613  /** \brief Normalized partial sum of histogram bins
614 
615  This vector has size \ref n plus one.
616  */
617  ubvector sum;
618 
619  /** \brief Vector specifying original histogram bins
620 
621  This vector has size \ref n plus one.
622  */
623  ubvector range;
624 
625  /// Random number generator
626  mutable rng_gsl rng;
627 
628  public:
629 
630  prob_dens_hist();
631 
632  ~prob_dens_hist();
633 
634  /// Initialize with histogram \c h
635  void init(hist &h);
636 
637  /// Get reference to partial sums
638  const ubvector &partial_sums() { return sum; }
639 
640  /// Get reference to bin ranges
641  const ubvector &bin_ranges() { return range; }
642 
643  /// Generate a sample
644  virtual double operator()() const;
645 
646  /// Lower limit of the range
647  virtual double lower_limit() const;
648 
649  /// Uower limit of the range
650  virtual double upper_limit() const;
651 
652  /// The normalized density
653  virtual double pdf(double x) const;
654 
655  /// The log of the normalized density
656  virtual double log_pdf(double x) const;
657 
658  /// Cumulative distribution function (from the lower tail)
659  virtual double cdf(double x) const;
660 
661  /// Inverse cumulative distribution function (from the lower tail)
662  virtual double invert_cdf(double x) const;
663 
664  /// Inverse cumulative distribution function (from the lower tail)
665  virtual double entropy() const {
666  return 0.0;
667  }
668 
669  };
670 
671  /** \brief A multi-dimensional probability density function
672 
673  This class is experimental.
674  */
675  template<class vec_t=boost::numeric::ublas::vector<double> >
677 
678  public:
679 
680  /// Return the dimensionality
681  virtual size_t dim() const {
682  return 0;
683  }
684 
685  /// The normalized density
686  virtual double pdf(const vec_t &x) const {
687  return 0.0;
688  }
689 
690  /// The log of the normalized density
691  virtual double log_pdf(const vec_t &x) const {
692  return 0.0;
693  }
694 
695  /// Sample the distribution
696  virtual void operator()(vec_t &x) const {
697  return;
698  }
699 
700  };
701 
702  /** \brief A multidimensional distribution formed by the product
703  of several one-dimensional distributions
704  */
705  template<class vec_t=boost::numeric::ublas::vector<double> >
706  class prob_dens_mdim_factor : public prob_dens_mdim<vec_t> {
707 
708  protected:
709 
710  /// Vector of one-dimensional distributions
711  std::vector<prob_dens_func> list;
712 
713  public:
714 
715  prob_dens_mdim_factor(std::vector<prob_dens_func> &p_list) {
716  list=p_list;
717  }
718 
719  /// Return the dimensionality
720  virtual size_t dim() const {
721  return list.size();
722  }
723 
724  /// The normalized density
725  virtual double pdf(const vec_t &x) const {
726  double ret=1.0;
727  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
728  return ret;
729  }
730 
731  /// The log of the normalized density
732  virtual double log_pdf(const vec_t &x) const {
733  double ret=1.0;
734  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
735  return log(ret);
736  }
737 
738  /// Sample the distribution
739  virtual void operator()(vec_t &x) const {
740  for(size_t i=0;i<list.size();i++) x[i]=list[i]();
741  return;
742  }
743 
744  };
745 
746  /** \brief A bivariate gaussian probability distribution
747 
748  For a two-dimensional gaussian, given a mean
749  \f$ ( \mu_x, \mu_y ) \f$ and a covariance matrix
750  \f[
751  \Sigma = \left(
752  \begin{array}{cc}
753  \sigma_x^2 & \rho \sigma_x \sigma_y \\
754  \rho \sigma_x \sigma_y & \sigma_y^2 \\
755  \end{array}
756  \right)
757  \f]
758  the PDF is
759  \f[
760  pdf(x,y) = \left(2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}\right)^{-1}
761  \exp \left\{ - \frac{1}{2 (1-\rho^2)}
762  \left[ \frac{(x-\mu_x)^2}{\sigma_x^2} +
763  \frac{(y-\mu_y)^2}{\sigma_y^2} -
764  \frac{2 \rho (x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \right]
765  \right\}
766  \f]
767  (taken from the Wikipedia page on the "Multivariate normal
768  distribution").
769 
770  The function \ref o2scl::prob_dens_mdim_biv_gaussian::contour()
771  gives a point on the contour line for a fixed value of the
772  PDF given an angle \f$ \theta \f$. In particular,
773  it solves
774  \f[
775  c = pdf(r \cos \theta + \mu_x, r \sin \theta + \mu_y )
776  \f]
777  for the radius \f$ r \f$ and then stores the values
778  \f$ r \cos \theta + \mu_x \f$ and \f$ r \sin \theta + \mu_y \f$
779  in the reference parameters named \c x and \c y . Thus
780  this function can be used to map out the full contour
781  by selecting values for \f$ \theta \in [0,2 \pi] \f$.
782 
783  The function \ref
784  o2scl::prob_dens_mdim_biv_gaussian::level_fixed_integral() gives
785  the value of the PDF for which the integral inside the
786  corresponding contour is some fraction of the total integral
787  (which is always 1). Given a fraction \f$ f \f$, the argument of
788  the exponential is related to the inverse of the cumulative
789  distribution function for the chi-squared probability
790  distribution for two degrees of freedom for \f$ 1-f \f$. For a
791  fraction \f$ f \f$, the value \f$ \chi^2 \f$ (i.e. the
792  Mahalanobis distance) is \f$ \chi^2 = -2 \log (1-f) \f$ and then
793  the value of the PDF for the corresponding contour is \f$
794  pdf(x,y) = \left(2 \pi \sigma_x \sigma_y
795  \sqrt{1-\rho^2}\right)^{-1} \exp (-\chi^2/2) \f$ .
796 
797  */
798  template<class vec_t=boost::numeric::ublas::vector<double> >
800 
801  private:
802 
803  /// The x coordinate of the centroid
804  double x0;
805 
806  /// The y coordinate of the centroid
807  double y0;
808 
809  /// The x standard deviation
810  double sig_x;
811 
812  /// The y standard deviation
813  double sig_y;
814 
815  /// The covariance
816  double rho;
817 
818  public:
819 
821  }
822 
823  /** \brief Set the properties of the distribution
824 
825  \note If \f$ |\rho|\geq 1 \f$ this function will
826  call the error handler.
827  */
828  void set(double x_cent, double y_cent, double x_std, double y_std,
829  double covar) {
830  if (fabs(covar)>=1.0) {
831  O2SCL_ERR2("Covariance cannot have magnitude equal or larger than ",
832  "1 in prob_dens_mdim_biv_gaussian::set().",
834  }
835  x0=x_cent;
836  y0=y_cent;
837  sig_x=x_std;
838  sig_y=y_std;
839  rho=covar;
840  return;
841  }
842 
843  /** \brief Compute the normalized probability density
844  */
845  virtual double pdf(const vec_t &v) const {
846  double x=v[0], y=v[1];
847  double arg=-((x-x0)*(x-x0)/sig_x/sig_x+
848  (y-y0)*(y-y0)/sig_y/sig_y-
849  2.0*rho*(x-x0)*(y-y0)/sig_x/sig_y)/2.0/(1.0-rho*rho);
850  double ret=exp(arg)/2.0/o2scl_const::pi/sig_x/sig_y/
851  sqrt(1.0-rho*rho);
852  return ret;
853  }
854 
855  /** \brief Return the contour level corresponding to a fixed
856  integral
857  */
858  virtual double level_fixed_integral(double integral) {
859  // This comes from inverting the cumulative distribution function
860  // for the chi-squared distribution for two degrees of of freedom,
861  // i.e. exp(-x/2)
862  double arg=-2.0*log(1.0-integral);
863  // Now compute the pdf for the fixed value of the
864  // squared Mahalanobis distance
865  return exp(-0.5*arg)/2.0/o2scl_const::pi/sig_x/
866  sig_y/sqrt(1.0-rho*rho);
867  }
868 
869  /** \brief Return a point on the contour for a specified level
870  given an angle
871  */
872  virtual void contour(double level, double theta, vec_t &x) {
873  if (level<0.0) {
874  O2SCL_ERR2("Cannot produce contours for negative values in ",
875  "prob_dens_mdim_biv_gaussian::contour().",
877  }
878  double max=0.5/sig_x/sig_y/o2scl_const::pi/sqrt(1.0-rho*rho);
879  if (level>max) {
880  O2SCL_ERR2("Cannot produce contours larger than maximum in ",
881  "prob_dens_mdim_biv_gaussian::contour().",
883  }
884  double arg=-log(level*2.0*o2scl_const::pi*sig_x*sig_y*
885  sqrt(1.0-rho*rho))*2.0*(1.0-rho*rho);
886  double r2=arg/(cos(theta)*cos(theta)/sig_x/sig_x+
887  sin(theta)*sin(theta)/sig_y/sig_y-
888  2.0*rho/sig_x/sig_y*cos(theta)*sin(theta));
889  x[0]=sqrt(r2)*cos(theta)+x0;
890  x[1]=sqrt(r2)*sin(theta)+y0;
891  return;
892  }
893 
894  };
895 
896  /** \brief A multi-dimensional Gaussian probability density function
897  using a Cholesky decomposition
898 
899  Given a (square) covariance matrix, \f$ \Sigma \f$, and a mean
900  vector \f$ \mu \f$ the PDF is
901  \f[
902  P(x) = \det \left( 2 \pi \Sigma \right)^{-1/2}
903  \exp \left[ -\frac{1}{2} (x-\mu)^T \Sigma^{-1} (x-\mu) \right]
904  \f]
905 
906  Given the Cholesky decomposition \f$ A A^{T} = \Sigma \f$,
907  and a vector, \f$ z \f$ of samples from the standard Gaussian
908  with 0 mean and unit variance, one can create a sample
909  \f$ x \f$ from \f$ x = \mu + A z \f$ .
910 
911  \note This class inverts the matrix, necessary for computing the
912  pdf, but not for sampling the distribution, so for large
913  matrices the inversion can be a waste of computation if the pdf
914  is not needed.
915 
916  A separate class for the two-dimensional case is \ref
917  prob_dens_mdim_biv_gaussian .
918 
919  \future Create alternate versions based on other
920  matrix decompositions?
921  */
922  template<class vec_t=boost::numeric::ublas::vector<double>,
923  class mat_t=boost::numeric::ublas::matrix<double> >
924  class prob_dens_mdim_gaussian : public prob_dens_mdim<vec_t> {
925 
926  protected:
927 
928  /// Cholesky decomposition
929  mat_t chol;
930 
931  /// Inverse of the covariance matrix
932  mat_t covar_inv;
933 
934  /// Location of the peak
935  vec_t peak;
936 
937  /// Normalization factor, \f$ \det ( 2 \pi \Sigma)^{-1/2} \f$
938  double norm;
939 
940  /// Number of dimensions
941  size_t ndim;
942 
943  /// Temporary storage 1
944  mutable vec_t q;
945 
946  /// Temporary storage 2
947  mutable vec_t vtmp;
948 
949  /// Standard normal
951 
952  public:
953 
954  /// The dimensionality
955  virtual size_t dim() const {
956  return ndim;
957  }
958 
959  /// Create an empty distribution
961  ndim=0;
962  }
963 
964  /** \brief Create a distribution from the covariance matrix
965  */
966  prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
967  set(p_ndim,p_peak,covar);
968  }
969 
970  /** \brief Set the peak and covariance matrix for the distribution
971  */
972  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
973  if (p_ndim==0) {
974  O2SCL_ERR("Zero dimension in prob_dens_mdim_gaussian::set().",
976  }
977  ndim=p_ndim;
978  norm=1.0;
979  peak.resize(ndim);
980  for(size_t i=0;i<ndim;i++) peak[i]=p_peak[i];
981  q.resize(ndim);
982  vtmp.resize(ndim);
983 
984  // Perform the Cholesky decomposition of the covariance matrix
985  chol=covar;
987 
988  // Find the inverse
989  covar_inv=chol;
990  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
991 
992  // Force chol to be lower triangular and compute the determinant
993  double det=1.0;
994  for(size_t i=0;i<ndim;i++) {
995  det*=chol(i,i);
996  for(size_t j=0;j<ndim;j++) {
997  if (i<j) chol(i,j)=0.0;
998  }
999  }
1000  det*=det;
1001 
1002  // Compute normalization
1003  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt(det);
1004  }
1005 
1006  /** \brief Alternate set function for use when covariance matrix
1007  has already been decomposed and inverted
1008  */
1009  void set_alt(size_t p_ndim, vec_t &p_peak, mat_t &p_chol,
1010  mat_t &p_covar_inv, double p_norm) {
1011  ndim=p_ndim;
1012  peak=p_peak;
1013  chol=p_chol;
1014  covar_inv=p_covar_inv;
1015  norm=p_norm;
1016  return;
1017  }
1018 
1019  /** \brief Given a data set and a covariance function, construct
1020  probability distribution based on a Gaussian process which
1021  includes noise
1022  */
1023  template<class vec_vec_t, class mat_col_t, class func_t>
1024  void set_gproc(size_t n_dim, size_t n_init,
1025  vec_vec_t &x, vec_t &y, func_t &fcovar) {
1026 
1027  // Construct the four covariance matrices
1028 
1029  mat_t KXsX(n_dim,n_init);
1030  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1031  for(size_t icol=0;icol<n_init;icol++) {
1032  KXsX(irow-n_init,icol)=fcovar(x[irow],x[icol]);
1033  }
1034  }
1035 
1036  mat_t KXXs=boost::numeric::ublas::trans(KXsX);
1037 
1038  mat_t KXX(n_init,n_init);
1039  for(size_t irow=0;irow<n_init;irow++) {
1040  for(size_t icol=0;icol<n_init;icol++) {
1041  if (irow>icol) {
1042  KXX(irow,icol)=KXX(icol,irow);
1043  } else {
1044  KXX(irow,icol)=fcovar(x[irow],x[icol]);
1045  }
1046  }
1047  }
1048 
1049  mat_t KXsXs(n_dim,n_dim);
1050  for(size_t irow=n_init;irow<n_dim+n_init;irow++) {
1051  for(size_t icol=n_init;icol<n_dim+n_init;icol++) {
1052  if (irow>icol) {
1053  KXsXs(irow-n_init,icol-n_init)=KXsXs(icol-n_init,irow-n_init);
1054  } else {
1055  KXsXs(irow-n_init,icol-n_init)=fcovar(x[irow],x[icol]);
1056  }
1057  }
1058  }
1059 
1060  // Construct the inverse of KXX
1061  mat_t inv_KXX(n_init,n_init);
1063  int signum;
1064  o2scl_linalg::LU_decomp(n_init,KXX,p,signum);
1065  if (o2scl_linalg::diagonal_has_zero(n_dim,KXX)) {
1066  O2SCL_ERR2("KXX matrix is singular in ",
1067  "prob_dens_mdim_gaussian::set_gproc().",
1069  }
1070  o2scl_linalg::LU_invert<mat_t,mat_t,mat_col_t>(n_init,KXX,p,inv_KXX);
1071 
1072  // Compute the mean vector
1073  vec_t prod(n_init), mean(n_dim);
1074  boost::numeric::ublas::axpy_prod(inv_KXX,y,prod,true);
1075  boost::numeric::ublas::axpy_prod(KXsX,prod,mean,true);
1076 
1077  // Compute the covariance matrix
1078  mat_t covar(n_dim,n_dim), prod2(n_init,n_dim), prod3(n_dim,n_dim);
1079  boost::numeric::ublas::axpy_prod(inv_KXX,KXXs,prod2,true);
1080  boost::numeric::ublas::axpy_prod(KXsX,prod2,prod3,true);
1081  covar=KXsXs-prod3;
1082 
1083  // Now use set() in the parent class
1084  this->set(n_dim,mean,covar);
1085 
1086  }
1087 
1088  /// The normalized density
1089  virtual double pdf(const vec_t &x) const {
1090  if (ndim==0) {
1091  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1092  "pdf().",o2scl::exc_einval);
1093  }
1094  double ret=norm;
1095  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
1096  vtmp=prod(covar_inv,q);
1097  ret*=exp(-0.5*inner_prod(q,vtmp));
1098  return ret;
1099  }
1100 
1101  /// The log of the normalized density
1102  virtual double log_pdf(const vec_t &x) const {
1103  if (ndim==0) {
1104  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1105  "pdf().",o2scl::exc_einval);
1106  }
1107  double ret=log(norm);
1108  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
1109  vtmp=prod(covar_inv,q);
1110  ret+=-0.5*inner_prod(q,vtmp);
1111  return ret;
1112  }
1113 
1114  /// Sample the distribution
1115  virtual void operator()(vec_t &x) const {
1116  if (ndim==0) {
1117  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
1118  "operator().",o2scl::exc_einval);
1119  }
1120  for(size_t i=0;i<ndim;i++) q[i]=pdg();
1121  vtmp=prod(chol,q);
1122  for(size_t i=0;i<ndim;i++) x[i]=peak[i]+vtmp[i];
1123  return;
1124  }
1125 
1126  };
1127 
1128  /** \brief A multi-dimensional conditional probability density function
1129 
1130  Note that conditional probabilities are typically written \f$
1131  P(A|B) \f$, i.e. the probability of \f$ A \f$ given \f$ B \f$ ,
1132  so when sampling the conditional probability density you specify
1133  a vector in \f$ B \f$ (the "input") and get a randomly chosen
1134  vector in \f$ A \f$ (the "output") . \o2 arranges
1135  function parameters so that input parameters are first and
1136  output parameters are last, thus
1137  the first argument to the functions \ref
1138  o2scl::prob_cond_mdim::pdf, \ref o2scl::prob_cond_mdim::log_pdf
1139  \ref o2scl::prob_cond_mdim::operator()(), and \ref
1140  o2scl::prob_cond_mdim::log_metrop_hast is a vector from \f$ B
1141  \f$ as denoted above.
1142 
1143  This class is experimental.
1144  */
1145  template<class vec_t=boost::numeric::ublas::vector<double> >
1147 
1148  public:
1149 
1150  /// The dimensionality
1151  virtual size_t dim() const {
1152  return 0;
1153  }
1154 
1155  /// The conditional probability
1156  virtual double pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1157 
1158  /// The log of the conditional probability
1159  virtual double log_pdf(const vec_t &x_B, const vec_t &x_A) const=0;
1160 
1161  /// Sample the distribution
1162  virtual void operator()(const vec_t &x_B, vec_t &x_A) const=0;
1163 
1164  /** \brief Sample the distribution and return the
1165  log of the Metropolis-Hastings ratio
1166 
1167  The Metropolis-Hastings ratio for a step beginning at \f$ x \f$
1168  and ending at \f$ x^{\prime} \f$ is
1169  obeys
1170  \f[
1171  \frac{P(x^{\prime})g(x|x^{\prime})}{P(x)g(x^{\prime}|x)}
1172  \f]
1173  thus this function computes
1174  \f[
1175  \log \left[ \frac{g(x|x^{\prime})}{g(x^{\prime}|x)} \right]
1176  \f]
1177  and thus this function takes <tt>x_B</tt> as input, obtains
1178  a sample <tt>x_A</tt> and returns the value
1179  \code
1180  log_pdf(x_A,x_B)-log_pdf(x_B,x_A);
1181  \endcode
1182 
1183  To check this, imagine that \f$ g \f$ is independent
1184  of the "input" argument, so we have
1185  \f[
1186  \log \left[ \frac{g(x)}{g(x^{\prime})} \right]
1187  \f]
1188  and then a Metropolis-Hastings step is more likely accepted
1189  if \f$ g(x^{\prime}) \f$ is an underestimate.
1190  */
1191  virtual double log_metrop_hast(const vec_t &x_B, vec_t &x_A) const {
1192  operator()(x_B,x_A);
1193  return log_pdf(x_A,x_B)-log_pdf(x_B,x_A);
1194  }
1195 
1196  };
1197 
1198  /** \brief Conditional probability for a random walk inside a
1199  hypercube
1200 
1201  \comment
1202  I had previously used std::uniform_real_distribution
1203  instead of rng_gsl, but this caused problems with
1204  intel compilers.
1205  \endcomment
1206 
1207  This conditional probability is most useful in providing a
1208  Metropolis-Hastings distribution with a fixed step size which
1209  properly handles a boundary. The Metropolis-Hastings step is
1210  accepted if \f$ r \in [0,1] \f$ obeys
1211  \f[
1212  r < \frac{P(x^{\prime})g(x|x^{\prime})}
1213  {P(x)g(x^{\prime}|x)}
1214  \f]
1215  The function \f$ g(x^{\prime}|x) = g(x^{\prime},x)/g(x) \f$, and
1216  because of the fixed step size these probabilities are just
1217  proportional to the inverse allowed volume, i.e. \f$ V(x)^{-1}
1218  V^{-1}(x^{\prime}) / V(x)^{-1} = V^{-1}(x^{\prime}) \f$ . If \f$
1219  x^{\prime} \f$ is near a boundary then \f$ V(x^{\prime}) \f$ is
1220  decreased and the conditional probability increases accordingly.
1221  If the distance between \f$ x \f$ and \f$ x^{\prime} \f$ is
1222  unreachable in a step, then the PDF is zero.
1223 
1224  \note This class is experimental.
1225 
1226  \comment
1227  If we do not include the g ratio, then the edges
1228  will be undersampled because we don't properly include
1229  the rejections
1230 
1231  For \f$ g(x^{\prime}|x) \f$, if x is near the edge, then the
1232  cond. prob. is larger, thus the g ratio is smaller than 1,
1233  encouraging more rejections.
1234  \endcomment
1235  */
1236  template<class vec_t=boost::numeric::ublas::vector<double> >
1238 
1239  protected:
1240 
1241  /** \brief Step sizes
1242  */
1243  std::vector<double> u_step;
1244 
1245  /** \brief Lower limits
1246  */
1247  std::vector<double> u_low;
1248 
1249  /** \brief Upper limits
1250  */
1251  std::vector<double> u_high;
1252 
1253  /** \brief Internal random number generator
1254  */
1255  mutable rng_gsl rg;
1256 
1257  /** \brief Internal set function
1258 
1259  \comment
1260  This can't be virtual because it needs to be called
1261  by the constructor
1262  \endcomment
1263  */
1264  int set_internal(size_t sz, vec_t &step, vec_t &low, vec_t &high) {
1265  u_step.resize(sz);
1266  u_low.resize(sz);
1267  u_high.resize(sz);
1268  for(size_t i=0;i<sz;i++) {
1269  u_step[i]=step[i];
1270 
1271  if (!std::isfinite(low[i]) || !std::isfinite(high[i])) {
1272  O2SCL_ERR2("Limit not finite in prob_cond_mdim_fixed_step::",
1273  "set_internal().",o2scl::exc_einval);
1274  }
1275 
1276  // Force low and high to be properly ordered
1277  if (low[i]>high[i]) {
1278  u_high[i]=low[i];
1279  u_low[i]=high[i];
1280  } else {
1281  u_low[i]=low[i];
1282  u_high[i]=high[i];
1283  }
1284  }
1285  return 0;
1286  }
1287 
1288  public:
1289 
1291  }
1292 
1293  /** \brief Set the random number generator seed
1294  */
1295  void set_seed(unsigned long int s) {
1296  rg.set_seed(s);
1297  return;
1298  }
1299 
1300  /** \brief Create a conditional probability object
1301  with specified step sizes and limits
1302  */
1303  template<class=vec_t> prob_cond_mdim_fixed_step
1304  (vec_t &step, vec_t &low, vec_t &high) {
1305  if (step.size()!=low.size()) {
1306  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1307  "prob_cond_mdim_fixed_step constructor.",
1309  }
1310  if (step.size()!=high.size()) {
1311  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1312  "prob_cond_mdim_fixed_step constructor.",
1314  }
1315  set_internal(step.size(),step,low,high);
1316  }
1317 
1318  /** \brief Set step sizes and limits
1319  */
1320  virtual int set(vec_t &step, vec_t &low, vec_t &high) {
1321  if (step.size()!=low.size()) {
1322  O2SCL_ERR2("Vectors 'step' and 'low' mismatched in ",
1323  "prob_cond_mdim_fixed_step::set().",
1325  }
1326  if (step.size()!=high.size()) {
1327  O2SCL_ERR2("Vectors 'step' and 'high' mismatched in ",
1328  "prob_cond_mdim_fixed_step::set().",
1330  }
1331  set_internal(step.size(),step,low,high);
1332  return 0;
1333  }
1334 
1335  /// The dimensionality
1336  virtual size_t dim() const {
1337  return u_step.size();
1338  }
1339 
1340  /// The conditional probability
1341  virtual double pdf(const vec_t &x, const vec_t &x2) const {
1342  double vol1=1.0;
1343  for(size_t i=0;i<u_step.size();i++) {
1344 
1345  if (fabs(x2[i]-x[i])>u_step[i]) return 0.0;
1346 
1347  if (x[i]-u_step[i]<u_low[i]) {
1348  if (x[i]+u_step[i]>u_high[i]) {
1349  // If x+step is too large and x-step is too small
1350  vol1*=u_high[i]-u_low[i];
1351  } else {
1352  // If x-step is too small
1353  vol1*=x[i]+u_step[i]-u_low[i];
1354  }
1355  } else {
1356  if (x[i]+u_step[i]>u_high[i]) {
1357  // If x+step is too large
1358  vol1*=u_high[i]-(x[i]-u_step[i]);
1359  } else {
1360  // The normal case, where the volumes are both inside
1361  // of the boundaries
1362  vol1*=2.0*u_step[i];
1363  }
1364  }
1365  }
1366  return 1.0/vol1;
1367  }
1368 
1369  /// The log of the conditional probability
1370  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
1371  return log(pdf(x,x2));
1372  }
1373 
1374  /// Sample the distribution
1375  virtual void operator()(const vec_t &x, vec_t &x2) const {
1376  size_t nv=u_step.size();
1377  for(size_t i=0;i<nv;i++) {
1378  if (x[i]<u_low[i] || x[i]>u_high[i]) {
1379  std::cout << "Herex: " << i << " " << x[i] << " "
1380  << u_low[i] << " " << u_high[i] << std::endl;
1381  O2SCL_ERR("Input out of bounds in fixed_step::operator().",
1383  }
1384  do {
1385  x2[i]=x[i]+u_step[i]*(rg.random()*2.0-1.0);
1386  } while (x2[i]<u_low[i] || x2[i]>u_high[i]);
1387  }
1388  return;
1389  }
1390 
1391  };
1392 
1393  /** \brief A multi-dimensional conditional probability density function
1394  independent of the input
1395 
1396  The conditional probability, \f$ P(A|B) = P(A,B)/P(B) \f$. If
1397  the joint probability is factorizable because the events \f$ A
1398  \f$ and \f$ B \f$ are independent, i.e. \f$ P(A,B) = P(A) P(B)
1399  \f$, then \f$ P(A|B) = P(A) \f$ and is independent of \f$ B \f$.
1400  This class handles that particular case.
1401 
1402  This class is experimental.
1403  */
1404  template<class vec_t=boost::numeric::ublas::vector<double> >
1405  class prob_cond_mdim_indep : public prob_cond_mdim<vec_t> {
1406 
1407  protected:
1408 
1409  prob_dens_mdim<vec_t> &base;
1410 
1411  public:
1412 
1413  /** \brief Create a conditional probability distribution
1414  based on the specified probability distribution
1415  */
1417  }
1418 
1419  /// The dimensionality
1420  virtual size_t dim() const {
1421  return base.dim();
1422  }
1423 
1424  /// The conditional probability
1425  virtual double pdf(const vec_t &x, const vec_t &x2) const {
1426  return base.pdf(x2);
1427  }
1428 
1429  /// The log of the conditional probability
1430  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
1431  return base.log_pdf(x2);
1432  }
1433 
1434  /// Sample the distribution
1435  virtual void operator()(const vec_t &x, vec_t &x2) const {
1436  return base(x2);
1437  }
1438 
1439  };
1440 
1441  /** \brief A multi-dimensional Gaussian conditional probability
1442  density function
1443 
1444  This class is experimental.
1445 
1446  \todo This should be a symmetric conditional probability,
1447  i.e. \f$ P(x|y) = P(y|x) \f$. Test this.
1448  */
1449  template<class vec_t=boost::numeric::ublas::vector<double>,
1450  class mat_t=boost::numeric::ublas::matrix<double> >
1451  class prob_cond_mdim_gaussian : public prob_cond_mdim<vec_t> {
1452 
1453  protected:
1454 
1455  /// Cholesky decomposition
1456  mat_t chol;
1457 
1458  /// Inverse of the covariance matrix
1459  mat_t covar_inv;
1460 
1461  /// Normalization factor
1462  double norm;
1463 
1464  /// Number of dimensions
1465  size_t ndim;
1466 
1467  /// Temporary storage 1
1468  mutable vec_t q;
1469 
1470  /// Temporary storage 2
1471  mutable vec_t vtmp;
1472 
1473  /// Standard normal
1475 
1476  public:
1477 
1478  /** \brief Create an empty distribution
1479  */
1481  ndim=0;
1482  }
1483 
1484  /** \brief Create a distribution from the covariance matrix
1485  */
1486  prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar) {
1487  set(p_ndim,covar);
1488  }
1489 
1490  /// The dimensionality
1491  virtual size_t dim() const {
1492  return ndim;
1493  }
1494 
1495  /** \brief Set the covariance matrix for the distribution
1496  */
1497  void set(size_t p_ndim, mat_t &covar) {
1498  if (p_ndim==0) {
1499  O2SCL_ERR("Zero dimension in prob_cond_mdim_gaussian::set().",
1501  }
1502  ndim=p_ndim;
1503  norm=1.0;
1504  q.resize(ndim);
1505  vtmp.resize(ndim);
1506 
1507  // Perform the Cholesky decomposition
1508  chol=covar;
1509  o2scl_linalg::cholesky_decomp(ndim,chol);
1510 
1511  // Find the inverse
1512  covar_inv=chol;
1513  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1514 
1515  // Force chol to be lower triangular and compute the determinant
1516  double det=1.0;
1517  for(size_t i=0;i<ndim;i++) {
1518  det*=chol(i,i);
1519  for(size_t j=0;j<ndim;j++) {
1520  if (i<j) chol(i,j)=0.0;
1521  }
1522  }
1523  det*=det;
1524 
1525  // Compute normalization
1526  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt(det);
1527  }
1528 
1529  /// The conditional probability
1530  virtual double pdf(const vec_t &x, const vec_t &x2) const {
1531  if (ndim==0) {
1532  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1533  "pdf().",o2scl::exc_einval);
1534  }
1535  double ret=norm;
1536  for(size_t i=0;i<ndim;i++) q[i]=x2[i]-x[i];
1537  vtmp=prod(covar_inv,q);
1538  ret*=exp(-0.5*inner_prod(q,vtmp));
1539  return ret;
1540  }
1541 
1542  /// The log of the conditional probability
1543  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
1544  if (ndim==0) {
1545  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1546  "pdf().",o2scl::exc_einval);
1547  }
1548  double ret=log(norm);
1549  for(size_t i=0;i<ndim;i++) q[i]=x2[i]-x[i];
1550  vtmp=prod(covar_inv,q);
1551  ret+=-0.5*inner_prod(q,vtmp);
1552  return ret;
1553  }
1554 
1555  /// Sample the distribution
1556  virtual void operator()(const vec_t &x, vec_t &x2) const {
1557  if (ndim==0) {
1558  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1559  "operator().",o2scl::exc_einval);
1560  }
1561  for(size_t i=0;i<ndim;i++) q[i]=pdg();
1562  vtmp=prod(chol,q);
1563  for(size_t i=0;i<ndim;i++) x2[i]=x[i]+vtmp[i];
1564  return;
1565  }
1566 
1567  };
1568 
1569 #ifdef O2SCL_NEVER_DEFINED
1570  /** \brief A multidimensional normal distribution from
1571  a Gaussian process
1572 
1573  \future The linear algebra only works with ublas and is
1574  not optimized.
1575  */
1576  template<class vec_t=boost::numeric::ublas::vector<double>,
1577  class mat_t=boost::numeric::ublas::matrix<double>,
1578  class mat_col_t=boost::numeric::ublas::matrix_column<mat_t> >
1579  class prob_dens_mdim_gproc :
1580  public o2scl::prob_dens_mdim_gaussian<vec_t> {
1581 
1582  public:
1583 
1584  };
1585 #endif
1586 
1587 #ifndef DOXYGEN_NO_O2NS
1588 }
1589 #endif
1590 
1591 #endif
vec_t peak
Location of the peak.
double sigma_
Width parameter.
virtual double entropy() const
The inverse cumulative distribution function.
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...
rng_gsl rg
Internal random number generator.
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
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 void operator()(const vec_t &x, vec_t &x2) const
Sample the distribution.
virtual void operator()(const vec_t &x, vec_t &x2) const
Sample the distribution.
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=.
o2scl::prob_dens_gaussian pdg
Standard normal.
virtual double log_pdf(const vec_t &x, const vec_t &x2) const
The log of the conditional probability.
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.
vec_t vtmp
Temporary storage 2.
virtual double entropy() const
Inverse cumulative distribution function (from the lower tail)
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
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.
double random()
Return a random number in .
Definition: rng_gsl.h:82
double x0
The x coordinate of the centroid.
virtual double pdf(const vec_t &x, const vec_t &x2) const
The conditional probability.
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 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.
prob_dens_lognormal(const prob_dens_lognormal &pdg)
Copy constructor.
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 double log_pdf(const vec_t &x, const vec_t &x2) const
The log of the conditional probability.
virtual double log_metrop_hast(const vec_t &x_B, vec_t &x_A) const
Sample the distribution and return the log of the Metropolis-Hastings ratio.
A multidimensional distribution formed by the product of several one-dimensional distributions.
double norm
Normalization factor.
A one-dimensional probability density function.
virtual double pdf(double x) const
The normalized density.
virtual double entropy() const
The inverse cumulative distribution function.
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
The inverse cumulative distribution function.
virtual void operator()(vec_t &x) const
Sample the distribution.
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
void cholesky_decomp(const size_t M, mat_t &A)
Compute the in-place Cholesky decomposition of a symmetric positive-definite square matrix...
Definition: cholesky_base.h:61
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
A multi-dimensional probability density function.
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)
ubvector range
Vector specifying original histogram bins.
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.
virtual double pdf(const vec_t &x, const vec_t &x2) const
The conditional probability.
vec_t q
Temporary storage 1.
virtual double log_pdf(double x) const
The log of the normalized density.
A one-dimensional probability density over the positive real numbers.
vec_t q
Temporary storage 1.
virtual void operator()(const vec_t &x, vec_t &x2) const
Sample the distribution.
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 double pdf(const vec_t &x, const vec_t &x2) const
The conditional probability.
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.
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.
prob_dens_uniform(const prob_dens_uniform &pdg)
Copy constructor.
void set_seed(unsigned long int s)
Set the seed.
void set_seed(unsigned long int s)
Set the seed.
std::vector< double > u_low
Lower limits.
static const double x2[5]
Definition: inte_qng_gsl.h:66
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.
const ubvector & partial_sums()
Get reference to partial sums.
double sig_y
The y standard deviation.
double ll
Lower limit.
Lognormal density function.
virtual double log_pdf(const vec_t &x, const vec_t &x2) const
The log of the conditional probability.
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).