# What are some common prior/likelihood choices for Bayesian logistic regression?

I'm not really clear on the Bayesian approach to logistic regression. From everything I've read, the prior and likelihood can be can be whatever you want them to be. Well, I've a couple things; namely, gaussian likelihood and p**y + (1-p)**(1-y) which I'm familiar with from the ML approach to logistic regression.

Long story short, both failed miserably. To answer this question, please address the following: (1) Why is this combination of likelihood and prior not adequate? (2) Is my synthetic data not appropriate for logistic regression? (3) Are there any obvious errors in my implementation (via MH algo)?

x_samples = []
y_samples = []
for i in range(1,500):
j = i -250
if j <-50:
y_samples.append(0)
elif j >50:
y_samples.append(1)

else:
p0 = i/500
p1 = 1 - p0
c = random.choices(population=[0,1],weights=[p0,p1])[0]
y_samples.append(c)
x_samples.append(i)


import numpy as np
import matplotlib.pyplot as plt
def sigmoid(z):
return 1/(1+np.exp(-z))

def normalPDF(x,mu,sigma):
num = np.exp(-1/2*((x-mu)/sigma)**2)
den = np.sqrt(2*np.pi)*sigma
return num/den

def invGamma(x,a,b):
non_zero = int(x>=0)
func = x**(a-1)*np.exp(-x/b)
return non_zero*func

def logreg_mcmc(X,Y,hops=10_000):
samples = []
curr_a = 1
curr_b = 1
curr_s = 1

prior_curr_a = normalPDF(x=curr_a,mu=1,sigma=1)
prior_curr_b = normalPDF(x=curr_b,mu=1,sigma=1)
prior_curr_s = invGamma(x=curr_s,a=3,b=1)

P = [sigmoid(z= curr_a * x + curr_b) for x in X]
curr_log_lik = np.sum([normalPDF(x=pi,mu=yi,sigma=curr_s) for pi,yi in zip(P,Y)])
current = curr_log_lik + np.log(prior_curr_a) + np.log(prior_curr_b) + np.log(prior_curr_s)

count = 0
for i in range(hops):
samples.append((curr_a,curr_b,curr_s))

if count == 0:
mov_a = curr_a + random.uniform(-0.25,0.25)
mov_b = curr_b
mov_s = curr_s
count +=1

elif count == 1:
mov_a = curr_a
mov_b = curr_b + random.uniform(-0.25,0.25)
mov_s = curr_s
count += 1

else:
mov_a = curr_a
mov_b = curr_b
mov_s = curr_s + random.uniform(-0.25,0.25)
count = 0
if mov_s <= 0:
continue # auto-reject

prior_mov_a = normalPDF(x=mov_a,mu=1,sigma=1)
prior_mov_b = normalPDF(x=mov_b,mu=1,sigma=1)
prior_mov_s = invGamma(x=mov_s,a=3,b=1)

P = [sigmoid(z= mov_a * x + mov_b) for x in X]
mov_log_lik = np.sum([normalPDF(x=pi,mu=yi,sigma=mov_s) for pi,yi in zip(P,Y)])

movement = mov_log_lik + np.log(prior_mov_a) + np.log(prior_mov_b) + np.log(prior_mov_s)

ratio = np.exp(movement - current)
event = random.uniform(0,1)
if event <= ratio:
curr_b = mov_b
curr_a = mov_a
curr_s = mov_s
current = movement

return samples


Trace plots to examine MH exploration.

t = logreg_mcmc(X=x_samples,Y=y_samples,hops=25_000)

def trace(data,index):
x = [i for i in range(len(data))]
y = [data[i][index] for i in range(len(data))]
plt.plot(x,y)

def postDensity(data,index1,index2):
g = sns.kdeplot([data[i][index1] for i in range(len(data))],[data[i][index2] for i in range(len(data))],
trace(t,0)
trace(t,1)
trace(t,2)


Trace plot, b0 (intercept)

Trace plot, b1 (slope)

Trace plot, s (sigma)

Performance on first two looks a bit messy. The third looks very skeptical. And the posterior density below:

From the plot (slope=x-axis, y-intercept=y-axis) you can see that density peaks around ~(1,1). Now to test performance:

def accuracy(X,C):
scores = []
for x,c in zip(X,C):
z = 1*x + 1
p = sigmoid(z)
pred = int(p>0.5)
if pred == c:
scores.append(1)
else:
scores.append(0)
return sum(scores)/len(scores)

accuracy(x_samples,y_samples)
>>> 0.503006012024048


So the accuracy is essentially random guessing. To recap, I'm not sure if the problem(s) are my synesthetic data, my priors and/or likelihoods, or my implementation. I did mention that I did I did a similar process with the likelihood p**y + (1-p)**(1-y) which also failed miserably. This post is getting quite long, so I'll only add that in upon request.

• Logistic regression implies a particular likelihood function (which is the same as used in non-Bayesian logistic regression). There's no choice about it. Naturally, you could use a different likelihood function, but this would no longer be logistic regression. But, you have complete freedom in choosing the prior. – user20160 Sep 17 '20 at 3:00
• Are you referring to the formula I’ve listed or not? – jbuddy_13 Sep 17 '20 at 3:02
• The likelihood function for logistic regression is $p(y \mid X, w) = \prod_{i=1}^n \sigma(w^T x_i)^{y_i} \big( 1 - \sigma(w^T x_i) \big)^{1-y_i}$ where $y = \{y_1, \dots, y_n\}$ are binary responses, $X = \{x_1, \dots, x_n\}$ are input points, $w$ is a weight vector, and $\sigma(\cdot)$ is the logistic sigmoid function. – user20160 Sep 17 '20 at 3:53
• Should elif j >50: be elif j >-50:? – Xi'an Sep 17 '20 at 7:12
• @user20160, ah yes - as I noted in my post, I did use that as well, with terrible results. I'll link the code. – jbuddy_13 Sep 17 '20 at 12:09

The likelihood function is the density of the data considered as a function of the parameter $$L:\,\theta\in\Theta\longmapsto L(\theta;x_1,\ldots,x_n)=\prod_{i=1}^n f_i(x_i;\theta)$$ assuming the random variables $$X_i$$ in the sample are independent. In the event the $$X_i$$'s have a finite or countable support, the density is usually defined wrt to the counting measure over that support and this means using the pmf. Note however that the likelihood function can be defined wrt any $$\theta$$-free dominating measure agreeing with the distribution of the rv's. There is no debate whatsoever about the "choice" of the likelihood function in the literature (and hence no such thing as a "good" likelihood).