Best practices for MAP estimation

Questions
1

I am trying to estimate the MAP for a logistic regression model. For simplicity, I’ve reduced the model down to only the intercept term:

import numpy as np
import pymc3 as pm


def sigmoid(x: np.array) -> np.array:
    """Return the sigmoid of X."""
    sigm = 1.0 / (1.0 + np.exp(-x))
    return sigm


data = np.array([0.05, 0.1, 0.1, 0.15, 0.2, 0.1])

with pm.Model():
    const = pm.Normal("const", mu=0, sd=10)
    pred = pm.Bernoulli("pred", p=sigmoid(const), observed=data)

    map_estimate, optresult = pm.find_MAP(
        progressbar=False,
        method="L-BFGS-B",
        return_raw=True,
    )
    
    print(optresult)
    
    map_estimate = float(map_estimate["const"])

print()
print("Avg: ", np.mean(data))
print("MAP: ", map_estimate)
print("sigmoid(MAP): ", sigmoid(map_estimate))

I feel like this should be an exceedingly “easy” problem to solve, but the optimizer returns without success and the returned parameter is way off.

      fun: 3.3127440492577924
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.1347064])
  message: b'ABNORMAL_TERMINATION_IN_LNSRCH'
     nfev: 45
      nit: 4
   status: 2
  success: False
        x: array([-3.53465344])

Avg:  0.11666666666666668
MAP:  -3.5346534415193562
sigmoid(MAP):  0.028342155477815184

Are there best practices for estimating MAP for these cases? Am I missing something obvious? Tinkering with the sd of the prior and optimization method have not proven fruitful.

2

Bernoulli only accept binary data like [0, 1, 0, 0, …]

1 Like
3

Ah, maybe that’s it! Thanks for the quick reply. What distribution would you suggest for modeling probabilities (real numbers in range 0-1)?

4

Solved it myself, switched to a beta likelihood function and it’s all working nicely. Thanks for your guidance @junpenglao!!

data = np.array([0.05, 0.1, 0.1, 0.15, 0.2, 0.1])

with pm.Model():
    const = pm.Normal("const", mu=0, sd=10)
    mu = sigmoid(const)
    sd = mu * (1 - mu)  # bernoulli SD per pymc3 docs
    pred = pm.Beta("pred", mu=mu, sd=sd, observed=data)

    map_estimate, optresult = pm.find_MAP(
        progressbar=False, method="L-BFGS-B", return_raw=True
    )

    print(optresult)

    map_estimate = float(map_estimate["const"])
1 Like