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:
  elif j >50:

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

binary data

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):
    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

      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))]

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 plot, b0 (intercept) trace1

Trace plot, b1 (slope) trace2

Trace plot, s (sigma) trace3

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:
  return sum(scores)/len(scores)

>>> 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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – user20160
    Sep 17 '20 at 3:00
  • $\begingroup$ Are you referring to the formula I’ve listed or not? $\endgroup$
    – jbuddy_13
    Sep 17 '20 at 3:02
  • 2
    $\begingroup$ 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. $\endgroup$
    – user20160
    Sep 17 '20 at 3:53
  • $\begingroup$ Should elif j >50: be elif j >-50:? $\endgroup$
    – Xi'an
    Sep 17 '20 at 7:12
  • $\begingroup$ @user20160, ah yes - as I noted in my post, I did use that as well, with terrible results. I'll link the code. $\endgroup$
    – 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).

  • $\begingroup$ I edited the post significantly to reflect my prior/likelihood selections, which have failed spectacularly. What's the naïve/generic prior/likelihood combination that you'd use for educational purposes? $\endgroup$
    – jbuddy_13
    Sep 17 '20 at 1:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.