File wellmet/f_models.py changed (mode: 100644) (index 8330910..6f430b0) |
... |
... |
class SNorm: |
407 |
407 |
|
|
408 |
408 |
|
|
409 |
409 |
|
|
|
410 |
|
class Norm: |
|
411 |
|
__slots__ = ('G', '_U', '_pdf_G', 'mean', 'std', 'alpha') |
|
412 |
|
""" |
|
413 |
|
Normal distribution |
|
414 |
|
""" |
|
415 |
|
@classmethod |
|
416 |
|
def _fromdata(cls, G, U, pdf_G, mean, std, alpha=1): |
|
417 |
|
f = super().__new__(cls) |
|
418 |
|
f.G = G |
|
419 |
|
f._U = U |
|
420 |
|
f._pdf_G = pdf_G |
|
421 |
|
f.mean = mean |
|
422 |
|
f.std = std |
|
423 |
|
f.alpha = alpha |
|
424 |
|
return f |
|
425 |
|
|
|
426 |
|
@classmethod |
|
427 |
|
def from_G(cls, G, mean, std, alpha=1): |
|
428 |
|
f = super().__new__(cls) |
|
429 |
|
f.G = G |
|
430 |
|
f.mean = mean |
|
431 |
|
f.std = std |
|
432 |
|
f._U = None |
|
433 |
|
f._pdf_G = None |
|
434 |
|
f.alpha = alpha |
|
435 |
|
return f |
|
436 |
|
|
|
437 |
|
@classmethod |
|
438 |
|
def from_U(cls, U, mean, std, alpha=1): |
|
439 |
|
f = super().__new__(cls) |
|
440 |
|
f._U = U |
|
441 |
|
f.G = stats.norm.ppf(U) |
|
442 |
|
f.mean = mean |
|
443 |
|
f.std = std |
|
444 |
|
f._pdf_G = None |
|
445 |
|
f.alpha = alpha |
|
446 |
|
return f |
|
447 |
|
|
|
448 |
|
def __init__(self, mean, std, alpha=1): |
|
449 |
|
assert len(mean) == len(std) |
|
450 |
|
self.mean = mean |
|
451 |
|
self.std = std |
|
452 |
|
|
|
453 |
|
self.G = np.empty((0, len(mean)), dtype=float) |
|
454 |
|
self._U = None |
|
455 |
|
self._pdf_G = None |
|
456 |
|
self.alpha = alpha |
|
457 |
|
|
|
458 |
|
# nemusím duplikovat a dědit |
|
459 |
|
set_alpha = set_alpha |
|
460 |
|
|
|
461 |
|
def __repr__(self): |
|
462 |
|
return "%s(%s, %s, %s)"%('Norm', self.mean, self.std, repr(self.alpha)) |
|
463 |
|
|
|
464 |
|
def __str__(f): |
|
465 |
|
return "Norm sample: " + str(f.R) |
|
466 |
|
|
|
467 |
|
def __call__(self, ns=0): |
|
468 |
|
if ns: |
|
469 |
|
return Norm.from_G(np.random.randn(ns, self.nvar), self.mean, self.std, self.alpha) |
|
470 |
|
#else: |
|
471 |
|
return Norm(self.mean, self.std, self.alpha) |
|
472 |
|
|
|
473 |
|
def __len__(self): |
|
474 |
|
return len(self.G) |
|
475 |
|
|
|
476 |
|
def __getitem__(self, slice): |
|
477 |
|
G = np.atleast_2d(self.G[slice]) |
|
478 |
|
|
|
479 |
|
if self._U is not None: |
|
480 |
|
U = np.atleast_2d(self._get_U()[slice]) |
|
481 |
|
else: |
|
482 |
|
U = None |
|
483 |
|
|
|
484 |
|
if self._pdf_G is not None: |
|
485 |
|
pdf_G = np.atleast_1d(self._get_pdf()[slice]) |
|
486 |
|
else: |
|
487 |
|
pdf_G = None |
|
488 |
|
|
|
489 |
|
return Norm._fromdata(G, U, pdf_G, self.mean, self.std, self.alpha) |
|
490 |
|
|
|
491 |
|
|
|
492 |
|
def __getattr__(f, attr): |
|
493 |
|
if attr == 'f_model': |
|
494 |
|
return f |
|
495 |
|
elif attr == 'R': |
|
496 |
|
sample = f.G * f.std |
|
497 |
|
sample += f.mean |
|
498 |
|
return sample |
|
499 |
|
elif attr in ('Rn', 'GK'): |
|
500 |
|
return f.G |
|
501 |
|
elif attr in ('aR', 'aRn', 'aGK', 'aG'): |
|
502 |
|
return f.G * f.alpha |
|
503 |
|
elif attr in ('P', 'U'): |
|
504 |
|
return f._get_U() |
|
505 |
|
elif attr in ('aP', 'aU'): |
|
506 |
|
return f.U * f.alpha |
|
507 |
|
elif attr == 'nvar': |
|
508 |
|
return len(f.mean) |
|
509 |
|
elif attr == 'nsim': |
|
510 |
|
return len(f.G) |
|
511 |
|
elif attr == 'marginals': |
|
512 |
|
return [stats.norm(m, s) for m, s in zip(self.mean, self.std)] |
|
513 |
|
elif attr == 'cor': |
|
514 |
|
return np.diag([1 for __ in range(f.nvar)]) |
|
515 |
|
|
|
516 |
|
# hustoty |
|
517 |
|
# I'm considering to deprecate attribute access |
|
518 |
|
elif attr[:4] == 'pdf_': |
|
519 |
|
return f.pdf(attr[4:]) |
|
520 |
|
|
|
521 |
|
raise AttributeError(attr) |
|
522 |
|
|
|
523 |
|
|
|
524 |
|
def _get_U(self): |
|
525 |
|
_U = self._U |
|
526 |
|
G = self.G |
|
527 |
|
if _U is None: |
|
528 |
|
U = stats.norm.cdf(G) |
|
529 |
|
self._U = U |
|
530 |
|
return U |
|
531 |
|
|
|
532 |
|
len_U = len(_U) |
|
533 |
|
if len_U == len(G): |
|
534 |
|
return _U |
|
535 |
|
else: |
|
536 |
|
U = np.empty_like(G) |
|
537 |
|
U[:len_U] = _U |
|
538 |
|
U[len_U:] = stats.norm.cdf(G[len_U:]) |
|
539 |
|
self._U = U |
|
540 |
|
return U |
|
541 |
|
|
|
542 |
|
|
|
543 |
|
def _get_pdf(self): |
|
544 |
|
_pdf_G = self._pdf_G |
|
545 |
|
G = self.G |
|
546 |
|
if _pdf_G is None: |
|
547 |
|
pdf_G = self._calculate_pdf(G) |
|
548 |
|
self._pdf_G = pdf_G |
|
549 |
|
return pdf_G |
|
550 |
|
|
|
551 |
|
len_pdf = len(_pdf_G) |
|
552 |
|
nsim = len(G) |
|
553 |
|
if len_pdf == nsim: |
|
554 |
|
return _pdf_G |
|
555 |
|
else: |
|
556 |
|
pdf_G = np.empty(nsim, dtype=float) |
|
557 |
|
pdf_G[:len_pdf] = _pdf_G |
|
558 |
|
pdf_G[len_pdf:] = self._calculate_pdf(G[len_pdf:]) |
|
559 |
|
self._pdf_G = pdf_G |
|
560 |
|
return pdf_G |
|
561 |
|
|
|
562 |
|
|
|
563 |
|
@staticmethod |
|
564 |
|
def _calculate_pdf(G): |
|
565 |
|
return np.prod(stats.norm.pdf(G), axis=1) |
|
566 |
|
|
|
567 |
|
|
|
568 |
|
def _add_U(self, U_to_add): |
|
569 |
|
"č method přeppokladá, že self.G ještě není dotčen" |
|
570 |
|
n = len(U_to_add) |
|
571 |
|
nsim = self.nsim |
|
572 |
|
U = np.empty((nsim + n, self.nvar), dtype=float) |
|
573 |
|
U[nsim:] = U_to_add |
|
574 |
|
|
|
575 |
|
_U = self._U |
|
576 |
|
G = self.G |
|
577 |
|
if _U is None: |
|
578 |
|
U[:nsim] = stats.norm.cdf(G) |
|
579 |
|
else: |
|
580 |
|
len_U = len(_U) |
|
581 |
|
if len_U == len(G): |
|
582 |
|
U[:nsim] = _U |
|
583 |
|
else: |
|
584 |
|
U[:len_U] = _U |
|
585 |
|
U[len_U:nsim] = stats.norm.cdf(G[len_U:]) |
|
586 |
|
|
|
587 |
|
self._U = U |
|
588 |
|
|
|
589 |
|
|
|
590 |
|
def add_sample(f, input_sample, space='R'): |
|
591 |
|
""" |
|
592 |
|
Adds coordinates from input sample |
|
593 |
|
""" |
|
594 |
|
try: # is sample another f_model object? |
|
595 |
|
sample = getattr(input_sample, space) |
|
596 |
|
except AttributeError: |
|
597 |
|
# no, it is just coordinates array |
|
598 |
|
sample = input_sample |
|
599 |
|
|
|
600 |
|
sample = np.atleast_2d(sample) |
|
601 |
|
isim, ivar = np.shape(sample) # input sim, var |
|
602 |
|
if ivar != f.nvar: |
|
603 |
|
raise WellMetError('%sD data expected, but %sD sample given'% (f.nvar, ivar)) |
|
604 |
|
|
|
605 |
|
if isim == 0: |
|
606 |
|
return None |
|
607 |
|
|
|
608 |
|
nsim = f.nsim |
|
609 |
|
G = np.empty((nsim + isim, f.nvar), dtype=float) |
|
610 |
|
G[:nsim] = f.G |
|
611 |
|
|
|
612 |
|
# new piece of code "alpha"-related |
|
613 |
|
if space in ('aR', 'aRn', 'aGK', 'aG', 'aP', 'aU'): |
|
614 |
|
sample = sample / f.alpha |
|
615 |
|
space = space[1:] |
|
616 |
|
|
|
617 |
|
if space in ('Rn', 'GK', 'G'): |
|
618 |
|
G[nsim:] = sample |
|
619 |
|
f.G = G |
|
620 |
|
elif space == 'R': |
|
621 |
|
G[nsim:] = sample |
|
622 |
|
G[nsim:] -= f.mean |
|
623 |
|
G[nsim:] /= f.std |
|
624 |
|
f.G = G |
|
625 |
|
elif space in ('P', 'U'): |
|
626 |
|
#č _add_U() přeppokladá, že self.G ještě není dotčen |
|
627 |
|
f._add_U(sample) |
|
628 |
|
G[nsim:] = stats.norm.ppf(sample) |
|
629 |
|
f.G = G |
|
630 |
|
else: |
|
631 |
|
raise WellMetError('SNorm: unknown space %s' % space) |
|
632 |
|
|
|
633 |
|
|
|
634 |
|
#č tohle už není "drobná pomucka" |
|
635 |
|
#č to je jedná z funkcí, která běží 30% času |
|
636 |
|
def new_sample(f, input_sample=None, space='R', extend=False): |
|
637 |
|
""" |
|
638 |
|
Returns new f_model object with the same distribution and with coordinates from 'sample' taken |
|
639 |
|
""" |
|
640 |
|
if input_sample is None: |
|
641 |
|
return Norm(f.mean, f.std, f.alpha) |
|
642 |
|
|
|
643 |
|
try: # does sample is another f_model object? |
|
644 |
|
sample = getattr(input_sample, space) |
|
645 |
|
except AttributeError: |
|
646 |
|
# no, it is just coordinates array |
|
647 |
|
sample = input_sample |
|
648 |
|
|
|
649 |
|
sample = np.atleast_2d(sample) |
|
650 |
|
isim, ivar = np.shape(sample) # input sim, var |
|
651 |
|
if extend: |
|
652 |
|
to_extend = sample |
|
653 |
|
sample = np.zeros((isim, f.nvar)) |
|
654 |
|
sample[:,:ivar] = to_extend |
|
655 |
|
elif ivar != f.nvar: |
|
656 |
|
raise WellMetError('%sD data expected, but %sD sample given'% (f.nvar, ivar)) |
|
657 |
|
|
|
658 |
|
|
|
659 |
|
# new piece of code "alpha"-related |
|
660 |
|
if space in ('aR', 'aRn', 'aGK', 'aG', 'aP', 'aU'): |
|
661 |
|
sample = sample / f.alpha |
|
662 |
|
space = space[1:] |
|
663 |
|
|
|
664 |
|
|
|
665 |
|
if space in ('Rn', 'GK', 'G'): |
|
666 |
|
return Norm.from_G(sample, f.mean, f.std, f.alpha) |
|
667 |
|
elif space == 'R': |
|
668 |
|
sample = sample - f.mean |
|
669 |
|
sample /= f.std |
|
670 |
|
return Norm.from_G(sample, f.mean, f.std, f.alpha) |
|
671 |
|
elif space in ('P', 'U'): |
|
672 |
|
return Norm.from_U(sample, f.mean, f.std, f.alpha) |
|
673 |
|
else: |
|
674 |
|
raise WellMetError('SNorm: unknown space %s' % space) |
|
675 |
|
|
|
676 |
|
|
|
677 |
|
|
|
678 |
|
|
|
679 |
|
|
|
680 |
|
def pdf(f, space='R'): |
|
681 |
|
""" |
|
682 |
|
Returns own probability densities in the given space |
|
683 |
|
""" |
|
684 |
|
if space in ('Rn', 'GK', 'G'): |
|
685 |
|
return f._get_pdf() |
|
686 |
|
elif space == 'R': |
|
687 |
|
return f._get_pdf() / np.prod(f.std) |
|
688 |
|
elif space in ('P', 'U'): |
|
689 |
|
ones = np.ones(f.nsim) |
|
690 |
|
mask = np.all((f.U >= 0) & (f.U <= 1), axis=1) |
|
691 |
|
if not np.all(mask): |
|
692 |
|
ones[~mask] = 0 |
|
693 |
|
return ones |
|
694 |
|
elif space in ('aR', 'aRn', 'aGK', 'aG', 'aP', 'aU'): |
|
695 |
|
return f.pdf(space[1:]) / np.prod(f.alpha) |
|
696 |
|
|
|
697 |
|
else: |
|
698 |
|
raise WellMetError('Unknown space %s' % space) |
|
699 |
|
|
|
700 |
|
|
|
701 |
|
|
|
702 |
|
def sample_pdf(f, input_sample, space='R'): |
|
703 |
|
""" |
|
704 |
|
Calculates probability density for the given external sample in the given space. |
|
705 |
|
Function intended for the case no transformation needed. |
|
706 |
|
Otherwise new f_sample should be instanciated. |
|
707 |
|
""" |
|
708 |
|
|
|
709 |
|
# does sample is another f_model object? |
|
710 |
|
# кинлы со zase кулэ? |
|
711 |
|
try: |
|
712 |
|
sample = getattr(input_sample, space) |
|
713 |
|
except AttributeError: |
|
714 |
|
# no |
|
715 |
|
sample = np.atleast_2d(input_sample) |
|
716 |
|
|
|
717 |
|
|
|
718 |
|
if space in ('Rn', 'GK', 'G'): |
|
719 |
|
pdfs_R = stats.norm.pdf(sample) |
|
720 |
|
return np.prod(pdfs_R, axis=pdfs_R.ndim-1) |
|
721 |
|
|
|
722 |
|
elif space == 'R': |
|
723 |
|
sample = sample - f.mean |
|
724 |
|
sample /= f.std |
|
725 |
|
pdfs_R = stats.norm.pdf(sample) |
|
726 |
|
return np.prod(pdfs_R, axis=pdfs_R.ndim-1) / np.prod(f.std) |
|
727 |
|
|
|
728 |
|
elif space in ('P', 'U'): |
|
729 |
|
ones = np.ones(len(sample)) |
|
730 |
|
mask = np.all((sample >= 0) & (sample <= 1), axis=1) |
|
731 |
|
if not np.all(mask): |
|
732 |
|
ones[~mask] = 0 |
|
733 |
|
return ones |
|
734 |
|
|
|
735 |
|
# new piece of code "alpha"-related |
|
736 |
|
elif space in ('aR', 'aRn', 'aGK', 'aG', 'aP', 'aU'): |
|
737 |
|
return f.sample_pdf(sample / f.alpha, space=space[1:]) / np.prod(f.alpha) |
|
738 |
|
|
|
739 |
|
else: |
|
740 |
|
raise WellMetError('Unknown space %s' % space) |
410 |
741 |
|
|
411 |
742 |
|
|
412 |
743 |
|
|