1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
################################################################################
# Copyright (C) 2014 Jaakko Luttinen
#
# This file is licensed under the MIT License.
################################################################################
"""
Module for the Poisson distribution node.
"""
import numpy as np
from scipy import special
from .expfamily import ExponentialFamily
from .expfamily import ExponentialFamilyDistribution
from .node import Moments
from .gamma import GammaMoments
from bayespy.utils import misc
class PoissonMoments(Moments):
"""
Class for the moments of Poisson variables
"""
dims = ( (), )
def compute_fixed_moments(self, x):
"""
Compute the moments for a fixed value
"""
# Make sure the values are integers in valid range
x = np.asanyarray(x)
if not misc.isinteger(x):
raise ValueError("Count not integer")
# Now, the moments are just the counts
return [x]
@classmethod
def from_values(cls, x):
"""
Return the shape of the moments for a fixed value.
The realizations are scalars, thus the shape of the moment is ().
"""
return cls()
class PoissonDistribution(ExponentialFamilyDistribution):
"""
Class for the VMP formulas of Poisson variables.
"""
def compute_message_to_parent(self, parent, index, u, u_lambda):
"""
Compute the message to a parent node.
"""
if index == 0:
m0 = -1
m1 = np.copy(u[0])
return [m0, m1]
else:
raise ValueError("Index out of bounds")
def compute_phi_from_parents(self, u_lambda, mask=True):
"""
Compute the natural parameter vector given parent moments.
"""
l = u_lambda[0]
logl = u_lambda[1]
phi0 = logl
return [phi0]
def compute_moments_and_cgf(self, phi, mask=True):
"""
Compute the moments and :math:`g(\phi)`.
"""
u0 = np.exp(phi[0])
u = [u0]
g = -u0
return (u, g)
def compute_cgf_from_parents(self, u_lambda):
"""
Compute :math:`\mathrm{E}_{q(p)}[g(p)]`
"""
l = u_lambda[0]
g = -l
return g
def compute_fixed_moments_and_f(self, x, mask=True):
"""
Compute the moments and :math:`f(x)` for a fixed value.
"""
# Check the validity of x
x = np.asanyarray(x)
if not misc.isinteger(x):
raise ValueError("Values must be integers")
if np.any(x < 0):
raise ValueError("Values must be positive")
# Compute moments
u0 = np.copy(x)
u = [u0]
# Compute f(x)
f = -special.gammaln(x+1)
return (u, f)
def random(self, *phi, plates=None):
"""
Draw a random sample from the distribution.
"""
return np.random.poisson(np.exp(phi[0]), size=plates)
class Poisson(ExponentialFamily):
"""
Node for Poisson random variables.
The node uses Poisson distribution:
.. math::
p(x) = \mathrm{Poisson}(x|\lambda)
where :math:`\lambda` is the rate parameter.
Parameters
----------
l : gamma-like node or scalar or array
:math:`\lambda`, rate parameter
See also
--------
Gamma, Exponential
"""
dims = ( (), )
_moments = PoissonMoments()
_parent_moments = [GammaMoments()]
_distribution = PoissonDistribution()
def __init__(self, l, **kwargs):
"""
Create Poisson random variable node
"""
super().__init__(l, **kwargs)
def __str__(self):
"""
Print the distribution using standard parameterization.
"""
l = self.u[0]
return ("%s ~ Categorical(lambda)\n"
" lambda =\n"
"%s\n"
% (self.name, l))
|
