-
Notifications
You must be signed in to change notification settings - Fork 3
/
likelihood.py
145 lines (127 loc) · 5.04 KB
/
likelihood.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
__author__ = "Johann Cohen Tanugi, Andrea Chiappo"
__email__ = "[email protected]"
#
# This module contains the classes to compute the (log)likelihood
# of the stellar kinematic data given a certain set of model
# parameter values entering Jeans equation formalism
#
# Currently, only one a Gaussian likelihood of the stellar
# line-of-sight velocities is implemented
#
import numpy as np
from sys import float_info
from exceptions import ValueError, Exception
from dispersion import SphericalJeansDispersion
class LogLikelihood(object):
"""
Base class of a (negative) LogLikelihood object
Defines the functions to set a specific model in the Jeans analysis
and to update the model parameter values explored
"""
def __init__(self, data, sigma, numprocs=1):
"""
To build an istance, you must provide
- data : list compraising (in this order) [R, v, dv], where
> R : projected radial distances
> v : line-of-sight stellar velocities
> dv : uncertainty on the line-of-sight velocities
- sigma : istance of SphericalJeansDispersion class
"""
if not isinstance(sigma, SphericalJeansDispersion):
raise Exception("sigma must be an istance of SphericalJeansDispersion")
self.data = data
self.sigma = sigma
self.sigma.data = data
self.free_pars = {'J':18}
self.numprocs = numprocs
self.cache = {}
def set_free(self, parname, **kwargs):
"""
function to let a parameter vary in the LogLikelihood optimisation
input : parname (parameter name)
"""
if parname not in self.sigma.params:
raise ValueError('%s not a parameter of logLike function.\n'%parname\
+'\t The parameters are %s'%self.sigma.params.keys())
if parname not in self.free_pars:
if kwargs == {}:
kwargs['val'] = self.sigma.params[parname]
self.free_pars[parname] = kwargs
def set_fixed(self, parname, **kwargs):
"""
function to fix a parameter to its default value
in the LogLikelihood optimisation
input : parname (parameter name)
"""
if parname not in self.sigma.params:
raise ValueError('%s not a parameter of logLike function.\n'%parname\
+'\t The parameters are %s'%self.sigma.params.keys())
if parname not in self.free_pars:
raise ValueError('%s is already fixed at value %g'%\
(parname, self.sigma.params[parname]))
else:
del self.free_pars[parname]
def set_value(self, parname, parvalue):
"""
function to set a parameter to a given value
input :
- parname : parameter name
- parvalue : parameter value
"""
if parname not in self.sigma.params:
raise ValueError('%s not a parameter of logLike function.\n'%parname\
+'\t The parameters are %s'%self.sigma.params.keys())
if parname in self.free_pars:
self.free_pars[parname] = parvalue
self.sigma.setparams(parname, parvalue)
def _retrieve(self, *freepars):
if freepars in self.cache.keys():
S = self.cache[freepars]
iscached = True
else:
S = 0 #dummy value
iscached = False
return iscached, S
def _store(self, S, *freepars):
self.cache[freepars] = S
def __call__(self, *par_array):
"""
function to evaluate the (negative) LogLikelihood
"""
if np.any(np.isnan(par_array)):
return float_info.max
else:
iscached, Scached = self._retrieve(*par_array)
if iscached:
S = Scached
else:
for i,key in enumerate(self.free_pars.keys()):
self.free_pars[key] = par_array[i]
self.sigma.setparams(key, par_array[i])
S = self.compute()
#self._store(S, *par_array)
return S
class GaussianLikelihood(LogLikelihood):
"""
Class to define a Gaussian loglikelihood of the observed stellar
line-of-sight velocities, given a set of model parameter values
"""
def __init__(self, *args):
"""
To build an instance, you must provide the same quantities as
in the __init__ method of the base class LogLikelihood
"""
super(GaussianLikelihood, self).__init__(*args)
self.R = self.data[0]
self.v = self.data[1]
self.dv2 = self.data[2]**2
self.vsys = self.data[3]
self.Dv2 = (self.v-self.vsys)**2
def compute(self):
"""
function to evaluate the (negative) Gaussian loglikelihood
of the stellar line-of-sight velocities
"""
S = self.dv2 + self.sigma.compute(self.R) #this is an array like R array
res = np.log(S) + self.Dv2/S
return res.sum() / 2.