-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlikelihood.pyx
212 lines (194 loc) · 12.6 KB
/
likelihood.pyx
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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
from libc.math cimport lgamma
from libc.math cimport exp
from libc.math cimport pow
from libc.math cimport log
from libc.math cimport sqrt
from scipy.optimize import minimize_scalar
from libc.math cimport pi
from scipy.integrate import quad
import cython
@cython.cdivision(True)
cdef double expit(double p):
return 1.0 / (1 + exp(-p))
@cython.cdivision(True)
cdef double second_order_derivative_nob(double alpha, double beta, double kappa, double tau, double theta_g,
double sigma_g, double mu_cg, long y_cg):
cdef double cg
if pi * sigma_g * sigma_g == 0:
return float('inf')
if y_cg == 0:
cg = ((0.2e1 * tau * tau * pow(exp(mu_cg * tau + kappa), 0.2e1) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.3e1) - tau * tau * exp(mu_cg * tau + kappa) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.2e1) + 0.2e1 * exp(-exp(beta * mu_cg + alpha)) * tau * tau * pow(exp(-mu_cg * tau - kappa), 0.2e1) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.3e1) - 0.2e1 * beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - exp(-exp(beta * mu_cg + alpha)) * tau * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - beta * beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa)) + beta * beta * pow(exp(beta * mu_cg + alpha), 0.2e1) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) / 0.2e1 - 0.2e1 * (-tau * exp(mu_cg * tau + kappa) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.2e1) + exp(-exp(beta * mu_cg + alpha)) * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * (mu_cg - theta_g) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) * pow(sigma_g, -0.2e1) - (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) * pow(sigma_g, -0.2e1) + 0.2e1 * (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * pow(mu_cg - theta_g, 0.2e1) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) * pow(sigma_g, -0.4e1)) * sqrt(0.2e1) * sqrt(pi * sigma_g * sigma_g) / (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) / exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) - ((-tau * exp(mu_cg * tau + kappa) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.2e1) + exp(-exp(beta * mu_cg + alpha)) * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) / 0.2e1 - (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * (mu_cg - theta_g) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) * pow(sigma_g, -0.2e1)) * sqrt(0.2e1) * sqrt(pi * sigma_g * sigma_g) * (-tau * exp(mu_cg * tau + kappa) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.2e1) + exp(-exp(beta * mu_cg + alpha)) * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * pow(0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa)), -0.2e1) / exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) + 0.2e1 * ((-tau * exp(mu_cg * tau + kappa) * pow(0.1e1 + exp(mu_cg * tau + kappa), -0.2e1) + exp(-exp(beta * mu_cg + alpha)) * tau * exp(-mu_cg * tau - kappa) * pow(0.1e1 + exp(-mu_cg * tau - kappa), -0.2e1) - beta * exp(beta * mu_cg + alpha) * exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) / 0.2e1 - (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) * sqrt(0.2e1) * (mu_cg - theta_g) * exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(pi * sigma_g * sigma_g, -0.1e1 / 0.2e1) * pow(sigma_g, -0.2e1)) * sqrt(0.2e1) * sqrt(pi * sigma_g * sigma_g) * (mu_cg - theta_g) / (0.1e1 / (0.1e1 + exp(mu_cg * tau + kappa)) + exp(-exp(beta * mu_cg + alpha)) / (0.1e1 + exp(-mu_cg * tau - kappa))) / exp(-pow(mu_cg - theta_g, 0.2e1) * pow(sigma_g, -0.2e1)) * pow(sigma_g, -0.2e1)
else:
cg = -(exp((beta * mu_cg - 2 * mu_cg * tau + alpha - 2 * kappa)) * (beta * beta) * sigma_g * sigma_g + 0.2e1 * (beta * beta) * exp((beta * mu_cg - mu_cg * tau + alpha - kappa)) * sigma_g * sigma_g + (tau * tau) * exp((-mu_cg * tau - kappa)) * sigma_g * sigma_g + (beta * beta) * exp((beta * mu_cg + alpha)) * sigma_g * sigma_g + 0.2e1 * exp((-2 * mu_cg * tau - 2 * kappa)) + 0.4e1 * exp((-mu_cg * tau - kappa)) + 0.2e1) * pow(sigma_g, -0.2e1) * pow(0.1e1 + exp((-mu_cg * tau - kappa)), -0.2e1)
return cg
cdef double log_dpois0(double log_mean):
return -exp(log_mean)
cdef double log_dpois(long count, double log_mean):
return count * log_mean - lgamma(count + 1) - exp(log_mean)
cdef double log_expit(double x):
return -log(1.0 + exp(-x))
cdef double log_sum_exp2(double a, double b):
cdef double max_el = max(a, b)
return max_el + log(exp(a - max_el) + exp(b - max_el))
@cython.cdivision(True)
cdef double log_dnorm(double x, double mu, double sigma):
if sigma == 0.0:
if x == mu:
return 0.0
else:
return -float('inf')
else:
return -0.918938533204672669540968854562379419803619384765625 - log(sigma) - (x - mu) * (
x - mu) / sigma / sigma / 2.0
cdef double neg_log_single_complete_likelihood_nob(double mu_cg, double a_c, double b_c, double k_c, double t_c,
double theta_g, double sigma_g, long y_cg):
if y_cg == 0:
return -(log_sum_exp2(log_expit(-(k_c + t_c * mu_cg)),
log_expit(k_c + t_c * mu_cg) + log_dpois0(a_c + b_c * mu_cg)) + log_dnorm(mu_cg, theta_g,
sigma_g))
else:
return -(log_expit(k_c + t_c * mu_cg) + log_dpois(y_cg, a_c + b_c * mu_cg) + log_dnorm(mu_cg, theta_g, sigma_g))
cdef double single_complete_likelihood_nob(double mu_cg, double a_c, double b_c, double k_c, double t_c, double theta_g,
double sigma_g, long y_cg, double scale_factor_cg):
return exp(
-neg_log_single_complete_likelihood_nob(mu_cg, a_c, b_c, k_c, t_c, theta_g, sigma_g, y_cg) + scale_factor_cg)
cdef double neg_log_single_marginal_likelihood_nob(double a_c, double b_c, double k_c, double t_c, double theta_g,
double sigma_g, long y_cg):
# first get the min of the neg log-likelihood
# use brent method
min_neg_log = minimize_scalar(neg_log_single_complete_likelihood_nob,
args=(a_c, b_c, k_c, t_c, theta_g, sigma_g, y_cg), method='brent')
cdef double min_val
cdef double hessian
cdef double lower_b
cdef double upper_b
cdef double arg_min
if min_neg_log.success:
arg_min = min_neg_log.x
min_val = min_neg_log.fun
hessian = second_order_derivative_nob(a_c, b_c, k_c, t_c, theta_g, sigma_g, arg_min, y_cg)
lower_b = arg_min - 20 / sqrt(abs(hessian))
upper_b = arg_min + 20 / sqrt(abs(hessian))
integral = quad(single_complete_likelihood_nob, lower_b, upper_b,
args=(a_c, b_c, k_c, t_c, theta_g, sigma_g, y_cg, min_val))
return -(log(integral[0]) - min_val)
else:
return float('nan')
def neg_log_sum_marginal_likelihood_nob(real_params_g, abkt, y_g):
cdef double sum_marginal_likelihood = 0
cdef double a_c
cdef double b_c
cdef double k_c
cdef double t_c
cdef double theta_g = real_params_g[0]
cdef double sigma_g = exp(real_params_g[1])
for i in range(len(y_g)):
a_c = abkt[i, 0]
b_c = abkt[i, 1]
k_c = abkt[i, 2]
t_c = abkt[i, 3]
sum_marginal_likelihood += neg_log_single_marginal_likelihood_nob(a_c, b_c, k_c, t_c, theta_g, sigma_g, y_g[i])
return sum_marginal_likelihood
cdef extern from "math.h":
bint isnan(double x)
def neg_log_sum_marginal_likelihood(real_params_g, abkt, y_g):
cdef double sum_marginal_likelihood = 0
cdef double a_c
cdef double b_c
cdef double k_c
cdef double t_c
cdef double theta_g = real_params_g[0]
cdef double sigma_g = exp(real_params_g[1])
cdef double p_g = expit(real_params_g[2])
cdef double t2 = 0
for i in range(len(y_g)):
a_c = abkt[i, 0]
b_c = abkt[i, 1]
k_c = abkt[i, 2]
t_c = abkt[i, 3]
t2 = neg_log_single_marginal_likelihood_nob(a_c,b_c,k_c,t_c,theta_g,sigma_g,y_g[i])
if isnan(t2):
return float('nan')
if y_g[i]>0:
sum_marginal_likelihood += (t2-log(p_g))
elif y_g[i]==0:
sum_marginal_likelihood += (-log_sum_exp2(log(1-p_g),-t2+log(p_g)))
return sum_marginal_likelihood
def neg_log_sum_marginal_likelihood_free_p(real_params_g, abkt, y_g, x_g):
cdef double sum_marginal_likelihood = 0
cdef double a_c
cdef double b_c
cdef double k_c
cdef double t_c
cdef double theta_g
cdef double sigma_g
cdef double p_g
cdef double t2 = 0
for i in range(len(y_g)):
a_c = abkt[i, 0]
b_c = abkt[i, 1]
k_c = abkt[i, 2]
t_c = abkt[i, 3]
theta_g = real_params_g[0]
sigma_g = exp(real_params_g[1])
p_g = expit(real_params_g[2]) * (1 - x_g[i]) + expit(real_params_g[3]) * x_g[i]
t2 = neg_log_single_marginal_likelihood_nob(a_c,b_c,k_c,t_c,theta_g,sigma_g,y_g[i])
if isnan(t2):
return float('nan')
if y_g[i]>0:
sum_marginal_likelihood += (t2-log(p_g))
elif y_g[i]==0:
sum_marginal_likelihood += (-log_sum_exp2(log(1-p_g),-t2+log(p_g)))
return sum_marginal_likelihood
def neg_log_sum_marginal_likelihood_free_theta(real_params_g, abkt, y_g, x_g):
cdef double sum_marginal_likelihood = 0
cdef double a_c
cdef double b_c
cdef double k_c
cdef double t_c
cdef double theta_g
cdef double sigma_g
cdef double p_g
cdef double t2 = 0
for i in range(len(y_g)):
a_c = abkt[i, 0]
b_c = abkt[i, 1]
k_c = abkt[i, 2]
t_c = abkt[i, 3]
theta_g = real_params_g[0] * (1 - x_g[i]) + real_params_g[1] * x_g[i]
sigma_g = exp(real_params_g[2])
p_g = expit(real_params_g[3])
t2 = neg_log_single_marginal_likelihood_nob(a_c,b_c,k_c,t_c,theta_g,sigma_g,y_g[i])
if isnan(t2):
return float('nan')
if y_g[i]>0:
sum_marginal_likelihood += (t2-log(p_g))
elif y_g[i]==0:
sum_marginal_likelihood += (-log_sum_exp2(log(1-p_g),-t2+log(p_g)))
return sum_marginal_likelihood
def neg_log_sum_marginal_likelihood_free_both(real_params_g, abkt, y_g, x_g):
cdef double sum_marginal_likelihood = 0
cdef double a_c
cdef double b_c
cdef double k_c
cdef double t_c
cdef double theta_g
cdef double sigma_g
cdef double p_g
cdef double t2 = 0
for i in range(len(y_g)):
a_c = abkt[i, 0]
b_c = abkt[i, 1]
k_c = abkt[i, 2]
t_c = abkt[i, 3]
theta_g = real_params_g[0] * (1 - x_g[i]) + real_params_g[1] * x_g[i]
sigma_g = exp(real_params_g[2])
p_g = expit(real_params_g[3]) * (1 - x_g[i]) + expit(real_params_g[4]) * x_g[i]
t2 = neg_log_single_marginal_likelihood_nob(a_c,b_c,k_c,t_c,theta_g,sigma_g,y_g[i])
if isnan(t2):
return float('nan')
if y_g[i]>0:
sum_marginal_likelihood += (t2-log(p_g))
elif y_g[i]==0:
sum_marginal_likelihood += (-log_sum_exp2(log(1-p_g),-t2+log(p_g)))
return sum_marginal_likelihood