Theano warnings when fitting a simple psychometric function (cumulative normal distribution)

Questions
1

I am getting the following theano warnings when fitting a psychometric function using a cumulative normal distribution to response data.

ERROR (theano.graph.opt): Optimization failure due to: local_grad_log_erfc_neg
ERROR (theano.graph.opt): node: Elemwise{true_div,no_inplace}(Elemwise{mul,no_inplace}.0, Elemwise{erfc,no_inplace}.0)
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "C:\ProgramData\Anaconda3\lib\site-packages\theano\graph\opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "C:\ProgramData\Anaconda3\lib\site-packages\theano\graph\opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "C:\ProgramData\Anaconda3\lib\site-packages\theano\tensor\opt.py", line 7307, in local_grad_log_erfc_neg
    if not exp_in.owner.inputs[0].owner:
AttributeError: 'NoneType' object has no attribute 'owner'

The tt.erf function seems to be causing the problem here, as the warnings are not present when using a logistic function instead (that does not require this function). This code snippet should reproduce the behaviour:

import numpy as np
import pymc3 as pm
import arviz as az
from theano import tensor as tt


def comulative_normal(x, alpha, beta, s=np.sqrt(2)):
    # Cumulative distribution function for the standard normal distribution
    return 0.5 * (1 + tt.erf((x-alpha)/(beta*s)))

xij = np.array([-40.5, -40. , -39.5, -39. , -38.5, -38. , -37.5, -37. , -36.5,
       -36. , -35.5, -35. , -34.5, -34. , -33.5, -33. , -32.5, -32. ,
       -31.5, -31. , -30.5, -30. , -29.5, -29. , -28.5, -28. , -27.5,
       -27. , -26.5, -26. , -25.5, -25. , -24.5, -24. , -23.5, -23. ,
       -22.5, -22. , -21.5, -21. , -20.5, -20. , -19.5, -19. , -18.5,
       -18. , -17.5, -17. , -16.5, -16. , -15.5, -15. , -14.5, -14. ,
       -13.5, -13. , -12.5, -12. , -11.5, -11. , -10.5, -10. ,  -9.5,
        -9. ,  -8.5,  -8. ,  -7.5,  -7. ,  -6.5,  -6. ,  -5.5,  -5. ,
        -4.5,  -4. ,  -3.5,  -3. ,  -2.5,  -2. ,  -1.5,  -1. ,  -0.5,
         0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. ,
         4.5,   5. ,   5.5,   6. ,   6.5,   7. ,   7.5,   8. ,   8.5,
         9. ,   9.5,  10. ,  10.5,  11. ,  11.5,  12. ,  12.5,  13. ,
        13.5,  14. ,  14.5,  15. ,  15.5,  16. ,  16.5,  17. ,  17.5,
        18. ,  18.5,  19. ,  19.5,  20. ,  20.5,  21. ,  21.5,  22. ,
        22.5,  23. ,  23.5,  24. ,  24.5,  25. ,  25.5,  26. ,  26.5,
        27. ,  27.5,  28. ,  28.5,  29. ,  29.5,  30. ,  30.5,  31. ,
        31.5,  32. ,  32.5,  33. ,  33.5,  34. ,  34.5,  35. ,  35.5,
        36. ,  36.5,  37. ,  37.5,  38. ,  38.5,  39. ,  39.5,  40. ,
        40.5])

nij = np.array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 1., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 3., 0., 5., 0., 8., 0.,
       5., 0., 4., 0., 2., 0., 1., 0., 1., 0., 2., 0., 3., 0., 4., 0., 3.,
       0., 2., 0., 2., 0., 3., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 4., 0., 3.,
       0., 9., 0., 5., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 1.])

rij = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2., 0.,
       2., 0., 1., 0., 0., 0., 1., 0., 1., 0., 1., 0., 0., 0., 2., 0., 1.,
       0., 1., 0., 1., 0., 2., 0., 0., 0., 0., 0., 2., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 4., 0., 3.,
       0., 9., 0., 5., 0., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 1.])

with pm.Model():

    alpha = pm.Uniform('alpha', lower=-40.5, upper=40.5)
    
    beta = pm.Uniform("beta", lower=10e-5, upper=40)

    thetaij = pm.Deterministic('thetaij', comulative_normal(xij, alpha, beta))

    rij_ = pm.Binomial('rij', p=thetaij, n=nij, observed=rij)
    trace = pm.sample(chains=4, cores=4, tune=2000)
```

Besides that, the model is still sampling "correctly" and returns credible results. I just want to make sure I am not introducing a weird behaviour with this custom function.

I am using pmc3 3.11.0 and theano-PYMC 1.1.0
2

@LegrandNico Thanks for your question. I was able to replicate the issue on my side. Funnily enough if you change your function to be:

def cumulative_normal(x, alpha, beta, s=np.sqrt(2)):
    # Cumulative distribution function for the standard normal distribution
    return 0.5 + 0.5 * tt.erf((x-alpha)/(beta*s))

It seems to work fine. I have opened an issue on the GH: ERROR (theano.graph.opt): Optimization failure due to: local_grad_log_erfc_neg · Issue #291 · pymc-devs/aesara · GitHub

1 Like
3

Thank you @ricardoV94 for your response.

I can confirm that changing the function works fine, and the results are comparable to what is obtained using the first definition.