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]) 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))], cmap="inferno",shade=True) trace(t,0) trace(t,1) trace(t,2)
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.