# Best practices for MAP estimation

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.

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

1 Like

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

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