kumaraswamy.py 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106
  1. # mypy: allow-untyped-defs
  2. from typing import Optional, Union
  3. import torch
  4. from torch import nan, Tensor
  5. from torch.distributions import constraints
  6. from torch.distributions.transformed_distribution import TransformedDistribution
  7. from torch.distributions.transforms import AffineTransform, PowerTransform
  8. from torch.distributions.uniform import Uniform
  9. from torch.distributions.utils import broadcast_all, euler_constant
  10. __all__ = ["Kumaraswamy"]
  11. def _moments(a, b, n):
  12. """
  13. Computes nth moment of Kumaraswamy using using torch.lgamma
  14. """
  15. arg1 = 1 + n / a
  16. log_value = torch.lgamma(arg1) + torch.lgamma(b) - torch.lgamma(arg1 + b)
  17. return b * torch.exp(log_value)
  18. class Kumaraswamy(TransformedDistribution):
  19. r"""
  20. Samples from a Kumaraswamy distribution.
  21. Example::
  22. >>> # xdoctest: +IGNORE_WANT("non-deterministic")
  23. >>> m = Kumaraswamy(torch.tensor([1.0]), torch.tensor([1.0]))
  24. >>> m.sample() # sample from a Kumaraswamy distribution with concentration alpha=1 and beta=1
  25. tensor([ 0.1729])
  26. Args:
  27. concentration1 (float or Tensor): 1st concentration parameter of the distribution
  28. (often referred to as alpha)
  29. concentration0 (float or Tensor): 2nd concentration parameter of the distribution
  30. (often referred to as beta)
  31. """
  32. arg_constraints = {
  33. "concentration1": constraints.positive,
  34. "concentration0": constraints.positive,
  35. }
  36. support = constraints.unit_interval
  37. has_rsample = True
  38. def __init__(
  39. self,
  40. concentration1: Union[Tensor, float],
  41. concentration0: Union[Tensor, float],
  42. validate_args: Optional[bool] = None,
  43. ) -> None:
  44. self.concentration1, self.concentration0 = broadcast_all(
  45. concentration1, concentration0
  46. )
  47. base_dist = Uniform(
  48. torch.full_like(self.concentration0, 0),
  49. torch.full_like(self.concentration0, 1),
  50. validate_args=validate_args,
  51. )
  52. transforms = [
  53. PowerTransform(exponent=self.concentration0.reciprocal()),
  54. AffineTransform(loc=1.0, scale=-1.0),
  55. PowerTransform(exponent=self.concentration1.reciprocal()),
  56. ]
  57. super().__init__(base_dist, transforms, validate_args=validate_args)
  58. def expand(self, batch_shape, _instance=None):
  59. new = self._get_checked_instance(Kumaraswamy, _instance)
  60. new.concentration1 = self.concentration1.expand(batch_shape)
  61. new.concentration0 = self.concentration0.expand(batch_shape)
  62. return super().expand(batch_shape, _instance=new)
  63. @property
  64. def mean(self) -> Tensor:
  65. return _moments(self.concentration1, self.concentration0, 1)
  66. @property
  67. def mode(self) -> Tensor:
  68. # Evaluate in log-space for numerical stability.
  69. log_mode = (
  70. self.concentration0.reciprocal() * (-self.concentration0).log1p()
  71. - (-self.concentration0 * self.concentration1).log1p()
  72. )
  73. log_mode[(self.concentration0 < 1) | (self.concentration1 < 1)] = nan
  74. return log_mode.exp()
  75. @property
  76. def variance(self) -> Tensor:
  77. return _moments(self.concentration1, self.concentration0, 2) - torch.pow(
  78. self.mean, 2
  79. )
  80. def entropy(self):
  81. t1 = 1 - self.concentration1.reciprocal()
  82. t0 = 1 - self.concentration0.reciprocal()
  83. H0 = torch.digamma(self.concentration0 + 1) + euler_constant
  84. return (
  85. t0
  86. + t1 * H0
  87. - torch.log(self.concentration1)
  88. - torch.log(self.concentration0)
  89. )