Is my Bayesian analysis correct? This is my first time doing a Bayesian analysis, so I'm not sure whether what I did makes perfect sense. I'm trying to tell if two samples come from the same distribution, more specifically, if they have the same mean.
I have very few samples, and no subjective beliefs about them. I am modelling the observed distribution as a Gaussian. Using pyMC:
import pymc as pm

sample1 = [8.254828927, 8.485694524,10.58058423,5.221356325]
sample2 = [7.921107724, 9.744301571, 5.750874082,8.421883012,11.09813068]

sample_mean = sum(sample1+sample2)

#Empirical prior for the mean
normal_mean_prior = pm.Normal("mean_prior", mu = sample_mean, tau = 1)

#Uninformative prior for the precision
normal_precision_prior = pm.Uniform("precision_prior", lower = 0, upper = 100)

#Normal observation
normal_observation_sample1 = pm.Normal("normal1", mu = normal_mean_prior, tau = normal_precision_prior, observed = True, value = sample1)
normal_observation_sample2 = pm.Normal("normal2", mu = normal_mean_prior, tau = normal_precision_prior, observed = True, value = sample2)

With this model, I would sample it using MCMC and check if the posterior distribution of the means are different or the same.
Is my analysis correct? Have I made any significant blunder? Thank you.
 A: Your prior on the mean is too narrow.  Also, you don't just want samples for the two means, you want samples of the difference between the means.  The posterior distribution of the difference is what you want to look at.  For more details, see the paper Bayesian estimation supersedes the t test.
A: I would say that you have an good start, but your analysis is incomplete.  How did you pick your hyper-parameters?  How will you check if the posterior distribution of the means are the same?  The devil is in the details.
@TomMinka has answered with a reference to a paper by Kruschke that proceeds in the general direction you have been going and fills in the details.  Here is the PyMC implementation of something along those lines (but using the normal instead of t for the likelihood to keep it simple):
def model():
    sample_0 = [8.254828927, 8.485694524, 10.58058423, 5.221356325]
    sample_1 = [7.921107724, 9.744301571, 5.750874082, 8.421883012, 11.09813068]
    pooled_data = np.hstack((sample_0, sample_1))

    # hyper-parameters
    S = 1000 * np.std(pooled_data)
    M = np.mean(pooled_data)

    L = .001 * np.std(pooled_data)
    H = 1000 * np.std(pooled_data)

    # prior
    mu = pm.Normal('mu_1', M, S**-2, value=[M,M])
    sigma = pm.Uniform('sigma_1', L, H, value=[np.std(pooled_data), np.std(pooled_data)])

    # likelihood
    obs_1 = pm.Normal('obs_1', mu=mu[0], tau=sigma[0]**-2, value=sample_0, observed=True)
    obs_2 = pm.Normal('obs_2', mu=mu[1], tau=sigma[1]**-2, value=sample_1, observed=True)

    # derived quantities of interest
    d_mu = mu[0] - mu[1]
    d_sigma = sigma[0] - sigma[1]
    @pm.deterministic
    def effect_size(d_mu=d_mu, sigma=sigma):
        return d_mu / np.sqrt(.5*(sigma[0]**2 + sigma[1]**2))
    @pm.deterministic
    def outside_rope(effect_size=effect_size):
        return -.1 < effect_size < .1

    return locals()

And here is an IPython Notebook to sample from the approximation posterior distribution with MCMC to see that the probability that the two distributions have different means is 0.1.
