# Sp.stats and PyMC3 logps different

**URL:** <https://discourse.pymc.io/t/sp-stats-and-pymc3-logps-different/8421>\
**Category:** Questions\
**Created:** [December 9, 2021, 9:06am UTC](https://discourse.pymc.io/t/sp-stats-and-pymc3-logps-different/8421 "2021-12-09T09:06:14Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![juststarted](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/juststarted/32/1997_2.png) [@juststarted](https://discourse.pymc.io/u/juststarted)\
**Post date:** [December 9, 2021, 9:06am UTC](https://discourse.pymc.io/t/sp-stats-and-pymc3-logps-different/8421/1 "2021-12-09T09:06:15Z")

</div>

Hi everyone,

I am fitting a geometric distribution to the following data:

```auto
[40000, 600, 1500, 30000, 12000, 25000, 65000, 1500, 40000, 10000000, 25000, 2000, 2000, 500, 800, 1500, 30000, 850, 25000, 1000, 15000, 40000, 9000, 3000, 12000, 1000, 1000, 1500, 10000, 25000, 7000, 35000, 30000, 25000, 750, 20000, 7000, 1500, 50000, 1250, 17000, 300, 15000, 1000, 50000, 3000, 35000, 8000, 41000]

```

Initially I just wanted to find the MLE, so I used sp.stats and sp.optimize:

```python
f_geom = lambda p: -np.sum(stats.geom.logpmf(k=data, p=p))
mle_geom_parameters = optimize.minimize_scalar(
    f_geom,
    bounds=(0,1),
    method='bounded'
)

```

Which results in:

```auto
     fun: 653.6801753844973
 message: 'Solution found.'
    nfev: 25
  status: 0
 success: True
       x: 5.9608609865491405e-06

```

So far so good.

Then, I wanted to do model comparison with other distributions on the same data, so I thought I would find the full posterior with PyMC3. Here’s the code:

```auto
with pm.Model() as model:
    p = pm.Uniform('p', lower=0, upper=1)
    observed = pm.Geometric('obs', p=p, observed=data)
    logpt = pm.Deterministic('logp', model.logpt)
    trace_geom = pm.sample(
        cores=1,
        tune=2000,
# return_inferencedata=True
    )

```

Just to check, I compared the predicted loglikelihoods on a point:

```python
x = 0.1
print(model.logp({
    'p_interval__': np.log(x / (1-x))
}))
print(np.sum(stats.geom.logpmf(k=data, p=x)))

```

Which prints out:

```python
-1134795.4133526431
-1134793.005407033

```

Is this just numerical error? There is a pattern in the difference between the two, namely the closer to the boundaries of the p parameter, the greater the difference:

```python
xs = np.linspace(0, 1, 1000)

logliks_geom_pymc3 = np.array([
    model.logp({
        'p': x,
        'p_interval__': np.log(x / (1-x))
    })
    for x in xs
])

logliks_geom_np = np.array([
    np.sum(stats.geom.logpmf(k=data, p=x))
    for x in xs
])
plt.plot(xs, np.abs(logliks_geom_pymc3-logliks_geom_np), color='blue')

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/3/39f467af0f140fe50aed4ba61b5f5bf2dd2e682b.png)

I was wondering what might be causing this difference? Thank you!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 9, 2021, 10:11am UTC](https://discourse.pymc.io/t/sp-stats-and-pymc3-logps-different/8421/2 "2021-12-09T10:11:41Z")

</div>

The difference probably comes from the Jacobian term to account for the parameter transformation of `p`

You can try to use `model.logp_nojac` which should be similar to the scipy results

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 9, 2021, 10:16am UTC](https://discourse.pymc.io/t/sp-stats-and-pymc3-logps-different/8421/3 "2021-12-09T10:16:01Z")

</div>

Or you can disable the transform of `p` by using

```python
p = pm.Uniform('p', lower=0, upper=1, transform=None)

```