@@ -379,374 +379,7 @@ def logpt_sum(rv_var: TensorVariable, rv_value: Optional[TensorVariable] = None,
379
379
return aet .sum (logpt (rv_var , rv_value , ** kwargs ))
380
380
381
381
382
- import aesara .tensor as aet
383
- import numpy as np
384
-
385
- from aesara import config
386
- from aesara .graph .basic import Variable , ancestors , clone_replace
387
- from aesara .graph .op import compute_test_value
388
- from aesara .tensor .random .op import Observed , RandomVariable
389
- from aesara .tensor .subtensor import AdvancedSubtensor , AdvancedSubtensor1 , Subtensor
390
- from aesara .tensor .var import TensorVariable
391
-
392
- from pymc3 .aesaraf import floatX
393
-
394
- PotentialShapeType = Union [
395
- int , np .ndarray , Tuple [Union [int , Variable ], ...], List [Union [int , Variable ]], Variable
396
- ]
397
-
398
-
399
- def _get_scaling (total_size , shape , ndim ):
400
- """
401
- Gets scaling constant for logp
402
-
403
- Parameters
404
- ----------
405
- total_size: int or list[int]
406
- shape: shape
407
- shape to scale
408
- ndim: int
409
- ndim hint
410
-
411
- Returns
412
- -------
413
- scalar
414
- """
415
- if total_size is None :
416
- coef = floatX (1 )
417
- elif isinstance (total_size , int ):
418
- if ndim >= 1 :
419
- denom = shape [0 ]
420
- else :
421
- denom = 1
422
- coef = floatX (total_size ) / floatX (denom )
423
- elif isinstance (total_size , (list , tuple )):
424
- if not all (isinstance (i , int ) for i in total_size if (i is not Ellipsis and i is not None )):
425
- raise TypeError (
426
- "Unrecognized `total_size` type, expected "
427
- "int or list of ints, got %r" % total_size
428
- )
429
- if Ellipsis in total_size :
430
- sep = total_size .index (Ellipsis )
431
- begin = total_size [:sep ]
432
- end = total_size [sep + 1 :]
433
- if Ellipsis in end :
434
- raise ValueError (
435
- "Double Ellipsis in `total_size` is restricted, got %r" % total_size
436
- )
437
- else :
438
- begin = total_size
439
- end = []
440
- if (len (begin ) + len (end )) > ndim :
441
- raise ValueError (
442
- "Length of `total_size` is too big, "
443
- "number of scalings is bigger that ndim, got %r" % total_size
444
- )
445
- elif (len (begin ) + len (end )) == 0 :
446
- return floatX (1 )
447
- if len (end ) > 0 :
448
- shp_end = shape [- len (end ) :]
449
- else :
450
- shp_end = np .asarray ([])
451
- shp_begin = shape [: len (begin )]
452
- begin_coef = [floatX (t ) / shp_begin [i ] for i , t in enumerate (begin ) if t is not None ]
453
- end_coef = [floatX (t ) / shp_end [i ] for i , t in enumerate (end ) if t is not None ]
454
- coefs = begin_coef + end_coef
455
- coef = aet .prod (coefs )
456
- else :
457
- raise TypeError (
458
- "Unrecognized `total_size` type, expected int or list of ints, got %r" % total_size
459
- )
460
- return aet .as_tensor (floatX (coef ))
461
-
462
-
463
- def change_rv_size (
464
- rv_var : TensorVariable ,
465
- new_size : PotentialShapeType ,
466
- expand : Optional [bool ] = False ,
467
- ) -> TensorVariable :
468
- """Change or expand the size of a `RandomVariable`.
469
-
470
- Parameters
471
- ==========
472
- rv_var
473
- The `RandomVariable` output.
474
- new_size
475
- The new size.
476
- expand:
477
- Whether or not to completely replace the `size` parameter in `rv_var`
478
- with `new_size` or simply prepend it to the existing `size`.
479
-
480
- """
481
- rv_node = rv_var .owner
482
- rng , size , dtype , * dist_params = rv_node .inputs
483
- name = rv_var .name
484
- tag = rv_var .tag
485
-
486
- if expand :
487
- new_size = tuple (np .atleast_1d (new_size )) + tuple (size )
488
-
489
- new_rv_node = rv_node .op .make_node (rng , new_size , dtype , * dist_params )
490
- rv_var = new_rv_node .outputs [- 1 ]
491
- rv_var .name = name
492
- for k , v in tag .__dict__ .items ():
493
- rv_var .tag .__dict__ .setdefault (k , v )
494
-
495
- if config .compute_test_value != "off" :
496
- compute_test_value (new_rv_node )
497
-
498
- return rv_var
499
-
500
-
501
- def rv_log_likelihood_args (
502
- rv_var : TensorVariable ,
503
- rv_value : Optional [TensorVariable ] = None ,
504
- transformed : Optional [bool ] = True ,
505
- ) -> Tuple [TensorVariable , TensorVariable ]:
506
- """Get a `RandomVariable` and its corresponding log-likelihood `TensorVariable` value.
507
-
508
- Parameters
509
- ==========
510
- rv_var
511
- A variable corresponding to a `RandomVariable`, whether directly or
512
- indirectly (e.g. an observed variable that's the output of an
513
- `Observed` `Op`).
514
- rv_value
515
- The measure-space input `TensorVariable` (i.e. "input" to a
516
- log-likelihood).
517
- transformed
518
- When ``True``, return the transformed value var.
519
-
520
- Returns
521
- =======
522
- The first value in the tuple is the `RandomVariable`, and the second is the
523
- measure-space variable that corresponds with the latter. The first is used
524
- to determine the log likelihood graph and the second is the "input"
525
- parameter to that graph. In the case of an observed `RandomVariable`, the
526
- "input" is actual data; in all other cases, it's just another
527
- `TensorVariable`.
528
-
529
- """
530
-
531
- if rv_value is None :
532
- if rv_var .owner and isinstance (rv_var .owner .op , Observed ):
533
- rv_var , rv_value = rv_var .owner .inputs
534
- elif hasattr (rv_var .tag , "value_var" ):
535
- rv_value = rv_var .tag .value_var
536
- else :
537
- return rv_var , None
538
-
539
- rv_value = aet .as_tensor_variable (rv_value )
540
-
541
- transform = getattr (rv_value .tag , "transform" , None )
542
- if transformed and transform :
543
- rv_value = transform .forward (rv_value )
544
-
545
- return rv_var , rv_value
546
-
547
-
548
- def rv_ancestors (graphs : List [TensorVariable ]) -> Generator [TensorVariable , None , None ]:
549
- """Yield the ancestors that are `RandomVariable` outputs for the given `graphs`."""
550
- for anc in ancestors (graphs ):
551
- if anc in graphs :
552
- continue
553
- if anc .owner and isinstance (anc .owner .op , RandomVariable ):
554
- yield anc
555
-
556
-
557
- def strip_observed (x : TensorVariable ) -> TensorVariable :
558
- """Return the `RandomVariable` term for an `Observed` node input; otherwise, return the input."""
559
- if x .owner and isinstance (x .owner .op , Observed ):
560
- return x .owner .inputs [0 ]
561
- else :
562
- return x
563
-
564
-
565
- def sample_to_measure_vars (graphs : List [TensorVariable ]) -> List [TensorVariable ]:
566
- """Replace `RandomVariable` terms in graphs with their measure-space counterparts."""
567
- replace = {}
568
- for anc in ancestors (graphs ):
569
- if anc .owner and isinstance (anc .owner .op , RandomVariable ):
570
- measure_var = getattr (anc .tag , "value_var" , None )
571
- if measure_var is not None :
572
- replace [anc ] = measure_var
573
-
574
- dist_params = clone_replace (graphs , replace = replace )
575
- return dist_params
576
-
577
-
578
- def logpt (
579
- rv_var : TensorVariable ,
580
- rv_value : Optional [TensorVariable ] = None ,
581
- jacobian : bool = True ,
582
- scaling : Optional [bool ] = True ,
583
- ** kwargs ,
584
- ) -> TensorVariable :
585
- """Create a measure-space (i.e. log-likelihood) graph for a random variable at a given point.
586
-
587
- The input `rv_var` determines which log-likelihood graph is used and
588
- `rv_value` is that graph's input parameter. For example, if `rv_var` is
589
- the output of a `NormalRV` `Op`, then the output is
590
- ``normal_log_pdf(rv_value)``.
591
-
592
- Parameters
593
- ==========
594
- rv_var
595
- The `RandomVariable` output that determines the log-likelihood graph.
596
- rv_value
597
- The input variable for the log-likelihood graph.
598
- jacobian
599
- Whether or not to include the Jacobian term.
600
- scaling
601
- A scaling term to apply to the generated log-likelihood graph.
602
-
603
- """
604
-
605
- rv_var , rv_value = rv_log_likelihood_args (rv_var , rv_value )
606
- rv_node = rv_var .owner
607
-
608
- if not rv_node :
609
- raise TypeError ("rv_var must be the output of a RandomVariable Op" )
610
-
611
- if not isinstance (rv_node .op , RandomVariable ):
612
-
613
- if isinstance (rv_node .op , (Subtensor , AdvancedSubtensor , AdvancedSubtensor1 )):
614
-
615
- raise NotImplementedError ("Missing value support is incomplete" )
616
-
617
- # "Flatten" and sum an array of indexed RVs' log-likelihoods
618
- rv_var , missing_values = rv_node .inputs
619
- rv_value = rv_var .tag .value_var
620
-
621
- missing_values = missing_values .data
622
- logp_var = aet .sum (
623
- [
624
- logpt (
625
- rv_var ,
626
- )
627
- for idx , missing in zip (
628
- np .ndindex (missing_values .shape ), missing_values .flatten ()
629
- )
630
- if missing
631
- ]
632
- )
633
- return logp_var
634
-
635
- return aet .zeros_like (rv_var )
636
-
637
- rng , size , dtype , * dist_params = rv_node .inputs
638
-
639
- dist_params = sample_to_measure_vars (dist_params )
640
-
641
- if jacobian :
642
- logp_var = _logp (rv_node .op , rv_value , * dist_params , ** kwargs )
643
- else :
644
- logp_var = _logp_nojac (rv_node .op , rv_value , * dist_params , ** kwargs )
645
-
646
- # Replace `RandomVariable` ancestors with their corresponding
647
- # log-likelihood input variables
648
- lik_replacements = [
649
- (v , v .tag .value_var )
650
- for v in ancestors ([logp_var ])
651
- if v .owner and isinstance (v .owner .op , RandomVariable ) and getattr (v .tag , "value_var" , None )
652
- ]
653
-
654
- (logp_var ,) = clone_replace ([logp_var ], replace = lik_replacements )
655
-
656
- if scaling :
657
- logp_var *= _get_scaling (
658
- getattr (rv_var .tag , "total_size" , None ), rv_value .shape , rv_value .ndim
659
- )
660
-
661
- if rv_var .name is not None :
662
- logp_var .name = "__logp_%s" % rv_var .name
663
-
664
- return logp_var
665
-
666
-
667
- @singledispatch
668
- def _logp (op , value , * dist_params , ** kwargs ):
669
- """Create a log-likelihood graph.
670
-
671
- This function dispatches on the type of `op`, which should be a subclass
672
- of `RandomVariable`. If you want to implement new log-likelihood graphs
673
- for a `RandomVariable`, register a new function on this dispatcher.
674
-
675
- """
676
- return aet .zeros_like (value )
677
-
678
-
679
- def logcdf (rv_var , rv_value , ** kwargs ):
680
- """Create a log-CDF graph."""
681
-
682
- rv_var , rv_value = rv_log_likelihood_args (rv_var , rv_value )
683
- rv_node = rv_var .owner
684
-
685
- if not rv_node :
686
- raise TypeError ()
687
-
688
- rng , size , dtype , * dist_params = rv_node .inputs
689
-
690
- dist_params = sample_to_measure_vars (dist_params )
691
-
692
- return _logcdf (rv_node .op , rv_value , * dist_params , ** kwargs )
693
-
694
-
695
- @singledispatch
696
- def _logcdf (op , value , * args , ** kwargs ):
697
- """Create a log-CDF graph.
698
-
699
- This function dispatches on the type of `op`, which should be a subclass
700
- of `RandomVariable`. If you want to implement new log-CDF graphs
701
- for a `RandomVariable`, register a new function on this dispatcher.
702
-
703
- """
704
- raise NotImplementedError ()
705
-
706
-
707
- def logp_nojac (rv_var , rv_value = None , ** kwargs ):
708
- """Create a graph of the log-likelihood that doesn't include the Jacobian."""
709
-
710
- rv_var , rv_value = rv_log_likelihood_args (rv_var , rv_value )
711
- rv_node = rv_var .owner
712
-
713
- if not rv_node :
714
- raise TypeError ()
715
-
716
- rng , size , dtype , * dist_params = rv_node .inputs
717
-
718
- dist_params = sample_to_measure_vars (dist_params )
719
-
720
- return _logp_nojac (rv_node .op , rv_value , ** kwargs )
721
-
722
-
723
- @singledispatch
724
- def _logp_nojac (op , value , * args , ** kwargs ):
725
- """Return the logp, but do not include a jacobian term for transforms.
726
-
727
- If we use different parametrizations for the same distribution, we
728
- need to add the determinant of the jacobian of the transformation
729
- to make sure the densities still describe the same distribution.
730
- However, MAP estimates are not invariant with respect to the
731
- parameterization, we need to exclude the jacobian terms in this case.
732
-
733
- This function should be overwritten in base classes for transformed
734
- distributions.
735
- """
736
- return logpt (op , value , * args , ** kwargs )
737
-
738
-
739
- def logpt_sum (rv_var : TensorVariable , rv_value : Optional [TensorVariable ] = None , ** kwargs ):
740
- """Return the sum of the logp values for the given observations.
741
-
742
- Subclasses can use this to improve the speed of logp evaluations
743
- if only the sum of the logp values is needed.
744
- """
745
- return aet .sum (logpt (rv_var , rv_value , ** kwargs ))
746
-
747
-
748
- # from pymc3.distributions import timeseries
749
- from pymc3 .distributions import shape_utils , transforms
382
+ from pymc3 .distributions import shape_utils , timeseries , transforms
750
383
from pymc3 .distributions .bart import BART
751
384
from pymc3 .distributions .bound import Bound
752
385
from pymc3 .distributions .continuous import (
0 commit comments