Bayesian equivalent of two sample t-test? I'm not looking for a plug and play method like BEST in R but rather a mathematical explanation of what are some Bayesian methods I can use to test the difference between the mean of two samples. 
 A: With a Bayesian analysis you have more things to specify (that is actually a good thing, since it gives much more flexibility and ability to model what you believe the truth to be).  Are you assuming normals for the likelihoods?  Will the 2 groups have the same variance?
One straight forward approach is to model the 2 means (and 1 or 2 variances/dispersions) then look at the posterior on the difference of the 2 means and/or the Credible Interval on the difference of the 2 means.
A: 
a mathematical explanation of what are some Bayesian methods I can use to test the difference between the mean of two samples.

There are several approaches to "testing" this. I'll mention a couple:

*

*If you want an explicit decision you could look at decision theory.


*A pretty simple thing that's sometimes done is to find an interval for the difference in the means and consider whether it includes 0 or not. That would involve starting with a model for the observations, priors on the parameters and computation of the posterior distribution of the difference in means conditional on the data.
You'd need to say what your model is (e.g. normal, constant variance), and then (at least) some prior for the difference in means and a prior for the variance. You might have priors on the parameters of those priors in turn. Or you might not assume constant variance. Or you might assume something other than normality.
A: The excellent answer by user1068430 implemented in Python
import numpy as np
from pylab import plt

def dnorm(x, mu, sig):
    return 1/(sig * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sig**2))

def dexp(x, l):
    return l * np.exp(- l*x)

def like(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(sample1, mu1, sig1).prod()*dnorm(sample2, mu2, sig2).prod()
    
def prior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return dnorm(mu1, pooled.mean(), 1000*pooled.std()) * dnorm(mu2, pooled.mean(), 1000*pooled.std()) * dexp(sig1, 0.1) * dexp(sig2, 0.1)

def posterior(parameters):
    [mu1, sig1, mu2, sig2] = parameters
    return like([mu1, sig1, mu2, sig2])*prior([mu1, sig1, mu2, sig2])


#create samples
sample1 = np.random.normal(100, 3, 8)
sample2 = np.random.normal(100, 7, 10)

pooled= np.append(sample1, sample2)

plt.figure(0)
plt.hist(sample1)
plt.hold(True)
plt.hist(sample2)
plt.show(block=False)

mu1 = 100 
sig1 = 10
mu2 = 100
sig2 = 10
parameters = np.array([mu1, sig1, mu2, sig2])

niter = 10000

results = np.zeros([niter, 4])
results[1,:] = parameters

for iteration in np.arange(2,niter):
    candidate = parameters + np.random.normal(0,0.5,4)
    ratio = posterior(candidate)/posterior(parameters)
    if np.random.uniform() < ratio:
        parameters = candidate
    results[iteration,:] = parameters

#burn-in
results = results[499:niter-1,:]

mu1 = results[:,1]
mu2 = results[:,3]

d = (mu1 - mu2)
p_value = np.mean(d > 0)

plt.figure(1)
plt.hist(d,normed = 1)
plt.show()
```

