1
$\begingroup$

I am thinking I can't be the only one encountering this. I am trying to do maximum likelihood estimation on a probit model, i.e. trying to find the most optimal fit for three parameters in my case through the following equation:

$$l\left(\beta; y, X\right) = \sum_{i=1}^N \left[y_i \ln\left(F\left(x_i \beta\right)\right) + \left(1 - y_i\right) \ln \left(1-F\left(x_i \beta\right)\right)\right]$$

$y_i$ can take the value 0 or 1 depending on whether there has been an event or not, and then some form of probability inside the $\ln$ function. So in my case I am just computing a large grid of parameters (which is used as input into the $F$ function), and then doing this computation over $i$ number of cases.

However, my problem is that some parameter sets (in the grid) will return either 1 or 0 as probability. That is unavoidable. And in those $i$-th cases (depending on whether the corresponding $y_i$ is 1 or 0) the result will be $\log(0)$, which is not something you would want I imagine.

How do one deal with this problem, or should it be dealt with at all ? My first assumption was to first calculate the $F$ function for all cases $i$ for all parameters in the grid. Then figure out which specific elements/indexes have the value 1 or 0, and subtract all combined indexes from all $i$ cases in the grid, so I wouldn't have 1 or 0 anywhere, and all cases has had those particular indexes removed all over the board. But that seems like some kind of weird solution to me, so that's what I am here for.

$\endgroup$
  • 1
    $\begingroup$ You can impose the rule that $0\times\log(0)=0$ and take $1\times\log(0)$ as a very large negative value. $\endgroup$ – Xi'an Sep 4 '18 at 11:40
  • 1
    $\begingroup$ Why are probabilities of 0 and 1 unavoidable in a probit model? Or are you talking about some problem with finite precision? $\endgroup$ – The Laconic Sep 4 '18 at 12:06
  • $\begingroup$ @TheLaconic I don't know if that is general, but in my case it is. Or at least the way I am using MLE on my data. I have like a 100x100x100 grid with 100 of each of the values I am trying to estimate (in a range I know is more probable than others). But obviously at some of the upper/lower bounds of the parameters they are so far from the true value that they in turn yield the probability 0 or 1 (that is the value $F$ can take in my case. Some value between 0 and 1). And if for that case $y_i$ is 0, and its probability is 1, then I will end up with $ln(0)$. $\endgroup$ – Denver Dang Sep 4 '18 at 13:09
  • $\begingroup$ @Xi'an It just seems weird that while large MLE values indicate a good agreement with parameters, i.e. if you have a lot of $i$-th cases where the MLE is large, that particular set of parameters will most likely be the "chosen set" since that would maximize the log-loglikelihood the most when summed together for all $i$. And now I impose either a $0 \times \log(0) = 0$ for a set of parameters that should not be very likely, i.e. a high value since most MLE values in this case are negative, or another one, equally unlikely, as very unlikely (i.e. very large negative number). It seems wrong. $\endgroup$ – Denver Dang Sep 4 '18 at 13:32
  • 1
    $\begingroup$ F is just the normal CDF. It doesn’t return 0 or 1 for a finite argument. So what’s going on here? Is it returning 0 or 1 to within numerical precision, for certain parameter values? $\endgroup$ – The Laconic Sep 4 '18 at 14:03
2
$\begingroup$

I agree with TheLaconic that this is a numerical precision issue; the easiest work-around would be to use the scipy.stats.norm.logcdf method, instead of computing the CDF and then taking the log, which sometimes creates exactly this kind of numerical precision issue.

For the component that depends on $1 - F(\cdot)$, consider that the so-called survival function $s(a)$ is defined as $s(a) = 1 - F(a)$, so we can use scipy.stats.norm.logsf to take care of that without even doing algebra.

$\endgroup$
  • $\begingroup$ Hmmm, I haven't noticed that function, and it do seem to work. However, I do run into some trouble when running into the second term of the function $\left(1 - y_i\right) \ln \left(1-F\left(x_i \beta\right)\right)$. For the first term ($y_i \ln\left(F\left(x_i \beta\right)\right)$) I could easily just use the logcdf, but in the case that $y_i = 0$ I would in principle be doing $\ln(1-\ln(F))$ which wouldn't be the same as the original. $\endgroup$ – Denver Dang Sep 4 '18 at 17:13
  • $\begingroup$ @DenverDang I've updated my answer. $\endgroup$ – Sycorax Sep 4 '18 at 17:32
  • $\begingroup$ Ah, that is brilliant. One question though, should that only be used when I know that it will be in the $1-F(x)$ term, or should it just be used in all cases ? For example, I am using if statements to say if y_i = 1, and then I would usually just do log(scipy.stats.norm(0, 1).cdf(x)) for that and log(1 - scipy.stats.norm(0, 1).cdf(x)) for $y_i = 0$. So should that instead be: scipy.stats.norm(0, 1).logcdf(x) for $y_i = 1$ and scipy.stats.norm(0, 1).logsf(x)* for $y_i = 0$ ? $\endgroup$ – Denver Dang Sep 4 '18 at 17:54
  • 1
    $\begingroup$ The way I would code it would be y * scipy.stats.norm(0, 1).logcdf(x) + (1 - y) * scipy.stats.norm(0, 1).logsf(x) $\endgroup$ – Sycorax Sep 4 '18 at 17:57
  • $\begingroup$ Thank you so much, I didn't know this existed :) $\endgroup$ – Denver Dang Sep 4 '18 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.