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
|
from cython cimport floating
cimport numpy as cnp
import numpy as np
cnp.import_array()
from libc.math cimport exp, fabs, log
def mean_change(const floating[:] arr_1, const floating[:] arr_2):
"""Calculate the mean difference between two arrays.
Equivalent to np.abs(arr_1 - arr2).mean().
"""
cdef cnp.float64_t total, diff
cdef cnp.npy_intp i, size
size = arr_1.shape[0]
total = 0.0
for i in range(size):
diff = fabs(arr_1[i] - arr_2[i])
total += diff
return total / size
def _dirichlet_expectation_1d(
floating[:] doc_topic,
floating doc_topic_prior,
floating[:] out
):
"""Dirichlet expectation for a single sample:
exp(E[log(theta)]) for theta ~ Dir(doc_topic)
after adding doc_topic_prior to doc_topic, in-place.
Equivalent to
doc_topic += doc_topic_prior
out[:] = np.exp(psi(doc_topic) - psi(np.sum(doc_topic)))
"""
cdef floating dt, psi_total, total
cdef cnp.npy_intp i, size
size = doc_topic.shape[0]
total = 0.0
for i in range(size):
dt = doc_topic[i] + doc_topic_prior
doc_topic[i] = dt
total += dt
psi_total = psi(total)
for i in range(size):
out[i] = exp(psi(doc_topic[i]) - psi_total)
def _dirichlet_expectation_2d(const floating[:, :] arr):
"""Dirichlet expectation for multiple samples:
E[log(theta)] for theta ~ Dir(arr).
Equivalent to psi(arr) - psi(np.sum(arr, axis=1))[:, np.newaxis].
Note that unlike _dirichlet_expectation_1d, this function doesn't compute
the exp and doesn't add in the prior.
"""
cdef floating row_total, psi_row_total
cdef floating[:, :] d_exp
cdef cnp.npy_intp i, j, n_rows, n_cols
n_rows = arr.shape[0]
n_cols = arr.shape[1]
d_exp = np.empty_like(arr)
for i in range(n_rows):
row_total = 0
for j in range(n_cols):
row_total += arr[i, j]
psi_row_total = psi(row_total)
for j in range(n_cols):
d_exp[i, j] = psi(arr[i, j]) - psi_row_total
return d_exp.base
# Psi function for positive arguments. Optimized for speed, not accuracy.
#
# After: J. Bernardo (1976). Algorithm AS 103: Psi (Digamma) Function.
# https://www.uv.es/~bernardo/1976AppStatist.pdf
cdef floating psi(floating x) noexcept nogil:
cdef double EULER = 0.577215664901532860606512090082402431
if x <= 1e-6:
# psi(x) = -EULER - 1/x + O(x)
return -EULER - 1. / x
cdef floating r, result = 0
# psi(x + 1) = psi(x) + 1/x
while x < 6:
result -= 1. / x
x += 1
# psi(x) = log(x) - 1/(2x) - 1/(12x**2) + 1/(120x**4) - 1/(252x**6)
# + O(1/x**8)
r = 1. / x
result += log(x) - .5 * r
r = r * r
result -= r * ((1./12.) - r * ((1./120.) - r * (1./252.)))
return result
|
