Coverage for pythia/sampler.py: 89%
517 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-08 17:13 +0000
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-08 17:13 +0000
1"""
2File: pythia/sampler.py
3Author: Nando Hegemann
4Gitlab: https://gitlab.com/Nando-Hegemann
5Description: Sampler classes for generating in random samples and PDF evaluations.
6SPDX-License-Identifier: LGPL-3.0-or-later OR Hippocratic-3.0-ECO-MEDIA-MIL
7"""
8from typing import Callable
9from abc import ABC, abstractmethod, abstractproperty
10import warnings
11import numpy as np
12import scipy.stats
13from scipy.integrate import quad
14from scipy.special import gamma
15import pythia as pt
18class Sampler(ABC):
19 """Base class for all continuous samplers."""
21 # set abstract attributes
22 domain: np.ndarray
24 @abstractproperty
25 def dimension(self) -> int:
26 """Dimension of the ambient space."""
27 raise NotImplementedError
29 @abstractproperty
30 def mass(self):
31 """Mass of the sampler distribution.
33 The integral of the sampler distribution over the domain of
34 definition. If the density is normalised this value should be one.
35 """
36 raise NotImplementedError
38 @abstractproperty
39 def maximum(self) -> float:
40 """Maximum of the pdf."""
41 raise NotImplementedError
43 @abstractproperty
44 def mean(self) -> float | np.ndarray:
45 """Mean value of the pdf."""
46 raise NotImplementedError
48 @abstractproperty
49 def cov(self) -> float | np.ndarray:
50 """(Co)Variance of the pdf."""
51 raise NotImplementedError
53 @abstractmethod
54 def pdf(self, x: np.ndarray) -> np.ndarray:
55 """Density of the samplers distribution.
57 Computes the density of the samplers underlying distribution at the
58 given points `x`.
60 Parameters
61 ----------
62 x : array_like of shape (..., D)
63 list of points or single point. `D` is the objects dimension.
65 Returns
66 -------
67 :
68 Density values at the points.
69 """
70 raise NotImplementedError
72 @abstractmethod
73 def log_pdf(self, x: np.ndarray) -> np.ndarray:
74 """Log-density of the samplers distribution.
76 Computes the log-density of the samplers underlying distribution at the
77 given points `x`.
79 Parameters
80 ----------
81 x : array_like of shape (..., D)
82 list of points or single point. `D` is the objects dimension.
84 Returns
85 -------
86 :
87 Log-density values at the points.
88 """
89 raise NotImplementedError
91 @abstractmethod
92 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
93 """Gradient of log-density of the samplers distribution.
95 Computes the gradient of the log-density of the samplers underlying
96 distribution at the given points `x`.
98 Parameters
99 ----------
100 x : array_like of shape (..., D)
101 list of points or single point. `D` is the objects dimension.
103 Returns
104 -------
105 :
106 Gradient values of the log-density at the points with shape
107 (..., D).
108 """
109 raise NotImplementedError
111 @abstractmethod
112 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
113 """Hessian of log-density of the samplers distribution.
115 Computes the Hessian of the log-density of the samplers underlying
116 distribution at the given points `x`.
118 Parameters
119 ----------
120 x : array_like of shape (..., D)
121 list of points or single point. `D` is the objects dimension.
123 Returns
124 -------
125 :
126 Hessian values of the log-density at the points with shape
127 (..., D, D).
128 """
129 raise NotImplementedError
131 @abstractmethod
132 def sample(self, shape: list | tuple | np.ndarray) -> np.ndarray:
133 """Random values in a given shape.
135 Create an array of the given shape and populate it with random samples
136 from the samplers distribution.
138 Parameters
139 ----------
140 shape : array_like, optional
141 The dimensions of the returned array, should all be positive.
142 If no argument is given a single Python float is returned.
144 Returns
145 -------
146 :
147 Random values of specified shape.
148 """
149 raise NotImplementedError
152class UniformSampler(Sampler):
153 """Sampler for univariate uniformly distributed samples on given domain.
155 Parameters
156 ----------
157 domain : array_like
158 Interval of support of distribution.
159 """
161 def __init__(self, domain: list | tuple | np.ndarray) -> None:
162 """Initiate UniformSampler object."""
163 self.domain = np.reshape(np.array(domain), (-1, 2))
164 assert self.domain.shape == (1, 2)
165 self._length = self.domain[0, 1] - self.domain[0, 0]
167 @property
168 def dimension(self) -> int:
169 """Parameter dimension."""
170 return self.domain.shape[0]
172 @property
173 def mass(self) -> float:
174 """Mass of the PDF."""
175 return 1
177 @property
178 def maximum(self) -> float:
179 """Maximum value of the PDF."""
180 return 1 / self._length
182 @property
183 def mean(self) -> float:
184 """Mean value of the distribution."""
185 return 0.5 * (self.domain[0, 0] + self.domain[0, 1])
187 @property
188 def cov(self) -> float:
189 """(Co)Variance of the distribution."""
190 return 1 / 12 * (self.domain[0, 1] - self.domain[0, 0]) ** 2
192 @property
193 def var(self) -> float:
194 """Variance of the distribution."""
195 return self.cov
197 @property
198 def std(self) -> float:
199 """Standard deviation of the distribution."""
200 return np.sqrt(self.var)
202 def pdf(self, x: np.ndarray) -> np.ndarray:
203 """Evaluate uniform PDF.
205 Parameters
206 ----------
207 x : array_like
208 Evaluation points.
210 Returns
211 -------
212 :
213 Values of PDF evaluated in `x`.
214 """
215 return scipy.stats.uniform.pdf(x, loc=self.domain[0, 0], scale=self._length)
217 def log_pdf(self, x: np.ndarray) -> np.ndarray:
218 """Evaluate uniform log-PDF.
220 Parameters
221 ----------
222 x : array_like
223 Evaluation points.
225 Returns
226 -------
227 :
228 Values of log-PDF evaluated in `x`.
229 """
230 return scipy.stats.uniform.logpdf(x, loc=self.domain[0, 0], scale=self._length)
232 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
233 """Evaluate gradient of uniform log-PDF.
235 Parameters
236 ----------
237 x : array_like
238 Evaluation points.
240 Returns
241 -------
242 :
243 Values of gradient (vector valued) of log-PDF evaluated in `x`.
244 """
245 return np.zeros_like(x)
247 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
248 """Evaluate Hessian of uniform log-PDF.
250 Parameters
251 ----------
252 x : array_like
253 Evaluation points.
255 Returns
256 -------
257 :
258 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
259 """
260 return np.zeros_like(x)
262 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
263 """Draw samples from uniform distribution.
265 Parameters
266 ----------
267 shape : array_like
268 Shape of the samples.
270 Returns
271 -------
272 :
273 Random samples of specified shape.
274 """
275 return np.random.uniform(self.domain[0, 0], self.domain[0, 1], shape)
278class NormalSampler(Sampler):
279 """Sampler for univariate normally distributed samples.
281 Parameters
282 ----------
283 mean : float
284 Mean of the Gaussian distribution.
285 var : float
286 Variance of the Gaussian distribution.
287 """
289 def __init__(self, mean: float, var: float) -> None:
290 """Initiate NormalSampler object."""
291 self.domain = np.array([-np.inf, np.inf], ndmin=2) # shape is (1,2)
292 self._mean = mean # work-around to comply with ABC Sampler class
293 assert var >= 0
294 self.var = var
296 @property
297 def mass(self) -> float:
298 """Mass of the PDF."""
299 return 1
301 @property
302 def dimension(self) -> float:
303 """Dimension of the parameters."""
304 return self.domain.shape[0]
306 @property
307 def maximum(self) -> float:
308 """Maximum value of the PDF."""
309 return 1 / np.sqrt(2 * np.pi * self.var)
311 @property
312 def mean(self) -> float:
313 """Mean value of the distribution."""
314 return self._mean
316 @property
317 def cov(self) -> float:
318 """(Co)Variance of the distribution."""
319 return self.var
321 @property
322 def std(self) -> float:
323 """Standard deviation."""
324 return np.sqrt(self.var)
326 def pdf(self, x: np.ndarray) -> np.ndarray:
327 """Evaluate PDF.
329 Parameters
330 ----------
331 x : array_like
332 Evaluation points.
334 Returns
335 -------
336 :
337 Values of PDF evaluated in `x`.
338 """
339 return scipy.stats.norm.pdf(x, loc=self.mean, scale=np.sqrt(self.var))
341 def log_pdf(self, x: np.ndarray) -> np.ndarray:
342 """Evaluate log-PDF.
344 Parameters
345 ----------
346 x : array_like
347 Evaluation points.
349 Returns
350 -------
351 :
352 Values of log-PDF evaluated in `x`.
353 """
354 return scipy.stats.norm.logpdf(x, loc=self.mean, scale=np.sqrt(self.var))
356 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
357 """Evaluate gradient of log-PDF.
359 Parameters
360 ----------
361 x : array_like
362 Evaluation points.
364 Returns
365 -------
366 :
367 Values of gradient (vector valued) of log-PDF evaluated in `x`.
368 """
369 return -(x - self.mean) / self.var
371 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
372 """Evaluate Hessian of log-PDF.
374 Parameters
375 ----------
376 x : array_like
377 Evaluation points.
379 Returns
380 -------
381 :
382 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
383 """
384 return -1 / self.var * np.ones_like(x)
386 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
387 """Draw samples from distribution.
389 Parameters
390 ----------
391 shape : array_like
392 Shape of the samples.
394 Returns
395 -------
396 :
397 Random samples of specified shape.
398 """
399 return np.random.normal(self.mean, np.sqrt(self.var), shape)
402class GammaSampler(Sampler):
403 """Sampler for univariate Gamma distributed samples on given domain.
405 Parameters
406 ----------
407 domain : array_like
408 Supported domain of distribution.
409 alpha : float
410 Parameter for Gamma distribution.
411 beta : float
412 Parameter for Gamma distribution.
413 """
415 def __init__(
416 self, domain: list | tuple | np.ndarray, alpha: float, beta: float
417 ) -> None:
418 """Initiate GammaSampler object."""
419 self.domain = np.reshape(np.array(domain), (-1, 2))
420 assert self.domain.shape == (1, 2)
421 assert self.domain[0, 1] == np.inf
422 self.alpha = alpha
423 self.beta = beta
424 assert self.alpha > 0 and self.beta > 0
426 @property
427 def dimension(self) -> float:
428 """Dimension of the parameters."""
429 return 1
431 @property
432 def mass(self) -> float:
433 """Mass of the PDF."""
434 return 1
436 @property
437 def maximum(self) -> float:
438 """Maximum value of the PDF.
440 The maximum of the Gamma distribution is given by
442 .. math::
443 \\max_{x\\in[a,\\infty)} f(x) =
444 \\begin{cases}
445 \\infty & \\mbox{if } 0 < \\alpha < 1\\\\
446 \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} & \\mbox{if } \\alpha = 1\\\\
447 \\frac{\\beta^\\alpha}{\\Gamma(\\alpha)} \\Bigl(\\frac{\\alpha-1}{\\beta} \\Bigr)^{\\alpha-1} e^{1-\\alpha} & \\mbox{if } \\alpha > 1\\\\
448 \\end{cases}
449 """
450 if self.alpha < 1:
451 return np.inf
453 if self.alpha == 1:
454 return self.beta**self.alpha / gamma(self.alpha)
456 return (
457 self.beta**self.alpha
458 / gamma(self.alpha)
459 * ((self.alpha - 1) / self.beta) ** (self.alpha - 1)
460 * np.exp(-(self.alpha - 1))
461 )
463 @property
464 def mean(self) -> float:
465 """Mean value of the distribution."""
466 return self.alpha / self.beta + self.domain[0, 0]
468 @property
469 def cov(self) -> float:
470 """(Co)Variance of the distribution."""
471 return self.alpha / self.beta**2
473 @property
474 def var(self) -> float:
475 """Variance of the distribution."""
476 return self.cov
478 @property
479 def std(self) -> float:
480 """Standard deviation of the distribution."""
481 return np.sqrt(self.var)
483 def pdf(self, x: np.ndarray) -> np.ndarray:
484 """Evaluate PDF.
486 Parameters
487 ----------
488 x : array_like
489 Evaluation points.
491 Returns
492 -------
493 :
494 Values of PDF evaluated in `x`.
495 """
496 y = x - self.domain[0, 0]
497 return scipy.stats.gamma.pdf(y, a=self.alpha, scale=1.0 / self.beta)
499 def log_pdf(self, x: np.ndarray) -> np.ndarray:
500 """Evaluate log-PDF.
502 Parameters
503 ----------
504 x : array_like
505 Evaluation points.
507 Returns
508 -------
509 :
510 Values of log-PDF evaluated in `x`.
511 """
512 y = x - self.domain[0, 0]
513 return scipy.stats.gamma.logpdf(y, a=self.alpha, scale=1.0 / self.beta)
515 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
516 """Evaluate gradient of log-PDF.
518 .. note::
519 Not yet implemented.
521 Parameters
522 ----------
523 x : array_like
524 Evaluation points.
526 Returns
527 -------
528 :
529 Values of gradient (vector valued) of log-PDF evaluated in `x`.
530 """
531 raise NotImplementedError
533 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
534 """Evaluate Hessian of log-PDF.
536 .. note::
537 Not yet implemented.
539 Parameters
540 ----------
541 x : array_like
542 Evaluation points.
544 Returns
545 -------
546 :
547 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
548 """
549 raise NotImplementedError
551 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
552 """Draw samples from distribution.
554 Parameters
555 ----------
556 shape : array_like
557 Shape of the samples.
559 Returns
560 -------
561 :
562 Random samples of specified shape.
563 """
564 ret = np.random.gamma(self.alpha, 1.0 / self.beta, shape)
565 return ret + self.domain[0, 0]
568class BetaSampler(Sampler):
569 """Sampler for univariate Beta distributed samples on given domain.
571 Parameters
572 ----------
573 domain : array_like
574 Supported domain of distribution.
575 alpha : float
576 Parameter for Beta distribution.
577 beta : float
578 Parameter for Beta distribution.
579 """
581 def __init__(
582 self, domain: list | tuple | np.ndarray, alpha: float, beta: float
583 ) -> None:
584 """Initiate BetaSampler object."""
585 self.domain = np.reshape(np.array(domain), (-1, 2))
586 assert self.domain.shape == (1, 2)
587 self.length = self.domain[0, 1] - self.domain[0, 0]
588 self.alpha = alpha
589 self.beta = beta
590 assert self.alpha > 0 and self.beta > 0
592 @property
593 def dimension(self):
594 """Dimension of the parameters."""
595 return 1
597 @property
598 def mass(self):
599 """Mass of the PDF."""
600 return 1
602 @property
603 def maximum(self):
604 """Maximum value of the PDF.
606 The maximum of the Beta distribution is given by
608 .. math::
609 \\max_{x\\in[a,b]} f(x) =
610 \\begin{cases}
611 \\infty & \\mbox{if } 0 < \\alpha < 1 \\mbox{ or } 0 < \\beta < 1,\\\\
612 \\frac{1}{(b-a)B(\\alpha,\\beta)} & \\mbox{if } \\alpha = 1 \\mbox{ or } \\beta = 1,\\\\
613 \\frac{(\\alpha-1)^{\\alpha-1}(\\beta-1)^{\\beta-1}}{(\\alpha+\\beta-2)^{\\alpha+\\beta-2}(b-a)B(\\alpha,\\beta)} & \\mbox{if } \\alpha > 1 \\mbox{ and } \\beta > 1, \\\\
614 \\end{cases}
616 where :math:`B(\\alpha,\\beta)` denotes the Beta-function.
617 """
618 if self.alpha < 1 or self.beta < 1:
619 return np.inf
621 Beta = gamma(self.alpha) * gamma(self.beta) / gamma(self.alpha + self.beta)
622 if self.alpha == 1 or self.beta == 1:
623 return 1 / (Beta * self.length)
625 val = (
626 (self.alpha - 1) ** (self.alpha - 1)
627 * (self.beta - 1) ** (self.beta - 1)
628 / (self.alpha + self.beta - 2) ** (self.alpha + self.beta - 2)
629 )
630 return val / (Beta * self.length)
632 @property
633 def mean(self) -> float:
634 """Mean value of the distribution."""
635 shift = self.domain[0, 0]
636 scale = self.length
637 return scale * self.alpha / (self.alpha + self.beta) + shift
639 @property
640 def cov(self) -> float:
641 """(Co)Variance of the distribution."""
642 a = self.alpha
643 b = self.beta
644 return a * b / ((a + b) ** 2 * (a + b + 1)) * self.length
646 @property
647 def var(self) -> float:
648 """Variance of the distribution."""
649 return self.cov
651 @property
652 def std(self) -> float:
653 """Standard deviation of the distribution."""
654 return np.sqrt(self.var)
656 def pdf(self, x: np.ndarray) -> np.ndarray:
657 """Evaluate PDF.
659 Parameters
660 ----------
661 x : array_like
662 Evaluation points.
664 Returns
665 -------
666 :
667 Values of PDF evaluated in `x`.
668 """
669 y = pt.misc.shift_coord(x, self.domain.flatten(), [0, 1])
670 ret = scipy.stats.beta.pdf(y, a=self.alpha, b=self.beta)
671 return ret / self.length
673 def log_pdf(self, x: np.ndarray) -> np.ndarray:
674 """Evaluate log-PDF.
676 Parameters
677 ----------
678 x : array_like
679 Evaluation points.
681 Returns
682 -------
683 :
684 Values of log-PDF evaluated in `x`.
685 """
686 y = pt.misc.shift_coord(x, self.domain.flatten(), [0, 1])
687 ret = scipy.stats.beta.logpdf(y, a=self.alpha, b=self.beta)
688 return ret - np.log(self.length)
690 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
691 """Evaluate gradient of log-PDF.
693 .. note::
694 Not yet implemented.
696 Parameters
697 ----------
698 x : array_like
699 Evaluation points.
701 Returns
702 -------
703 :
704 Values of gradient (vector valued) of log-PDF evaluated in `x`.
705 """
706 raise NotImplementedError
708 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
709 """Evaluate Hessian of log-PDF.
711 .. note::
712 Not yet implemented.
714 Parameters
715 ----------
716 x : array_like
717 Evaluation points.
719 Returns
720 -------
721 :
722 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
723 """
724 raise NotImplementedError
726 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
727 """Draw samples from distribution.
729 Parameters
730 ----------
731 shape : array_like
732 Shape of the samples.
734 Returns
735 -------
736 :
737 Random samples of specified shape.
738 """
739 SAMPLE = np.random.beta(self.alpha, self.beta, shape)
740 return pt.misc.shift_coord(SAMPLE, [0, 1], self.domain.flatten())
743class WLSUnivariateSampler(Sampler):
744 """Sampler for univariate optimally distributed samples on given domain.
746 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}`
747 with :math:`y\\sim\\pi` and a finite subset
748 :math:`\\{P_j\\}_{j=0}^{d-1}` of an orthonormal polynomial basis of
749 :math:`L^2(\\Gamma,\\pi)`, the optimal weighted least-squares sampling
750 distribution for a function
751 :math:`u\\in\\operatorname{span}\\{P_j\\ \\vert\\ j=0,\\dots,d-1 \\}` reads
753 .. math::
754 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi
755 \\qquad\\mbox{with weight}\\qquad
756 w^{-1}(y) = \\frac{1}{d}\\sum_{j=0}^{d-1}\\vert P_j(y)\\vert^2.
758 Parameters
759 ----------
760 domain : array_like
761 Interval of support of distribution.
763 Notes
764 -----
765 To generate samples from the weighted least-squares distribution rejection
766 sampling is used. For certain basis functions it is possible to choose a
767 well-suited trial sampler for the rejection sampling, which can be enabled
768 via setting ``tsa=True``.
770 See Also
771 --------
772 pythia.sampler.WLSTensorSampler
774 References
775 ----------
776 The optimal weighted least-squares sampling is based on the results in
777 Cohen & Migliorati [1]_.
779 .. [1] Cohen, A. and Migliorati, G.,
780 “Optimal weighted least-squares methods”,
781 SMAI Journal of Computational Mathematics 3, 181-203 (2017).
782 """
784 def __init__(
785 self, param: pt.parameter.Parameter, deg: int, tsa: bool = True
786 ) -> None:
787 """Initiate WLSUnivariateSampler object."""
788 self.parameter = param
789 self.deg = deg
790 self._base_sampler = assign_sampler(param)
791 self.domain = self._base_sampler.domain
792 self._basis = pt.basis.univariate_basis([self.parameter], [self.deg])[0]
793 self._tsa = tsa
794 self._trial_sampler, self._bulk = self._compute_trial_sampler()
795 assert self.domain.shape == (1, 2)
796 # self._length = self.domain[0, 1] - self.domain[0, 0]
798 @property
799 def dimension(self) -> int:
800 """Parameter dimension."""
801 return self.domain.shape[0]
803 @property
804 def mass(self) -> float:
805 """Mass of the PDF."""
806 return 1
808 @property
809 def maximum(self) -> float:
810 """Maximum value of the PDF."""
811 if self.parameter.distribution == "uniform":
812 # alternatively return pdf in domain[0, 1]
813 return self.pdf(self.domain[0, 0])
814 elif self.parameter.distribution == "normal":
815 if self.deg % 2 == 0:
816 return self.pdf(self.parameter.mean)
817 else:
818 # Note: this can be improved
819 return get_maximum(self.pdf, self.domain)
820 elif self.parameter.distribution == "gamma":
821 # Note: this is probably not going to work for domain = [a, inf]
822 return get_maximum(self.pdf, self.domain)
823 elif self.parameter.distribution == "beta":
824 # Note: this can be improved
825 return get_maximum(self.pdf, self.domain)
826 return get_maximum(self.pdf, self.domain)
828 @property
829 def mean(self) -> float:
830 """Mean value of the distribution."""
831 mean, err = quad(
832 lambda x: x * self.pdf(x), self.domain[0, 0], self.domain[0, 1]
833 )
834 if err > 1e-8:
835 warnings.warn(f"quadrature error large: ({err})")
836 return mean
838 @property
839 def cov(self) -> float:
840 """(Co)Variance of the distribution."""
841 moment, err = quad(
842 lambda x: x**2 * self.pdf(x), self.domain[0, 0], self.domain[0, 1]
843 )
844 if err > 1e-8:
845 warnings.warn(f"quadrature error large: ({err})")
846 return moment - self.mean**2
848 @property
849 def var(self) -> float:
850 """Variance of the distribution."""
851 return self.cov
853 @property
854 def std(self) -> float:
855 """Standard deviation of the distribution."""
856 return np.sqrt(self.var)
858 def weight(self, x: np.ndarray | float | int) -> np.ndarray:
859 """Weights for the pdf.
861 Parameters
862 ----------
863 x : np.ndarray
864 Points the weight function is evaluated in.
866 Returns
867 -------
868 w : array_like
869 Weights of evaluation points `x`.
870 """
871 if isinstance(x, (float, int)):
872 x = np.array(x).reshape(1)
873 if x.ndim == 2 and x.shape[1] == 1:
874 x = np.array(x).reshape(-1)
875 assert x.ndim < 2
876 b_eval = np.sum([np.abs(p(x)) ** 2 for p in self._basis], axis=0)
877 return (self.deg + 1) / b_eval
879 def pdf(self, x: np.ndarray | float | int) -> np.ndarray:
880 """Evaluate uniform PDF.
882 Parameters
883 ----------
884 x : array_like
885 Evaluation points.
887 Returns
888 -------
889 :
890 Values of PDF evaluated in `x`.
891 """
892 if isinstance(x, (float, int)):
893 x = np.array(x).reshape(1)
894 if x.ndim == 2 and x.shape[1] == 1:
895 x = np.array(x).reshape(-1)
896 return self._base_sampler.pdf(x) / self.weight(x)
898 def log_pdf(self, x: np.ndarray | float | int) -> np.ndarray:
899 """Evaluate uniform log-PDF.
901 Parameters
902 ----------
903 x : array_like
904 Evaluation points.
906 Returns
907 -------
908 :
909 Values of log-PDF evaluated in `x`.
910 """
911 if isinstance(x, (float, int)):
912 x = np.array(x).reshape(1)
913 if x.ndim == 2 and x.shape[1] == 1:
914 x = np.array(x).reshape(-1)
915 return self._base_sampler.log_pdf(x) - np.log(self.weight(x))
917 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
918 """Evaluate gradient of uniform log-PDF.
920 Parameters
921 ----------
922 x : array_like
923 Evaluation points.
925 Returns
926 -------
927 :
928 Values of gradient (vector valued) of log-PDF evaluated in `x`.
929 """
930 raise NotImplementedError
932 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
933 """Evaluate Hessian of uniform log-PDF.
935 Parameters
936 ----------
937 x : array_like
938 Evaluation points.
940 Returns
941 -------
942 :
943 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
944 """
945 raise NotImplementedError
947 def _compute_trial_sampler(self) -> tuple[Sampler, float]:
948 """Trial sampler adaptation.
950 .. note::
951 TSA currently only available for uniform parameter distribution.
953 Parameters
954 ----------
955 tsa : bool
956 Adapt trial sampler or simply use uniform product sampler.
958 Returns
959 -------
960 trial_sampler : pt.sampler.Sampler
961 Trial sampler.
962 bulk : array_like
963 Domain estimate of the mass of the distribution.
964 """
965 if self._tsa is True and self.parameter.distribution == "uniform":
966 trial_sampler = BetaSampler(self.parameter.domain, 0.5, 0.5)
967 else:
968 trial_sampler = UniformSampler(self.parameter.domain)
969 bulk = get_maximum(lambda x: self.pdf(x) / trial_sampler.pdf(x), self.domain)
970 return trial_sampler, bulk
972 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
973 """Draw samples from weighted least-squares parameter distribution.
975 Parameters
976 ----------
977 shape : array_like
978 Shape of the samples.
980 Returns
981 -------
982 :
983 Random samples of specified shape.
984 """
985 if isinstance(shape, int):
986 shape = (shape,)
987 samples = rejection_sampling(
988 self.pdf, self._trial_sampler, self._bulk, self.dimension, shape
989 )
990 return samples
993class ProductSampler(Sampler):
994 """Tensor sampler for independent parameters.
996 Sampler for cartesian product samples of a list of (independent) univariate
997 samplers.
999 Parameters
1000 ----------
1001 sampler_list : list of `pythia.sampler.Sampler`
1002 list of (univariate) Sampler objects.
1003 """
1005 def __init__(self, sampler_list: list[Sampler]) -> None:
1006 """Initiate ProductSampler object."""
1007 self.samplers = list(sampler_list)
1008 self.domain = np.squeeze(np.array([s.domain for s in self.samplers]))
1009 # Add dimension if len(sampler_list) == 1.
1010 if self.domain.ndim < 2:
1011 self.domain.shape = 1, 2
1012 # sampler_list should contain 1D samplers only
1013 assert self.domain.shape == (len(sampler_list), 2)
1015 @property
1016 def dimension(self) -> int:
1017 """Dimension of the parameters."""
1018 return self.domain.shape[0]
1020 @property
1021 def mass(self) -> float:
1022 """Mass of the PDF."""
1023 return np.prod(np.array([s.mass for s in self.samplers]))
1025 @property
1026 def maximum(self) -> float:
1027 """Maximum value of the PDF."""
1028 return np.prod(np.array([s.maximum for s in self.samplers]))
1030 @property
1031 def mean(self) -> np.ndarray:
1032 """Mean of the PDF."""
1033 return np.array([s.mean for s in self.samplers])
1035 @property
1036 def cov(self) -> np.ndarray:
1037 """Covariance of the PDF."""
1038 return np.diag([s.cov for s in self.samplers])
1040 def weight(self, x: np.ndarray) -> np.ndarray:
1041 """Weights of the product PDF.
1043 Parameters
1044 ----------
1045 x : np.ndarray
1046 Evaluation points.
1048 Returns
1049 -------
1050 :
1051 Array of uniform weights for samples.
1052 """
1053 if x.ndim == 1:
1054 x.shape = 1, -1
1055 assert x.ndim == 2
1056 return np.ones(x.shape[0])
1058 def pdf(self, x: np.ndarray) -> np.ndarray:
1059 """Evaluate PDF.
1061 The PDF is given by the product of the univariate PDFs.
1063 Parameters
1064 ----------
1065 x : array_like
1066 Evaluation points.
1068 Returns
1069 -------
1070 :
1071 Values of PDF evaluated in `x`.
1072 """
1073 assert x.shape[-1] == self.dimension
1074 densities = [s.pdf for s in self.samplers]
1075 val = np.array([densities[jj](x[..., jj]) for jj in range(self.dimension)])
1076 return np.prod(val, axis=0)
1078 def log_pdf(self, x: np.ndarray) -> np.ndarray:
1079 """Evaluate log-PDF.
1081 The log-PDF is given by the sum of the univariate log-PDFs.
1083 Parameters
1084 ----------
1085 x : array_like
1086 Evaluation points.
1088 Returns
1089 -------
1090 :
1091 Values of log-PDF evaluated in `x`.
1092 """
1093 assert x.shape[-1] == self.dimension
1094 densities = [s.log_pdf for s in self.samplers]
1095 val = np.array([densities[jj](x[..., jj]) for jj in range(self.dimension)])
1096 return np.sum(val, axis=0)
1098 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1099 """Evaluate gradient of log-PDF.
1101 Parameters
1102 ----------
1103 x : array_like
1104 Evaluation points.
1106 Returns
1107 -------
1108 :
1109 Values of gradient (vector valued) of log-PDF evaluated in `x`.
1110 """
1111 grad_densities = [s.grad_x_log_pdf for s in self.samplers]
1112 return np.array(
1113 [grad_densities[jj](x[..., jj]) for jj in range(self.dimension)]
1114 ).T
1116 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1117 """Evaluate Hessian of log-PDF.
1119 Parameters
1120 ----------
1121 x : array_like
1122 Evaluation points.
1124 Returns
1125 -------
1126 :
1127 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
1128 """
1129 # create an (x.shape[0], self._dim, self._dim) tensor where each
1130 # (self._dim, self._dim) matrix is the identity
1131 eye = np.tile(np.expand_dims(np.eye(self.dimension), 0), (x.shape[0], 1, 1))
1133 # create an (x.shape[0],self._dim,1) tensor where each (x.shape[0],1,1)
1134 # subtensor is a diagonal entry of the hessian
1135 hess_densities = [s.hess_x_log_pdf for s in self.samplers]
1136 hess = np.expand_dims(
1137 np.array(
1138 [hess_densities[jj](x[..., jj]) for jj in range(self.dimension)]
1139 ).T,
1140 2,
1141 )
1142 return eye * hess
1144 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
1145 """Draw samples from distribution.
1147 Parameters
1148 ----------
1149 shape : array_like
1150 Shape of the samples.
1152 Returns
1153 -------
1154 :
1155 Random samples of specified shape.
1156 """
1157 if isinstance(shape, int):
1158 shape = (shape,)
1159 samples = [s.sample(shape) for s in self.samplers]
1160 return np.stack(samples, -1)
1163class ParameterSampler(Sampler):
1164 """Product sampler of given parameters.
1166 Parameters
1167 ----------
1168 params : list of `pythia.parameter.Parameter`
1169 list containing information of parameters.
1170 """
1172 def __init__(self, params: list[pt.parameter.Parameter]) -> None:
1173 """Initiate ParameterSampler object."""
1174 self.parameter = params
1175 assert isinstance(self.parameter, list)
1176 self.domain = np.array([param.domain for param in self.parameter])
1177 self._product_sampler = ProductSampler(
1178 [assign_sampler(param) for param in self.parameter]
1179 )
1181 @property
1182 def dimension(self) -> int:
1183 """Dimension of the parameters."""
1184 return self._product_sampler.dimension
1186 @property
1187 def mass(self) -> float:
1188 """Mass of the PDF."""
1189 return self._product_sampler.mass
1191 @property
1192 def maximum(self) -> float:
1193 """Maximum value of the PDF."""
1194 return self._product_sampler.maximum
1196 @property
1197 def mean(self) -> np.ndarray:
1198 """Mean of the PDF."""
1199 return self._product_sampler.mean
1201 @property
1202 def cov(self) -> np.ndarray:
1203 """Covariance of the PDF."""
1204 return self._product_sampler.cov
1206 def weight(self, x: np.ndarray) -> np.ndarray:
1207 """Weights of the parameter product PDF.
1209 Parameters
1210 ----------
1211 x : np.ndarray
1212 Evaluation points.
1214 Returns
1215 -------
1216 :
1217 Array of uniform weights for samples.
1218 """
1219 return self._product_sampler.weight(x)
1221 def pdf(self, x: np.ndarray) -> np.ndarray:
1222 """Evaluate PDF.
1224 Parameters
1225 ----------
1226 x : array_like
1227 Evaluation points.
1229 Returns
1230 -------
1231 :
1232 Values of PDF evaluated in `x`.
1233 """
1234 return self._product_sampler.pdf(x)
1236 def log_pdf(self, x: np.ndarray) -> np.ndarray:
1237 """Evaluate log-PDF.
1239 The log-PDF is given by the sum of the univariate log-PDFs.
1241 Parameters
1242 ----------
1243 x : array_like
1244 Evaluation points.
1246 Returns
1247 -------
1248 :
1249 Values of log-PDF evaluated in `x`.
1250 """
1251 return self._product_sampler.log_pdf(x)
1253 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1254 """Evaluate gradient of log-PDF.
1256 Parameters
1257 ----------
1258 x : array_like
1259 Evaluation points.
1261 Returns
1262 -------
1263 :
1264 Values of gradient (vector valued) of log-PDF evaluated in `x`.
1265 """
1266 return self._product_sampler.grad_x_log_pdf(x)
1268 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1269 """Evaluate Hessian of log-PDF.
1271 Parameters
1272 ----------
1273 x : array_like
1274 Evaluation points.
1276 Returns
1277 -------
1278 :
1279 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
1280 """
1281 return self._product_sampler.hess_x_log_pdf(x)
1283 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
1284 """Draw samples from distribution.
1286 Parameters
1287 ----------
1288 shape : array_like
1289 Shape of the samples.
1291 Returns
1292 -------
1293 :
1294 Random samples of specified shape.
1295 """
1296 return self._product_sampler.sample(shape)
1299class WLSSampler(Sampler):
1300 """Weighted Least-Squares sampler.
1302 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}^M`
1303 with :math:`y\\sim\\pi`, a set of multiindices
1304 :math:`\\Lambda\\subset\\mathbb{N}_0^M` and a finite subset
1305 :math:`\\{P_\\alpha\\}_{\\alpha\\in\\Lambda}` of an orthonormal polynomial
1306 basis of :math:`L^2(\\Gamma,\\pi)`, the optimal weighted least-squares
1307 sampling distribution for a function
1308 :math:`u\\in\\operatorname{span}\\{P_\\alpha\\ \\vert\\ \\alpha\\in\\Lambda\\}`
1309 reads
1311 .. math::
1312 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi
1313 \\qquad\\mbox{with weight}\\qquad
1314 w^{-1}(y) = \\frac{1}{\\vert\\Lambda\\vert}\\sum_{\\alpha\\in\\Lambda}\\vert P_\\alpha(y)\\vert^2,
1316 where :math:`\\vert\\Lambda\\vert` denotes the number of elements in
1317 :math:`\\Lambda`.
1319 Parameters
1320 ----------
1321 params : list of `pythia.parameter.Parameter`
1322 list of parameters.
1323 basis : list
1324 list of basis functions.
1325 tsa : bool, default=False
1326 Trial sampler adaptation. If True, a trial sampler is chosen on the
1327 distributions of parameters, if false a uniform trial sampler is used.
1329 Other Parameters
1330 ----------------
1331 trial_sampler : pythia.sampler.Sampler, default=None
1332 Trial sampler for rejection sampling. If `tsa` is true and either
1333 `trial_sampler` or `bulk` are `None`, the trial sampler is chosen
1334 automatically.
1335 bulk : float, defaul=None
1336 Scaling for trial sampler. If `tsa` is true and either
1337 `trial_sampler` or `bulk` are `None`, the trial sampler is chosen
1338 automatically.
1340 Notes
1341 -----
1342 To generate samples from the weighted least-squares distribution rejection
1343 sampling is used. For certain basis functions it is possible to choose a
1344 well-suited trial sampler for the rejection sampling, which can be enabled
1345 via setting ``tsa=True``.
1347 See Also
1348 --------
1349 pythia.sampler.WLSUnivariateSampler, pythia.sampler.WLSTensorSampler
1351 References
1352 ----------
1353 The optimal weighted least-squares sampling is based on the results of
1354 Cohen & Migliorati [1]_.
1355 """
1357 def __init__(
1358 self,
1359 params: list[pt.parameter.Parameter],
1360 basis: list[Callable],
1361 tsa: bool = True,
1362 trial_sampler: Sampler | None = None,
1363 bulk: float | None = None,
1364 ) -> None:
1365 """Initiate WLSSampler object."""
1366 self.parameter = params
1367 self._param_sampler = ParameterSampler(self.parameter)
1368 self.domain = self._param_sampler.domain
1369 self.basis = basis
1370 self._tsa = tsa
1371 if trial_sampler is None or bulk is None:
1372 self._trial_sampler, self._bulk = self._compute_trial_sampler()
1373 else:
1374 self._trial_sampler = trial_sampler
1375 self._bulk = bulk
1377 @property
1378 def dimension(self) -> int:
1379 """Dimension of the parameters."""
1380 return self._param_sampler.dimension
1382 @property
1383 def mass(self):
1384 """Mass of the PDF."""
1385 return 1
1387 @property
1388 def maximum(self) -> float:
1389 """Maximum value of the PDF."""
1390 return get_maximum(self.pdf, self.domain)
1392 @property
1393 def mean(self) -> np.ndarray:
1394 """Mean of the PDF."""
1395 raise NotImplementedError
1397 @property
1398 def cov(self) -> np.ndarray:
1399 """Covariance of the PDF."""
1400 raise NotImplementedError
1402 def weight(self, x: np.ndarray) -> np.ndarray:
1403 """Weights for the PDF.
1405 Parameters
1406 ----------
1407 x : array_like
1408 Points the weight function is evaluated in.
1410 Returns
1411 -------
1412 :
1413 weights of evaluation points `x`.
1414 """
1415 km = np.sum([abs(p(x)) ** 2 for p in self.basis], axis=0)
1416 return len(self.basis) / km
1418 def pdf(self, x: np.ndarray) -> np.ndarray:
1419 """Evaluate PDF.
1421 Parameters
1422 ----------
1423 x : array_like
1424 Evaluation points.
1426 Returns
1427 -------
1428 :
1429 Values of PDF evaluated in `x`.
1430 """
1431 assert x.shape[-1] == self.dimension
1432 return self._param_sampler.pdf(x) / self.weight(x)
1434 def log_pdf(self, x: np.ndarray) -> np.ndarray:
1435 """Evaluate log-PDF.
1437 The log-PDF is given by the sum of the univariate log-PDFs.
1439 Parameters
1440 ----------
1441 x : array_like
1442 Evaluation points.
1444 Returns
1445 -------
1446 :
1447 Values of log-PDF evaluated in `x`.
1448 """
1449 assert x.shape[-1] == self.dimension
1450 return self._param_sampler.log_pdf(x) - np.log(self.weight(x))
1452 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1453 """Evaluate gradient of log-PDF.
1455 Parameters
1456 ----------
1457 x : array_like
1458 Evaluation points.
1460 Returns
1461 -------
1462 :
1463 Values of gradient (vector valued) of log-PDF evaluated in `x`.
1464 """
1465 raise NotImplementedError
1467 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1468 """Evaluate Hessian of log-PDF.
1470 Parameters
1471 ----------
1472 x : array_like
1473 Evaluation points.
1475 Returns
1476 -------
1477 :
1478 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
1479 """
1480 raise NotImplementedError
1482 def _compute_trial_sampler(self) -> tuple[Sampler, float]:
1483 """Trial sampler adaptation.
1485 Parameters
1486 ----------
1487 tsa : bool
1488 Adapt trial sampler or simply use uniform product sampler.
1490 Returns
1491 -------
1492 trial_sampler : pt.sampler.Sampler
1493 Trial sampler.
1494 bulk : array_like
1495 Domain estimate of the mass of the distribution.
1496 """
1497 sampler_list = []
1498 for param in self.parameter:
1499 if self._tsa is True and param.distribution == "uniform":
1500 sampler_list.append(BetaSampler(param.domain, 0.5, 0.5))
1501 else:
1502 sampler_list.append(UniformSampler(param.domain))
1503 trial_sampler = ProductSampler(sampler_list)
1504 bulk = get_maximum(lambda x: self.pdf(x) / trial_sampler.pdf(x), self.domain)
1505 return trial_sampler, bulk
1507 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
1508 """Draw samples from distribution.
1510 Parameters
1511 ----------
1512 shape : array_like
1513 Shape of the samples.
1515 Returns
1516 -------
1517 :
1518 Random samples of specified shape.
1519 """
1520 if isinstance(shape, int):
1521 shape = (shape,)
1522 samples = rejection_sampling(
1523 self.pdf, self._trial_sampler, self._bulk, self.dimension, shape
1524 )
1525 return samples
1528class WLSTensorSampler(Sampler):
1529 """Weighted least-squares sampler for tensor multivariate basis.
1531 Given a stochastic variable :math:`y\\in\\Gamma\\subset\\mathbb{R}^M`
1532 with :math:`y\\sim\\pi=\\prod_{m=1}^M\\pi_m` for one dimensional densities
1533 :math:`\\pi_m`, a tensor set of multiindices
1534 :math:`\\Lambda=[d_1]\\times\\dots\\times[d_M]\\subset\\mathbb{N}_0^M`,
1535 where :math:`[d_m]=\\{0,\\dots,d_m-1\\}`, and a finite subset
1536 :math:`\\{P_\\alpha\\}_{\\alpha\\in\\Lambda}` of an orthonormal product
1537 polynomial basis of :math:`L^2(\\Gamma,\\pi)`, i.e.,
1538 :math:`P_\\alpha(y) = \\prod_{m=1}^M P_{\\alpha_m}(y_m)`, the optimal
1539 weighted least-squares sampling distribution for a function
1540 :math:`u\\in\\operatorname{span}\\{P_\\alpha\\ \\vert\\ \\alpha\\in\\Lambda\\}`
1541 reads
1543 .. math::
1544 \\mathrm{d}\\mu = w^{-1} \\mathrm{d}\\pi
1545 \\qquad\\mbox{with weight}\\qquad
1546 w^{-1}(y) = \\prod_{m=1}^M\\frac{1}{d_m}\\sum_{\\alpha_m\\in[d_m]}\\vert P_{\\alpha_m}(y_m)\\vert^2.
1548 Parameters
1549 ----------
1550 params : list of `pythia.parameter.Parameter`
1551 Parameter list.
1552 deg : list of int
1553 Polynomial degree of each component (same for all).
1554 tsa : bool, default=True
1555 Trial sampler adaptation. If True, a trial sampler is chosen on the
1556 distributions of parameters, if false a uniform trial sampler is
1557 used.
1559 Notes
1560 -----
1561 To generate samples from the weighted least-squares distribution rejection
1562 sampling is used. For certain basis functions it is possible to choose a
1563 well-suited trial sampler for the rejection sampling, which can be enabled
1564 via setting ``tsa=True``.
1566 See Also
1567 --------
1568 pythia.sampler.WLSUnivariateSampler
1570 References
1571 ----------
1572 The optimal weighted least-squares sampling is based on the results in
1573 Cohen & Migliorati [1]_.
1574 """
1576 def __init__(
1577 self, params: list[pt.parameter.Parameter], deg: list[int], tsa: bool = True
1578 ) -> None:
1579 """Initiate WLSTensorSampler object."""
1580 self.parameter = params
1581 self.deg = deg
1582 self.domain = np.array([param.domain for param in self.parameter])
1583 self._tsa = tsa
1584 # self._uBasis = pt.basis.univariate_basis(self.parameter, self.deg)
1585 # self._univariate_samplers = [assign_sampler(param)
1586 # for param in self.parameter]
1587 self._univariate_samplers = [
1588 WLSUnivariateSampler(p, d, self._tsa)
1589 for (p, d) in zip(self.parameter, self.deg)
1590 ]
1591 self._product_sampler = ProductSampler(self._univariate_samplers)
1593 @property
1594 def dimension(self) -> int:
1595 """Dimension of the parameters."""
1596 return self._product_sampler.dimension
1598 @property
1599 def mass(self) -> float:
1600 """Mass of the PDF."""
1601 return self._product_sampler.mass
1603 @property
1604 def maximum(self) -> float:
1605 """Maximum value of the PDF."""
1606 return self._product_sampler.maximum
1608 @property
1609 def mean(self) -> np.ndarray:
1610 """Mean of the PDF."""
1611 return self._product_sampler.mean
1613 @property
1614 def cov(self) -> np.ndarray:
1615 """Covariance of the PDF."""
1616 return self._product_sampler.cov
1618 def weight(self, x: np.ndarray) -> np.ndarray:
1619 """Weights for the PDF.
1621 Parameters
1622 ----------
1623 x : array_like
1624 Points the weight function is evaluated in.
1626 Returns
1627 -------
1628 :
1629 Weights of evaluation points `x`.
1630 """
1631 assert x.shape[-1] == self.dimension
1632 weights = [s.weight for s in self._univariate_samplers]
1633 val = np.array([weights[jj](x[..., jj]) for jj in range(self.dimension)])
1634 return np.prod(val, axis=0)
1636 def pdf(self, x: np.ndarray) -> np.ndarray:
1637 """Evaluate PDF.
1639 Parameters
1640 ----------
1641 x : array_like
1642 Evaluation points.
1644 Returns
1645 -------
1646 :
1647 Values of PDF evaluated in `x`.
1648 """
1649 return self._product_sampler.pdf(x)
1651 def log_pdf(self, x: np.ndarray) -> np.ndarray:
1652 """Evaluate log-PDF.
1654 The log-PDF is given by the sum of the univariate log-PDFs.
1656 Parameters
1657 ----------
1658 x : array_like
1659 Evaluation points.
1661 Returns
1662 -------
1663 :
1664 Values of log-PDF evaluated in `x`.
1665 """
1666 return self._product_sampler.log_pdf(x)
1668 def grad_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1669 """Evaluate gradient of log-PDF.
1671 Parameters
1672 ----------
1673 x : array_like
1674 Evaluation points.
1676 Returns
1677 -------
1678 :
1679 Values of gradient (vector valued) of log-PDF evaluated in `x`.
1680 """
1681 return self._product_sampler.grad_x_log_pdf(x)
1683 def hess_x_log_pdf(self, x: np.ndarray) -> np.ndarray:
1684 """Evaluate Hessian of log-PDF.
1686 Parameters
1687 ----------
1688 x : array_like
1689 Evaluation points.
1691 Returns
1692 -------
1693 :
1694 Values of Hessian (matrix valued) of log-PDF evaluated in `x`.
1695 """
1696 return self._product_sampler.hess_x_log_pdf(x)
1698 def sample(self, shape: int | list | tuple | np.ndarray) -> np.ndarray:
1699 """Draw samples from distribution.
1701 Parameters
1702 ----------
1703 shape : array_like
1704 Shape of the samples.
1706 Returns
1707 -------
1708 :
1709 Random samples of specified shape.
1710 """
1711 return self._product_sampler.sample(shape)
1714def rejection_sampling(
1715 pdf: Callable,
1716 trial_sampler: Sampler,
1717 scale: float,
1718 dimension: int,
1719 shape: int | list | tuple | np.ndarray,
1720) -> np.ndarray:
1721 """Draw samples from pdf by rejection sampling.
1723 Parameters
1724 ----------
1725 pdf : Callable
1726 Probability density the samples are generated from.
1727 trial_sampler : Sampler
1728 Trial sampler proposal samples are drawn from.
1729 scale : float
1730 Threshold parameter with ``pdf <= scale * trialSampler.pdf``
1731 dimension : int
1732 Dimension of the (input of the) pdf.
1733 shape : array_like
1734 Shape of the samples.
1736 Returns
1737 -------
1738 :
1739 Random samples of specified shape.
1740 """
1741 if isinstance(shape, int):
1742 shape = (shape,)
1743 size = np.prod(np.array(shape))
1744 trial_samples = trial_sampler.sample(size)
1745 samples = np.empty_like(trial_samples)
1746 pointer = 0
1747 while pointer < size:
1748 trial_samples = trial_sampler.sample(size)
1749 checker = np.random.rand(trial_samples.shape[0])
1750 is_valid = checker * max(scale, 1) * trial_sampler.pdf(trial_samples) <= pdf(
1751 trial_samples
1752 )
1753 valid_samples = trial_samples[is_valid]
1754 start = pointer
1755 end = min(pointer + valid_samples.shape[0], size)
1756 samples[start:end] = valid_samples[
1757 0 : min(valid_samples.shape[0], size - pointer)
1758 ]
1759 pointer = end
1760 if dimension > 1:
1761 return np.moveaxis(samples.T.reshape(dimension, *shape), 0, -1)
1762 return samples.reshape(*shape)
1765def constraint_sampling(
1766 sampler: Sampler,
1767 constraints: list[Callable],
1768 shape: int | list | tuple | np.ndarray,
1769) -> np.ndarray:
1770 """Draw samples according to algebraic constraints.
1772 Draw samples from target distribution and discard samples that do not
1773 satisfy the constraints.
1775 Parameters
1776 ----------
1777 sampler : Sampler
1778 Sampler to sample from.
1779 constraints : list of callable
1780 list of functions that return True if sample point satisfies the
1781 constraint.
1783 Returns
1784 -------
1785 :
1786 Samples drawn from sampler satisfying the constraints.
1788 Notes
1789 -----
1790 The constaints may lead to a non-normalized density function.
1791 """
1792 if isinstance(shape, int):
1793 shape = (shape,)
1794 samples = np.empty(sampler.sample(np.prod(shape)).shape)
1795 for j_sample in range(samples.shape[0]):
1796 proposal = sampler.sample((1,))
1797 met_constraints = np.all([c(proposal) for c in constraints])
1798 while not met_constraints:
1799 proposal = sampler.sample((1,))
1800 met_constraints = np.all([c(proposal) for c in constraints])
1801 samples[j_sample] = proposal
1802 return np.moveaxis(samples.T.reshape(sampler.dimension, *shape), 0, -1)
1805def assign_sampler(param: pt.parameter.Parameter) -> Sampler:
1806 """Assign a univariate sampler to the given parameter.
1808 Parameters
1809 ----------
1810 param : pythia.parameter.Parameter
1812 Returns
1813 -------
1814 :
1815 Univariate sampler.
1816 """
1817 if param.distribution == "uniform":
1818 return UniformSampler(param.domain)
1819 elif param.distribution == "normal":
1820 return NormalSampler(param.mean, param.var)
1821 elif param.distribution == "gamma":
1822 return GammaSampler(param.domain, param.alpha, param.beta)
1823 elif param.distribution == "beta":
1824 return BetaSampler(param.domain, param.alpha, param.beta)
1825 raise ValueError(f"unknown distribution '{param.distribution}'")
1828def get_maximum(f: Callable, domain: list | tuple | np.ndarray, n: int = 1000) -> float:
1829 """Compute essential maximum of function by point evaluations.
1831 Parameters
1832 ----------
1833 f : callable
1834 Function to evaluate. Needs to map from n-dim space to 1-dim space.
1835 domain : array_like
1836 Domain to evaluate function on.
1837 n : int, default=1000
1838 Number of function evaluations. Evaluations are done on a uniform grid
1839 in domain. Actual number of points may thus be a little greater.
1841 Returns
1842 -------
1843 :
1844 Approximation of maximum of function `f`.
1845 """
1846 if not isinstance(domain, np.ndarray):
1847 domain = np.array(domain)
1848 assert domain.shape[-1] == 2 and domain.ndim == 2
1849 n_points = int(np.ceil(np.power(n, 1 / domain.shape[0])))
1850 if domain.shape[0] > 1:
1851 eps = np.finfo(float).eps # circumvent inf on bdry
1852 x = pt.misc.cart_prod(
1853 [np.linspace(dom[0] + eps, dom[1] - eps, n_points) for dom in domain]
1854 )
1855 else:
1856 x = np.linspace(domain[0, 0], domain[0, 1], n_points).reshape(-1, 1)
1857 return np.max(np.abs(f(x)))
1860if __name__ == "__main__":
1861 pass