Coverage for brainbox/behavior/pyschofit.py: 0%
92 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-10-11 11:13 +0100
1"""
2(DEPRECATED) The psychofit toolbox contains tools to fit two-alternative psychometric
3data. The fitting is done using maximal likelihood estimation: one
4assumes that the responses of the subject are given by a binomial
5distribution whose mean is given by the psychometric function.
6The data can be expressed in fraction correct (from .5 to 1) or in
7fraction of one specific choice (from 0 to 1). To fit them you can use
8these functions:
9 - weibull50: Weibull function from 0.5 to 1, with lapse rate
10 - weibull: Weibull function from 0 to 1, with lapse rate
11 - erf_psycho: erf function from 0 to 1, with lapse rate
12 - erf_psycho_2gammas: erf function from 0 to 1, with two lapse rates
13Functions in the toolbox are:
14 - mle_fit_psycho: Maximumum likelihood fit of psychometric function
15 - neg_likelihood: Negative likelihood of a psychometric function
16For more info, see:
17 - Examples: Examples of use of psychofit toolbox
18Matteo Carandini, 2000-2015
20NB: USE THE PSYCHOFIT PIP PACKAGE INSTEAD.
21"""
23import functools
24import warnings
25import traceback
26import logging
28import numpy as np
29import scipy.optimize
30from scipy.special import erf
33for line in traceback.format_stack():
34 print(line.strip())
36msg = 'brainbox.behavior.pyschofit has been deprecated. Install psychofit via pip. See stack above'
37warnings.warn(msg, DeprecationWarning)
38logging.getLogger(__name__).warning(msg)
41def mle_fit_psycho(data, P_model='weibull', parstart=None, parmin=None, parmax=None, nfits=5):
42 """
43 Maximumum likelihood fit of psychometric function.
44 Args:
45 data: 3 x n matrix where first row corresponds to stim levels,
46 the second to number of trials for each stim level (int),
47 the third to proportion correct / proportion rightward (float between 0 and 1)
48 P_model: The psychometric function. Possibilities include 'weibull'
49 (DEFAULT), 'weibull50', 'erf_psycho' and 'erf_psycho_2gammas'
50 parstart: Non-zero starting parameters, used to try to avoid local minima.
51 The parameters are [threshold, slope, gamma], or if using the
52 'erf_psycho_2gammas' model append a second gamma value.
53 Recommended to use a value > 1. If None, some reasonable defaults are used.
54 parmin: Minimum parameter values. If None, some reasonable defaults are used
55 parmax: Maximum parameter values. If None, some reasonable defaults are used
56 nfits: The number of fits
57 Returns:
58 pars: The parameters from the best of the fits
59 L: The likelihood of the best fit
60 Raises:
61 TypeError: data must be a list or numpy array
62 ValueError: data must be m by 3 matrix
63 Examples:
64 Below we fit a Weibull function to some data:
65 >>> import numpy as np
66 >>> import matplotlib.pyplot as plt
67 >>> cc = np.array([-8., -6., -4., -2., 0., 2., 4., 6., 8.]) # contrasts
68 >>> nn = np.full((9,), 10) # number of trials at each contrast
69 >>> pp = np.array([5., 8., 20., 41., 54., 59., 79., 92., 96])/100 # proportion "rightward"
70 >>> pars, L = mle_fit_psycho(np.vstack((cc, nn, pp)), 'erf_psycho')
71 >>> plt.plot(cc, pp, 'bo', mfc='b')
72 >>> plt.plot(np.arange(-8, 8, 0.1), erf_psycho(pars, np.arange(-8, 8, 0.1)), '-b')
73 Information:
74 1999-11 FH wrote it
75 2000-01 MC cleaned it up
76 2000-04 MC took care of the 50% case
77 2009-12 MC replaced fmins with fminsearch
78 2010-02 MC, AZ added nfits
79 2013-02 MC+MD fixed bug with dealing with NaNs
80 2018-08 MW ported to Python
81 """
82 # Input validation
83 if isinstance(data, (list, tuple)):
84 data = np.array(data)
85 elif not isinstance(data, np.ndarray):
86 raise TypeError('data must be a list or numpy array')
88 if data.shape[0] != 3:
89 raise ValueError('data must be m by 3 matrix')
91 rep = lambda x: (x, x) if P_model.endswith('2gammas') else (x,) # noqa
92 if parstart is None:
93 parstart = np.array([np.mean(data[0, :]), 3., *rep(.05)])
94 if parmin is None:
95 parmin = np.array([np.min(data[0, :]), 0., *rep(0.)])
96 if parmax is None:
97 parmax = np.array([np.max(data[0, :]), 10., *rep(.4)])
99 # find the good values in pp (conditions that were effectively run)
100 ii = np.isfinite(data[2, :])
102 likelihoods = np.zeros(nfits,)
103 pars = np.empty((nfits, parstart.size))
105 f = functools.partial(neg_likelihood, data=data[:, ii],
106 P_model=P_model, parmin=parmin, parmax=parmax)
107 for ifit in range(nfits):
108 pars[ifit, :] = scipy.optimize.fmin(f, parstart, disp=False)
109 parstart = parmin + np.random.rand(parmin.size) * (parmax - parmin)
110 likelihoods[ifit] = -neg_likelihood(pars[ifit, :], data[:, ii], P_model, parmin, parmax)
112 # the values to be output
113 L = likelihoods.max()
114 iBestFit = likelihoods.argmax()
115 return pars[iBestFit, :], L
118def neg_likelihood(pars, data, P_model='weibull', parmin=None, parmax=None):
119 """
120 Negative likelihood of a psychometric function.
121 Args:
122 pars: Model parameters [threshold, slope, gamma], or if
123 using the 'erf_psycho_2gammas' model append a second gamma value.
124 data: 3 x n matrix where first row corresponds to stim levels,
125 the second to number of trials for each stim level (int),
126 the third to proportion correct / proportion rightward (float between 0 and 1)
127 P_model: The psychometric function. Possibilities include 'weibull'
128 (DEFAULT), 'weibull50', 'erf_psycho' and 'erf_psycho_2gammas'
129 parmin: Minimum bound for parameters. If None, some reasonable defaults are used
130 parmax: Maximum bound for parameters. If None, some reasonable defaults are used
131 Returns:
132 ll: The likelihood of the parameters. The equation is:
133 - sum(nn.*(pp.*log10(P_model)+(1-pp).*log10(1-P_model)))
134 See the the appendix of Watson, A.B. (1979). Probability
135 summation over time. Vision Res 19, 515-522.
136 Raises:
137 ValueError: invalid model, options are "weibull",
138 "weibull50", "erf_psycho" and "erf_psycho_2gammas"
139 TypeError: data must be a list or numpy array
140 ValueError data must be m by 3 matrix
141 Information:
142 1999-11 FH wrote it
143 2000-01 MC cleaned it up
144 2000-07 MC made it indep of Weibull and added parmin and parmax
145 2018-08 MW ported to Python
146 """
147 # Validate input
148 if isinstance(data, (list, tuple)):
149 data = np.array(data)
150 elif not isinstance(data, np.ndarray):
151 raise TypeError('data must be a list or numpy array')
153 if parmin is None:
154 parmin = np.array([.005, 0., 0.])
155 if parmax is None:
156 parmax = np.array([.5, 10., .25])
158 if data.shape[0] == 3:
159 xx = data[0, :]
160 nn = data[1, :]
161 pp = data[2, :]
162 else:
163 raise ValueError('data must be m by 3 matrix')
165 # here is where you effectively put the constraints.
166 if (any(pars < parmin)) or (any(pars > parmax)):
167 ll = 10000000
168 return ll
170 dispatcher = {
171 'weibull': weibull,
172 'weibull50': weibull50,
173 'erf_psycho': erf_psycho,
174 'erf_psycho_2gammas': erf_psycho_2gammas
175 }
176 try:
177 probs = dispatcher[P_model](pars, xx)
178 except KeyError:
179 raise ValueError('invalid model, options are "weibull", ' +
180 '"weibull50", "erf_psycho" and "erf_psycho_2gammas"')
182 assert (max(probs) <= 1) or (min(probs) >= 0), 'At least one of the probabilities is not ' \
183 'between 0 and 1'
185 probs[probs == 0] = np.finfo(float).eps
186 probs[probs == 1] = 1 - np.finfo(float).eps
188 ll = - sum(nn * (pp * np.log(probs) + (1 - pp) * np.log(1 - probs)))
189 return ll
192def weibull(pars, xx):
193 """
194 Weibull function from 0 to 1, with lapse rate.
195 Args:
196 pars: Model parameters [alpha, beta, gamma].
197 xx: vector of stim levels.
198 Returns:
199 A vector of length xx
200 Raises:
201 ValueError: pars must be of length 3
202 TypeError: pars must be list-like or numpy array
203 Information:
204 1999-11 FH wrote it
205 2000-01 MC cleaned it up
206 2018-08 MW ported to Python
207 """
208 # Validate input
209 if not isinstance(pars, (list, tuple, np.ndarray)):
210 raise TypeError('pars must be list-like or numpy array')
212 if len(pars) != 3:
213 raise ValueError('pars must be of length 3')
215 alpha, beta, gamma = pars
216 return (1 - gamma) - (1 - 2 * gamma) * np.exp(-((xx / alpha) ** beta))
219def weibull50(pars, xx):
220 """
221 Weibull function from 0.5 to 1, with lapse rate.
222 Args:
223 pars: Model parameters [alpha, beta, gamma].
224 xx: vector of stim levels.
225 Returns:
226 A vector of length xx
227 Raises:
228 ValueError: pars must be of length 3
229 TypeError: pars must be list-like or numpy array
230 Information:
231 2000-04 MC wrote it
232 2018-08 MW ported to Python
233 """
234 # Validate input
235 if not isinstance(pars, (list, tuple, np.ndarray)):
236 raise TypeError('pars must be list-like or numpy array')
238 if len(pars) != 3:
239 raise ValueError('pars must be of length 3')
241 alpha, beta, gamma = pars
242 return (1 - gamma) - (.5 - gamma) * np.exp(-((xx / alpha) ** beta))
245def erf_psycho(pars, xx):
246 """
247 erf function from 0 to 1, with lapse rate.
248 Args:
249 pars: Model parameters [bias, slope, lapse].
250 xx: vector of stim levels.
251 Returns:
252 ff: A vector of length xx
253 Examples:
254 >>> import numpy as np
255 >>> import matplotlib.pyplot as plt
256 >>> xx = np.arange(-50, 50)
257 >>> ff = erf_psycho(np.array([-10., 10., 0.1]), xx)
258 >>> plt.plot(xx, ff)
259 Raises:
260 ValueError: pars must be of length 3
261 TypeError: pars must be a list or numpy array
262 Information:
263 2000 MC wrote it
264 2018-08 MW ported to Python
265 """
266 # Validate input
267 if not isinstance(pars, (list, tuple, np.ndarray)):
268 raise TypeError('pars must be list-like or numpy array')
270 if len(pars) != 3:
271 raise ValueError('pars must be of length 4')
273 (bias, slope, gamma) = pars
274 return gamma + (1 - 2 * gamma) * (erf((xx - bias) / slope) + 1) / 2
277def erf_psycho_2gammas(pars, xx):
278 """
279 erf function from 0 to 1, with two lapse rates.
280 Args:
281 pars: Model parameters [bias, slope, gamma].
282 xx: vector of stim levels (%)
283 Returns:
284 ff: A vector of length xx
285 Examples:
286 >>> import numpy as np
287 >>> import matplotlib.pyplot as plt
288 >>> xx = np.arange(-50, 50)
289 >>> ff = erf_psycho_2gammas(np.array([-10., 10., 0.2, 0.]), xx)
290 >>> plt.plot(xx, ff)
291 Raises:
292 ValueError: pars must be of length 4
293 TypeError: pars must be list-like or numpy array
294 Information:
295 2000 MC wrote it
296 2018-08 MW ported to Python
297 """
298 # Validate input
299 if not isinstance(pars, (list, tuple, np.ndarray)):
300 raise TypeError('pars must be a list-like or numpy array')
302 if len(pars) != 4:
303 raise ValueError('pars must be of length 4')
305 (bias, slope, gamma1, gamma2) = pars
306 return gamma1 + (1 - gamma1 - gamma2) * (erf((xx - bias) / slope) + 1) / 2