-
Notifications
You must be signed in to change notification settings - Fork 0
/
DF_mle_minimize_mvMU.py
248 lines (203 loc) · 8.68 KB
/
DF_mle_minimize_mvMU.py
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
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue May 1 20:27:24 2018
Author Dimitrios Fafalis
Import this package as: DFmle
Contents of the package:
Log-likelihood functions of mixture distributions:
--------------------------------------------------
1. logLik_1vM
2. logLik_2vM
3. logLik_3vM
4. logLik_2vM1U
5. logLik_1vM1U
@author: df
"""
import numpy as np
from scipy import stats
# ---------------------------------------------------------------------- #
def logLik_1vM( theta, *args ):
"""
The negative log-likelihood function 'l' to be minimized
in order to estimate the von Mises distribution parameters:
kappa1: the concentration of the von Mises distribution
mu1: the location of the von Mises distribution
theta:= [ kappa1, mu1 ]
args[0]:= X_samples: the observations sample
args[1]:= I_samples: the intensities of the observations sample
The function returns the negative log-likelihood function.
This function is to be called with optimize.minimize function of scipy:
results = optimize.minimize( logLik_2vM, in_guess, args=(r_X,Int) ...)
"""
scal_ = 0.5
x_ = np.array(args[0]).T
i_ = np.array(args[1]).T
# the unknown parameters:
kappa = theta[0]
mu = theta[1]
# the 1st von Mises distribution on semi-circle:
fvm_ = stats.vonmises.pdf( x_, kappa, mu, scal_ )
# Calculate the negative log-likelihood as the negative sum of the log
# of the PDF of a von Mises distributions, with location mu and
# concentration kappa:
logLik = -np.sum( np.multiply( np.log( fvm_ ), i_ ) )
return logLik
# ---------------------------------------------------------------------- #
def logLik_2vM( theta, *args ):
"""
The negative log-likelihood function 'l' to be minimized
in order to estimate the mixture parameters:
p1, p2: the weights of the 2 von Mises distributions on semi-circle
kappa1, kappa2: the concentrations of the 2 von Mises distributions
mu1, mu2: the locations of the 2 von Mises distributions
theta:= [ p1, kappa1, mu1, p2, kappa2, mu2 ]
args[0]:= X_samples: the observations sample
args[1]:= I_samples: the intensities of the observations sample
The function returns the negative log-likelihood function.
This function is to be called with optimize.minimize function of scipy:
results = optimize.minimize( logLik_2vM, in_guess, args=(r_X,Int) ...)
"""
# the scale is 1.0 for the circle and 0.5 for the semi-circle:
scal_ = 0.5
x_ = np.array(args[0]).T
i_ = np.array(args[1]).T
# the unknown parameters:
p1_ = theta[0]
kappa1_ = theta[1]
m1_ = theta[2]
p2_ = theta[3]
kappa2_ = theta[4]
m2_ = theta[5]
# the 1st von Mises distribution on semi-circle:
fvm1_ = stats.vonmises.pdf( x_, kappa1_, m1_, scal_ )
# the 2nd von Mises distribution on semi-circle:
fvm2_ = stats.vonmises.pdf( x_, kappa2_, m2_, scal_ )
# mixture distribution:
fm_ = p1_*fvm1_ + p2_*fvm2_
logMix = -np.sum( np.multiply( np.log( fm_ ), i_ ) )
return logMix
# ---------------------------------------------------------------------- #
def logLik_3vM( theta, *args ):
"""
The negative log-likelihood function 'l' to be minimized
in order to estimate the mixture parameters:
p1, p2, p3: the weights of the 3 von Mises distributions on semi-circle
kappa1, kappa2, kappa3: the concentrations of the 3 von Mises
mu1, mu2, mu3: the locations of the 3 von Mises distributions
theta:= [ p1, kappa1, mu1, p2, kappa2, mu2, p3, kappa3, mu3 ]
args[0]:= X_samples: the observations sample
args[1]:= I_samples: the intensities of the observations sample
The function returns the negative log-likelihood function.
This function is to be called with optimize.minimize function of scipy:
results = optimize.minimize( logLik_3vM, in_guess, args=(r_X,Int) ...)
"""
# the scale is 1.0 for the circle and 0.5 for the semi-circle:
scal_ = 0.5
x_ = np.array(args[0]).T
i_ = np.array(args[1]).T
# print(type(x_)), print('x_=', x_.shape)
# print(type(i_)), print('i_=', i_.shape)
# the unknown parameters:
p1_ = theta[0]
kappa1_ = theta[1]
m1_ = theta[2]
p2_ = theta[3]
kappa2_ = theta[4]
m2_ = theta[5]
p3_ = theta[6]
kappa3_ = theta[7]
m3_ = theta[8]
# the 1st von Mises distribution on semi-circle:
fvm1_ = stats.vonmises.pdf( x_, kappa1_, m1_, scal_ )
# logvM1 = -np.sum( stats.vonmises.logpdf( x_, kappa1_, m1_, scal_ ) )
# the 2nd von Mises distribution on semi-circle:
fvm2_ = stats.vonmises.pdf( x_, kappa2_, m2_, scal_ )
# logvM2 = -np.sum( stats.vonmises.logpdf( x_, kappa2_, m2_, scal_ ) )
# the 3rd von Mises distribution on semi-circle:
fvm3_ = stats.vonmises.pdf( x_, kappa3_, m3_, scal_ )
# mixture distribution:
fm_ = p1_*fvm1_ + p2_*fvm2_ + p3_*fvm3_
# logMix = logvM1 + logvM2 + logU - len(x_)*(np.log(p1_*p2_*pu_))
# if the r.s. is the x_ only:
# logMix = -np.sum( np.log( fm_ ) )
# if the r.s. is equally-spaced angles with varying intensity i_:
logMix = -np.sum( np.multiply( np.log( fm_ ), i_ ) )
return logMix
# ---------------------------------------------------------------------- #
def logLik_2vM1U( theta, *args ):
"""
The negative log-likelihood function 'l' to be minimized
in order to estimate the mixture parameters:
p1, p2: the weights of the 2 von Mises distributions on semi-circle
pu: the weight of the uniform distribution on the semi-circle
kappa1, kappa2: the concentrations of the 2 von Mises distributions
mu1, mu2: the locations of the 2 von Mises distributions
theta:= [ p1, kappa1, mu1, p2, kappa2, mu2, pu ]
args[0]:= X_samples: the observations sample
args[1]:= I_samples: the intensities of the observations sample
The function returns a vector F with the derivatives of 'l' wrt the
components of theta.
The function returns the negative log-likelihood function.
This function is to be called with optimize.minimize function of scipy:
results = optimize.minimize( logLik_2vM1U, in_guess, args=(r_X,Int) ..)
"""
scal_ = 0.5
x_ = np.array(args[0]).T
i_ = np.array(args[1]).T
# the unknown parameters:
p1_ = theta[0]
kappa1_ = theta[1]
m1_ = theta[2]
p2_ = theta[3]
kappa2_ = theta[4]
m2_ = theta[5]
pu_ = theta[6]
# the 1st von Mises distribution on semi-circle:
fvm1_ = stats.vonmises.pdf( x_, kappa1_, m1_, scal_ )
# the 2nd von Mises distribution on semi-circle:
fvm2_ = stats.vonmises.pdf( x_, kappa2_, m2_, scal_ )
# the uniform distribution:
xu_1 = min(x_)
xu_2 = (max(x_) - min(x_))
fu_ = stats.uniform.pdf( x_, loc=xu_1, scale=xu_2 )
# mixture distribution:
fm_ = p1_*fvm1_ + p2_*fvm2_ + pu_*fu_
logMix = -np.sum( np.multiply( np.log( fm_ ), i_ ) )
return logMix
# ---------------------------------------------------------------------- #
def logLik_1vM1U( theta, *args ):
"""
The negative log-likelihood function 'l' to be minimized
in order to estimate the mixture parameters:
p1: the weight of the von Mises distribution on semi-circle
pu: the weight of the uniform distribution on the semi-circle
kappa1: the concentration of the von Mises distribution
mu1: the location of the von Mises distributions
theta:= [ p1, kappa1, mu1, pu ]
args[0]:= X_samples: the observations sample (the angles in our case)
args[1]:= I_samples: the intensities of the observations sample
The function returns a vector F with the derivatives of 'l' wrt the
components of theta.
The function returns the negative log-likelihood function.
This function is to be called with optimize.minimize function of scipy:
results = optimize.minimize( logLik_1vM1U, in_guess, args=(r_X,Int) ..)
"""
scal_ = 0.5
x_ = np.array(args[0]).T
i_ = np.array(args[1]).T
# the unknown parameters:
p1_ = theta[0]
kappa1_ = theta[1]
m1_ = theta[2]
pu_ = theta[3]
# the von Mises distribution on semi-circle:
fvm1_ = stats.vonmises.pdf( x_, kappa1_, m1_, scal_ )
# the uniform distribution:
xu_1 = min(x_)
xu_2 = (max(x_) - min(x_))
fu_ = stats.uniform.pdf( x_, loc=xu_1, scale=xu_2 )
# mixture distribution:
fm_ = p1_*fvm1_ + pu_*fu_
logMix = -np.sum( np.multiply( np.log( fm_ ), i_ ) )
return logMix