prob_dens_func.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2012-2017, 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 #include <random>
34 
35 #include <boost/numeric/ublas/vector.hpp>
36 
37 #include <o2scl/hist.h>
38 #include <o2scl/rng_gsl.h>
39 #include <o2scl/search_vec.h>
40 #include <o2scl/cholesky.h>
41 
42 #ifndef DOXYGEN_NO_O2NS
43 namespace o2scl {
44 #endif
45 
46  /** \brief A one-dimensional probability density function
47 
48  This class is experimental.
49 
50  \future Give functions for mean, median, mode, variance, etc?
51 
52  \comment
53  For now, there aren't any pure virtual functions,
54  since this causes problems in creating an
55  std::vector<prob_dens_func> object below (especially
56  with intel compilers)
57  \endcomment
58  */
60 
61  public:
62 
63  /// Sample from the specified density
64  virtual double operator()() const {
65  return 0.0;
66  }
67 
68  /// The normalized density
69  virtual double pdf(double x) const {
70  return 0.0;
71  }
72 
73  /// The log of the normalized density
74  virtual double log_pdf(double x) const {
75  return 0.0;
76  }
77 
78  /// The cumulative distribution function (from the lower tail)
79  virtual double cdf(double x) const {
80  return 0.0;
81  }
82 
83  /// The inverse cumulative distribution function
84  virtual double invert_cdf(double cdf) const {
85  return 0.0;
86  }
87 
88  /// Entropy of the distribution (\f$ - \int f \ln f \f$ )
89  virtual double entropy() const {
90  return 0.0;
91  }
92 
93  };
94 
95  /** \brief A one-dimensional Gaussian probability density
96 
97  The distribution
98  \f[
99  P(x)=\frac{1}{\sigma \sqrt{2 \pi}}
100  e^{-\frac{\left(x-x_0\right)^2}{2\sigma^2}}
101  \f]
102 
103  This class is experimental.
104  */
106 
107  protected:
108 
109  /** \brief Central value
110  */
111  double cent_;
112 
113  /** \brief Width parameter
114 
115  A value of -1 indicates it is yet unspecified.
116  */
117  double sigma_;
118 
119  /// Base GSL random number generator
121 
122  public:
123 
124  /** \brief Create a standard normal distribution
125  */
127  cent_=0.0;
128  sigma_=1.0;
129  }
130 
131  /** \brief Create a Gaussian distribution with width \c sigma
132 
133  The value of \c sigma must be larger than zero.
134  */
135  prob_dens_gaussian(double cent, double sigma) {
136  if (sigma<0.0) {
137  O2SCL_ERR2("Tried to create a Gaussian dist. with sigma",
138  "<0 in prob_dens_gaussian::prob_dens_gaussian().",
139  exc_einval);
140  }
141  cent_=cent;
142  sigma_=sigma;
143  }
144 
145  virtual ~prob_dens_gaussian() {
146  }
147 
148  /// Copy constructor
150  cent_=pdg.cent_;
151  sigma_=pdg.sigma_;
152  r=pdg.r;
153  }
154 
155  /// Copy constructor with operator=
157  // Check for self-assignment
158  if (this!=&pdg) {
159  cent_=pdg.cent_;
160  sigma_=pdg.sigma_;
161  r=pdg.r;
162  }
163  return *this;
164  }
165 
166  /// Set the seed
167  void set_seed(unsigned long int s) {
168  r.set_seed(s);
169  return;
170  }
171 
172  /// Set the center
173  void set_center(double cent) {
174  cent_=cent;
175  return;
176  }
177 
178  /// Set the Gaussian width (must be positive)
179  void set_sigma(double sigma) {
180  if (sigma<0.0) {
181  O2SCL_ERR2("Tried to set negative sigma ",
182  "in prob_dens_gaussian::prob_dens_gaussian().",
183  exc_einval);
184  }
185  sigma_=sigma;
186  return;
187  }
188 
189  /// Get the center
190  double mean() {
191  return cent_;
192  }
193 
194  /// Get the Gaussian width
195  double stddev() {
196  if (sigma_<0.0) {
197  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
198  "get_sigma().",exc_einval);
199  }
200  return sigma_;
201  }
202 
203  /// Sample from the specified density
204  virtual double operator()() const {
205  if (sigma_<0.0) {
206  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
207  "operator().",exc_einval);
208  }
209  return cent_+gsl_ran_gaussian(&r,sigma_);
210  }
211 
212  /// The normalized density
213  virtual double pdf(double x) const {
214  if (sigma_<0.0) {
215  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
216  "pdf().",exc_einval);
217  }
218  return gsl_ran_gaussian_pdf(x-cent_,sigma_);
219  }
220 
221  /// The log of the normalized density
222  virtual double log_pdf(double x) const {
223  if (sigma_<0.0) {
224  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
225  "pdf().",exc_einval);
226  }
227  return log(gsl_ran_gaussian_pdf(x-cent_,sigma_));
228  }
229 
230  /// The cumulative distribution function (from the lower tail)
231  virtual double cdf(double x) const {
232  if (sigma_<0.0) {
233  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
234  "cdf().",exc_einval);
235  }
236  return gsl_cdf_gaussian_P(x-cent_,sigma_);
237  }
238 
239  /// The inverse cumulative distribution function
240  virtual double invert_cdf(double in_cdf) const {
241  if (sigma_<0.0) {
242  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
243  "invert_cdf().",exc_einval);
244  }
245  if (in_cdf<0.0 || in_cdf>1.0) {
246  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
247  "prob_dens_gaussian::invert_cdf().",exc_einval);
248  }
249  return gsl_cdf_gaussian_Pinv(in_cdf,sigma_)+cent_;
250  }
251 
252  /// The inverse cumulative distribution function
253  virtual double entropy() const {
254  if (sigma_<0.0) {
255  O2SCL_ERR2("Width not set in prob_dens_gaussian::",
256  "invert_cdf().",exc_einval);
257  }
258  return log(2.0*o2scl_const::pi*exp(1.0)*sigma_*sigma_);
259  }
260 
261  };
262 
263  /** \brief A one-dimensional probability density over
264  a finite range
265 
266  This class is experimental.
267  */
269 
270  public:
271 
272  /// Lower limit of the range
273  virtual double lower_limit() const=0;
274 
275  /// Uower limit of the range
276  virtual double upper_limit() const=0;
277 
278  };
279 
280  /** \brief A uniform one-dimensional probability density
281  over a finite range
282 
283  A flat distribution given by \f$ P(x)=1/(b-a) \f$ for \f$ a<x<b
284  \f$, where \f$ a \f$ is the lower limit and \f$ b \f$ is the
285  upper limit.
286 
287  This class is experimental.
288  */
290 
291  protected:
292 
293  /// Lower limit
294  double ll;
295 
296  /// Upper limit
297  double ul;
298 
299  /// The GSL random number generator
301 
302  public:
303 
304  /** \brief Create a blank uniform distribution
305  */
307  ll=1.0;
308  ul=0.0;
309  }
310 
311  /** \brief Create a uniform distribution from \f$ a<x<b \f$
312  */
313  prob_dens_uniform(double a, double b) {
314  // Ensure a<b
315  if (a>b) {
316  double tmp=a;
317  a=b;
318  b=tmp;
319  }
320  ll=a;
321  ul=b;
322  }
323 
324  virtual ~prob_dens_uniform() {
325  }
326 
327  /// Copy constructor
329  ll=pdg.ll;
330  ul=pdg.ul;
331  }
332 
333  /// Copy constructor with operator=
335  // Check for self-assignment
336  if (this==&pdg) return *this;
337  ll=pdg.ll;
338  ul=pdg.ul;
339  return *this;
340  }
341 
342  /// Set the seed
343  void set_seed(unsigned long int s) {
344  r.set_seed(s);
345  return;
346  }
347 
348  /** \brief Set the limits of the uniform distribution
349  */
350  void set_limits(double a, double b) {
351  // Ensure a<b
352  if (a>b) {
353  double tmp=a;
354  a=b;
355  b=tmp;
356  }
357  ll=a;
358  ul=b;
359  return;
360  }
361 
362  /// Lower limit of the range
363  virtual double lower_limit() const {
364  if (ll>ul) {
365  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
366  "lower_limit().",exc_einval);
367  }
368  return ll;
369  }
370 
371  /// Uower limit of the range
372  virtual double upper_limit() const {
373  if (ll>ul) {
374  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
375  "upper_limit().",exc_einval);
376  }
377  return ul;
378  }
379 
380  /// Operator from the specified density
381  virtual double operator()() const {
382  if (ll>ul) {
383  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
384  "operator().",exc_einval);
385  }
386  return gsl_ran_flat(&r,ll,ul);
387  }
388 
389  /// The normalized density
390  virtual double pdf(double x) const {
391  if (ll>ul) {
392  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
393  "pdf().",exc_einval);
394  }
395  if (x<ll || x>ul) return 0.0;
396  return gsl_ran_flat_pdf(x,ll,ul);
397  }
398 
399  /// The log of the normalized density
400  virtual double log_pdf(double x) const {
401  if (ll>ul) {
402  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
403  "pdf().",exc_einval);
404  }
405  if (x<ll || x>ul) return 0.0;
406  return log(gsl_ran_flat_pdf(x,ll,ul));
407  }
408 
409  /// The cumulative distribution function (from the lower tail)
410  virtual double cdf(double x) const {
411  if (ll>ul) {
412  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
413  "cdf().",exc_einval);
414  }
415  if (x<ll) return 0.0;
416  if (x>ul) return 1.0;
417  return gsl_cdf_flat_P(x,ll,ul);
418  }
419 
420  /// The inverse cumulative distribution function
421  virtual double invert_cdf(double in_cdf) const {
422  if (ll>ul) {
423  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
424  "invert_cdf().",exc_einval);
425  }
426  if (in_cdf<0.0 || in_cdf>1.0) {
427  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
428  "prob_dens_uniform::invert_cdf().",exc_einval);
429  }
430  return gsl_cdf_flat_Pinv(in_cdf,ll,ul);
431  }
432 
433  /// The inverse cumulative distribution function
434  virtual double entropy() const {
435  if (ll>ul) {
436  O2SCL_ERR2("Limits not set in prob_dens_uniform::",
437  "entropy().",exc_einval);
438  }
439  return log(ul-ll);
440  }
441 
442  };
443 
444  /** \brief A one-dimensional probability density over the
445  positive real numbers
446 
447  This class is experimental.
448  */
450 
451  };
452 
453  /** \brief Lognormal density function
454 
455  The distribution
456  \f[
457  P(x)=\frac{1}{x \sigma \sqrt{2 \pi}}
458  \exp \left[-\frac{\left(\ln x-\mu\right)^2}{2\sigma^2}\right]
459  \f]
460 
461  This class is experimental.
462  */
464 
465  protected:
466 
467  /** \brief Width parameter
468 
469  A value of -1 indicates it is yet unspecified.
470  */
471  double sigma_;
472 
473  /** \brief Central value
474 
475  A value of -1 indicates it is yet unspecified.
476  */
477  double mu_;
478 
479  /// The GSL random number generator
481 
482  public:
483 
484  /** \brief Create a blank lognormal distribution
485  */
487  sigma_=-1.0;
488  mu_=0.0;
489  }
490 
491  /** \brief Create lognormal distribution with mean parameter \c mu
492  and width parameter \c sigma
493 
494  The value of \c sigma must be larger than zero.
495  */
496  prob_dens_lognormal(double mu, double sigma) {
497  if (sigma<0.0) {
498  O2SCL_ERR2("Tried to create log normal dist. with mu or sigma",
499  "<0 in prob_dens_lognormal::prob_dens_lognormal().",
500  exc_einval);
501  }
502  mu_=mu;
503  sigma_=sigma;
504  }
505 
506  virtual ~prob_dens_lognormal() {
507  }
508 
509  /// Copy constructor
511  mu_=pdg.mu_;
512  sigma_=pdg.sigma_;
513  }
514 
515  /// Copy constructor with operator=
517  // Check for self-assignment
518  if (this==&pdg) return *this;
519  mu_=pdg.mu_;
520  sigma_=pdg.sigma_;
521  return *this;
522  }
523 
524  /** \brief Set the maximum and width of the lognormal distribution
525  */
526  void set_mu_sigma(double mu, double sigma) {
527  if (sigma<0.0) {
528  O2SCL_ERR2("Tried to set mu or sigma negative",
529  "in prob_dens_lognormal::prob_dens_lognormal().",
530  exc_einval);
531  }
532  mu_=mu;
533  sigma_=sigma;
534  }
535 
536  /// Set the seed
537  void set_seed(unsigned long int s) {
538  r.set_seed(s);
539  }
540 
541  /// Sample from the specified density
542  virtual double operator()() const {
543  return gsl_ran_lognormal(&r,mu_,sigma_);
544  }
545 
546  /// The normalized density
547  virtual double pdf(double x) const {
548  if (x<0.0) {
549  return 0.0;
550  }
551  return gsl_ran_lognormal_pdf(x,mu_,sigma_);
552  }
553 
554  /// The log of the normalized density
555  virtual double log_pdf(double x) const {
556  if (x<0.0) {
557  return 0.0;
558  }
559  return log(gsl_ran_lognormal_pdf(x,mu_,sigma_));
560  }
561 
562  /// The cumulative distribution function (from the lower tail)
563  virtual double cdf(double x) const {
564  if (x<0.0) {
565  return 0.0;
566  }
567  return gsl_cdf_lognormal_P(x,mu_,sigma_);
568  }
569 
570  /// The inverse cumulative distribution function
571  virtual double invert_cdf(double in_cdf) const {
572  if (in_cdf<0.0 || in_cdf>1.0) {
573  O2SCL_ERR2("Requested cdf inverse outside of [0,1] in ",
574  "prob_dens_lognormal::invert_cdf().",exc_einval);
575  }
576  return gsl_cdf_lognormal_Pinv(in_cdf,mu_,sigma_);
577  }
578 
579  /// The inverse cumulative distribution function
580  virtual double entropy() const {
581  if (sigma_<0.0) {
582  O2SCL_ERR2("Parameters not set in prob_dens_lognormal::",
583  "entropy().",exc_einval);
584  }
585  return 0.5+0.5*log(2.0*o2scl_const::pi*sigma_*sigma_)+mu_;
586  }
587 
588  };
589 
590  /** \brief Probability density function based on a histogram
591 
592  This class is experimental.
593  */
595 
596  public:
597 
599 
600  protected:
601 
602  /// Search through the partial sums
604 
605  /// Number of original histogram bins
606  size_t n;
607 
608  /** \brief Normalized partial sum of histogram bins
609 
610  This vector has size \ref n plus one.
611  */
612  ubvector sum;
613 
614  /** \brief Vector specifying original histogram bins
615 
616  This vector has size \ref n plus one.
617  */
618  ubvector range;
619 
620  /// Random number generator
621  mutable rng_gsl rng;
622 
623  public:
624 
625  prob_dens_hist();
626 
627  ~prob_dens_hist();
628 
629  /// Initialize with histogram \c h
630  void init(hist &h);
631 
632  /// Generate a sample
633  virtual double operator()() const;
634 
635  /// Lower limit of the range
636  virtual double lower_limit() const;
637 
638  /// Uower limit of the range
639  virtual double upper_limit() const;
640 
641  /// The normalized density
642  virtual double pdf(double x) const;
643 
644  /// The log of the normalized density
645  virtual double log_pdf(double x) const;
646 
647  /// Cumulative distribution function (from the lower tail)
648  virtual double cdf(double x) const;
649 
650  /// Inverse cumulative distribution function (from the lower tail)
651  virtual double invert_cdf(double x) const;
652 
653  /// Inverse cumulative distribution function (from the lower tail)
654  virtual double entropy() const {
655  return 0.0;
656  }
657 
658  };
659 
660  /** \brief A multi-dimensional probability density function
661 
662  This class is experimental.
663  */
664  template<class vec_t=boost::numeric::ublas::vector<double> >
666 
667  public:
668 
669  /// Return the dimensionality
670  virtual size_t dim() const {
671  return 0;
672  }
673 
674  /// The normalized density
675  virtual double pdf(const vec_t &x) const {
676  return 0.0;
677  }
678 
679  /// The log of the normalized density
680  virtual double log_pdf(const vec_t &x) const {
681  return 0.0;
682  }
683 
684  /// Sample the distribution
685  virtual void operator()(vec_t &x) const {
686  return;
687  }
688 
689  };
690 
691  /** \brief A multidimensional distribution formed by the product
692  of several one-dimensional distributions
693  */
694  template<class vec_t=boost::numeric::ublas::vector<double> >
695  class prob_dens_mdim_factor : public prob_dens_mdim<vec_t> {
696 
697  protected:
698 
699  /// Vector of one-dimensional distributions
700  std::vector<prob_dens_func> list;
701 
702  public:
703 
704  prob_dens_mdim_factor(std::vector<prob_dens_func> &p_list) {
705  list=p_list;
706  }
707 
708  /// Return the dimensionality
709  virtual size_t dim() const {
710  return list.size();
711  }
712 
713  /// The normalized density
714  virtual double pdf(const vec_t &x) const {
715  double ret=1.0;
716  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
717  return ret;
718  }
719 
720  /// The log of the normalized density
721  virtual double log_pdf(const vec_t &x) const {
722  double ret=1.0;
723  for(size_t i=0;i<list.size();i++) ret*=list[i].pdf(x[i]);
724  return log(ret);
725  }
726 
727  /// Sample the distribution
728  virtual void operator()(vec_t &x) const {
729  for(size_t i=0;i<list.size();i++) x[i]=list[i]();
730  return;
731  }
732 
733  };
734 
735  /** \brief A multi-dimensional Gaussian probability density function
736 
737  This class is experimental.
738  */
739  template<class vec_t=boost::numeric::ublas::vector<double>,
740  class mat_t=boost::numeric::ublas::matrix<double> >
741  class prob_dens_mdim_gaussian : public prob_dens_mdim<vec_t> {
742 
743  protected:
744 
745  /// Cholesky decomposition
746  mat_t chol;
747 
748  /// Inverse of the covariance matrix
749  mat_t covar_inv;
750 
751  /// Location of the peak
752  vec_t peak;
753 
754  /// Normalization factor
755  double norm;
756 
757  /// Number of dimensions
758  size_t ndim;
759 
760  /// Temporary storage 1
761  mutable vec_t q;
762 
763  /// Temporary storage 2
764  mutable vec_t vtmp;
765 
766  /// Standard normal
768 
769  public:
770 
771  /// The dimensionality
772  virtual size_t dim() const {
773  return ndim;
774  }
775 
777  ndim=0;
778  }
779 
780  /** \brief Create a distribution from the covariance matrix
781  */
782  prob_dens_mdim_gaussian(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
783  set(p_ndim,p_peak,covar);
784  }
785 
786  /** \brief Set the peak and covariance matrix for the distribution
787  */
788  void set(size_t p_ndim, vec_t &p_peak, mat_t &covar) {
789  if (p_ndim==0) {
790  O2SCL_ERR("Zero dimension in prob_dens_mdim_gaussian::set().",
792  }
793  ndim=p_ndim;
794  norm=1.0;
795  peak.resize(ndim);
796  for(size_t i=0;i<ndim;i++) peak[i]=p_peak[i];
797  q.resize(ndim);
798  vtmp.resize(ndim);
799 
800  // Perform the Cholesky decomposition of the covariance matrix
801  chol=covar;
803 
804  // Find the inverse
805  covar_inv=chol;
806  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
807 
808  // Force chol to be lower triangular and compute the determinant
809  double det=1.0;
810  for(size_t i=0;i<ndim;i++) {
811  det*=chol(i,i);
812  for(size_t j=0;j<ndim;j++) {
813  if (i<j) chol(i,j)=0.0;
814  }
815  }
816  det*=det;
817 
818  // Compute normalization
819  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt(det);
820  }
821 
822  /// The normalized density
823  virtual double pdf(const vec_t &x) const {
824  if (ndim==0) {
825  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
826  "pdf().",o2scl::exc_einval);
827  }
828  double ret=norm;
829  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
830  vtmp=prod(covar_inv,q);
831  ret*=exp(-0.5*inner_prod(q,vtmp));
832  return ret;
833  }
834 
835  /// The log of the normalized density
836  virtual double log_pdf(const vec_t &x) const {
837  if (ndim==0) {
838  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
839  "pdf().",o2scl::exc_einval);
840  }
841  double ret=log(norm);
842  for(size_t i=0;i<ndim;i++) q[i]=x[i]-peak[i];
843  vtmp=prod(covar_inv,q);
844  ret+=-0.5*inner_prod(q,vtmp);
845  return ret;
846  }
847 
848  /// Sample the distribution
849  virtual void operator()(vec_t &x) const {
850  if (ndim==0) {
851  O2SCL_ERR2("Distribution not set in prob_dens_mdim_gaussian::",
852  "operator().",o2scl::exc_einval);
853  }
854  for(size_t i=0;i<ndim;i++) q[i]=pdg();
855  vtmp=prod(chol,q);
856  for(size_t i=0;i<ndim;i++) x[i]=peak[i]+vtmp[i];
857  return;
858  }
859 
860  };
861 
862  /** \brief A multi-dimensional conditional probability density function
863 
864  This class is experimental.
865  */
866  template<class vec_t=boost::numeric::ublas::vector<double> >
868 
869  public:
870 
871  /// The dimensionality
872  virtual size_t dim() const {
873  return 0;
874  }
875 
876  /// The normalized density
877  virtual double pdf(const vec_t &x, const vec_t &x2) const=0;
878 
879  /// The log of the normalized density
880  virtual double log_pdf(const vec_t &x, const vec_t &x2) const=0;
881 
882  /// Sample the distribution
883  virtual void operator()(const vec_t &x, vec_t &x2) const=0;
884 
885  /** \brief Sample the distribution and return the
886  log of the Metropolis-Hastings ratio
887  */
888  virtual double metrop_hast(const vec_t &x, vec_t &x2) const {
889  operator()(x,x2);
890  return log_pdf(x,x2)-log_pdf(x2,x);
891  }
892 
893  };
894 
895  /** \brief A constrained random walk in the shape of
896  a hypercube
897 
898  \comment
899  I had previously used std::uniform_real_distribution
900  instead of rng_gsl, but this caused problems with
901  intel compilers.
902  \endcomment
903  */
904  template<class vec_t=boost::numeric::ublas::vector<double> >
905  class prob_cond_mdim_rand_walk : public prob_cond_mdim<vec_t> {
906 
907  protected:
908 
909  /** \brief Desc
910  */
911  std::random_device rd;
912 
913  /** \brief Desc
914  */
915  std::vector<double> u_step;
916 
917  /** \brief Desc
918  */
919  std::vector<double> u_low;
920 
921  /** \brief Desc
922  */
923  std::vector<double> u_high;
924 
925  double d_pdf;
926 
927  /** \brief Desc
928  */
929  mutable rng_gsl rg;
930 
931  public:
932 
934  }
935 
936  template<class=vec_t> prob_cond_mdim_rand_walk
937  (vec_t &step, vec_t &low, vec_t &high) {
938  d_pdf=1.0;
939  for(size_t i=0;i<step.size();i++) {
940  u_step.push_back(step[i]);
941  if (low[i]>high[i]) {
942  double dtemp=low[i];
943  low[i]=high[i];
944  high[i]=dtemp;
945  }
946  u_low.push_back(low[i]);
947  u_high.push_back(high[i]);
948  d_pdf/=high[i]-low[i];
949  }
950  }
951 
952  /// The dimensionality
953  virtual size_t dim() const {
954  return u_step.size();
955  }
956 
957  /// The normalized density
958  virtual double pdf(const vec_t &x, const vec_t &x2) const {
959  return d_pdf;
960  }
961 
962  /// The log of the normalized density
963  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
964  return log(d_pdf);
965  }
966 
967  /// Sample the distribution
968  virtual void operator()(const vec_t &x, vec_t &x2) const {
969  size_t nv=u_step.size();
970  for(size_t i=0;i<nv;i++) {
971  while (x2[i]<u_low[i] || x2[i]>u_high[i]) {
972  x2[i]=x[i]+u_step[i]*(rg.random()*2.0-1.0);
973  }
974  }
975  return;
976  }
977 
978  };
979 
980  /** \brief A multi-dimensional conditional probability density function
981  independent of the input
982 
983  This class is experimental.
984  */
985  template<class vec_t=boost::numeric::ublas::vector<double> >
987 
988  protected:
989 
990  prob_dens_mdim<vec_t> &base;
991 
992  public:
993 
994  prob_cond_mdim_invar(prob_dens_mdim<vec_t> &out) : base(out) {
995  }
996 
997  /// The dimensionality
998  virtual size_t dim() const {
999  return base.dim();
1000  }
1001 
1002  /// The normalized density
1003  virtual double pdf(const vec_t &x, const vec_t &x2) const {
1004  return base.pdf(x2);
1005  }
1006 
1007  /// The log of the normalized density
1008  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
1009  return base.log_pdf(x2);
1010  }
1011 
1012  /// Sample the distribution
1013  virtual void operator()(const vec_t &x, vec_t &x2) const {
1014  return base(x2);
1015  }
1016 
1017  };
1018 
1019  /** \brief A multi-dimensional Gaussian conditional probability
1020  density function
1021 
1022  This class is experimental.
1023  */
1024  template<class vec_t=boost::numeric::ublas::vector<double>,
1025  class mat_t=boost::numeric::ublas::matrix<double> >
1026  class prob_cond_mdim_gaussian : public prob_cond_mdim<vec_t> {
1027 
1028  protected:
1029 
1030  /// Cholesky decomposition
1031  mat_t chol;
1032 
1033  /// Inverse of the covariance matrix
1034  mat_t covar_inv;
1035 
1036  /// Normalization factor
1037  double norm;
1038 
1039  /// Number of dimensions
1040  size_t ndim;
1041 
1042  /// Temporary storage 1
1043  mutable vec_t q;
1044 
1045  /// Temporary storage 2
1046  mutable vec_t vtmp;
1047 
1048  /// Standard normal
1050 
1051  public:
1052 
1053  /** \brief Create an empty distribution
1054  */
1056  ndim=0;
1057  }
1058 
1059  /** \brief Create a distribution from the covariance matrix
1060  */
1061  prob_cond_mdim_gaussian(size_t p_ndim, mat_t &covar) {
1062  set(p_ndim,covar);
1063  }
1064 
1065  /// The dimensionality
1066  virtual size_t dim() const {
1067  return ndim;
1068  }
1069 
1070  /** \brief Set the covariance matrix for the distribution
1071  */
1072  void set(size_t p_ndim, mat_t &covar) {
1073  if (p_ndim==0) {
1074  O2SCL_ERR("Zero dimension in prob_cond_mdim_gaussian::set().",
1076  }
1077  ndim=p_ndim;
1078  norm=1.0;
1079  q.resize(ndim);
1080  vtmp.resize(ndim);
1081 
1082  // Perform the Cholesky decomposition
1083  chol=covar;
1084  o2scl_linalg::cholesky_decomp(ndim,chol);
1085 
1086  // Find the inverse
1087  covar_inv=chol;
1088  o2scl_linalg::cholesky_invert<mat_t>(ndim,covar_inv);
1089 
1090  // Force chol to be lower triangular and compute the determinant
1091  double det=1.0;
1092  for(size_t i=0;i<ndim;i++) {
1093  det*=chol(i,i);
1094  for(size_t j=0;j<ndim;j++) {
1095  if (i<j) chol(i,j)=0.0;
1096  }
1097  }
1098  det*=det;
1099 
1100  // Compute normalization
1101  norm=pow(2.0*o2scl_const::pi,-((double)ndim)/2.0)/sqrt(det);
1102  }
1103 
1104  /// The normalized density
1105  virtual double pdf(const vec_t &x, const vec_t &x2) const {
1106  if (ndim==0) {
1107  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1108  "pdf().",o2scl::exc_einval);
1109  }
1110  double ret=norm;
1111  for(size_t i=0;i<ndim;i++) q[i]=x2[i]-x[i];
1112  vtmp=prod(covar_inv,q);
1113  ret*=exp(-0.5*inner_prod(q,vtmp));
1114  return ret;
1115  }
1116 
1117  /// The log of the normalized density
1118  virtual double log_pdf(const vec_t &x, const vec_t &x2) const {
1119  if (ndim==0) {
1120  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1121  "pdf().",o2scl::exc_einval);
1122  }
1123  double ret=log(norm);
1124  for(size_t i=0;i<ndim;i++) q[i]=x2[i]-x[i];
1125  vtmp=prod(covar_inv,q);
1126  ret+=-0.5*inner_prod(q,vtmp);
1127  return ret;
1128  }
1129 
1130  /// Sample the distribution
1131  virtual void operator()(const vec_t &x, vec_t &x2) const {
1132  if (ndim==0) {
1133  O2SCL_ERR2("Distribution not set in prob_cond_mdim_gaussian::",
1134  "operator().",o2scl::exc_einval);
1135  }
1136  for(size_t i=0;i<ndim;i++) q[i]=pdg();
1137  vtmp=prod(chol,q);
1138  for(size_t i=0;i<ndim;i++) x2[i]=x[i]+vtmp[i];
1139  return;
1140  }
1141 
1142  };
1143 
1144 
1145 #ifndef DOXYGEN_NO_O2NS
1146 }
1147 #endif
1148 
1149 #endif
vec_t peak
Location of the peak.
double sigma_
Width parameter.
virtual double entropy() const
The inverse cumulative distribution function.
virtual double log_pdf(const vec_t &x, const vec_t &x2) const
The log of the normalized density.
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.
const double pi
Definition: constants.h:62
std::vector< double > u_step
Desc.
virtual double lower_limit() const
Lower limit of the range.
virtual double pdf(const vec_t &x) const
The normalized density.
virtual size_t dim() const
The dimensionality.
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 pdf(const vec_t &x, const vec_t &x2) const
The normalized density.
virtual double log_pdf(const vec_t &x, const vec_t &x2) const
The log of the normalized density.
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.
virtual size_t dim() const
The dimensionality.
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)
prob_dens_lognormal & operator=(const prob_dens_lognormal &pdg)
Copy constructor with operator=.
virtual void operator()(vec_t &x) const
Sample the distribution.
prob_dens_lognormal()
Create a blank lognormal distribution.
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
std::random_device rd
Desc.
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 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.
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.
virtual void operator()(const vec_t &x, vec_t &x2) const
Sample the distribution.
vec_t vtmp
Temporary storage 2.
virtual double invert_cdf(double in_cdf) const
The inverse cumulative distribution function.
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.
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
virtual double metrop_hast(const vec_t &x, vec_t &x2) const
Sample the distribution and return the log of the Metropolis-Hastings ratio.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
A multi-dimensional probability density function.
virtual double cdf(double x) const
The cumulative distribution function (from the lower tail)
ubvector range
Vector specifying original histogram bins.
A multi-dimensional conditional probability density function independent of the input.
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.
#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 log_pdf(const vec_t &x, const vec_t &x2) const
The log of the normalized density.
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 size_t dim() const
The dimensionality.
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.
virtual double pdf(const vec_t &x, const vec_t &x2) const
The normalized density.
virtual size_t dim() const
The dimensionality.
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.
virtual double pdf(const vec_t &x, const vec_t &x2) const
The normalized density.
Searching class for monotonic data with caching.
Definition: search_vec.h:79
prob_dens_lognormal(double mu, double sigma)
Create lognormal distribution with mean parameter mu and width parameter sigma.
std::vector< double > u_high
Desc.
double norm
Normalization factor.
A constrained random walk in the shape of a hypercube.
rng_gsl rng
Random number generator.
Probability density function based on a histogram.
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 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.
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.
double ll
Lower limit.
Lognormal density function.
std::vector< double > u_low
Desc.
double stddev()
Get the Gaussian width.
mat_t chol
Cholesky decomposition.
void set_center(double cent)
Set the center.
virtual void operator()(const vec_t &x, vec_t &x2) const
Sample the distribution.
prob_dens_uniform & operator=(const prob_dens_uniform &pdg)
Copy constructor with operator=.

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