Combining information from multiple studies to estimate the mean and variance of normally distributed data - Bayesian vs meta-analytic approaches

I have reviewed a set of papers, each reporting the observed mean and SD of a measurement of $X$ in its respective sample of known size, $n$. I want to make the best possible guess about the likely distribution of the same measure in a new study that I am designing, and how much uncertainty is in that guess. I am happy to assume $X \sim N(\mu, \sigma^2$).

My first thought was meta-analysis, but the models typically employed focus on point estimates and corresponding confidence intervals. However, I want to say something about the full distribution of $X$, which in this case would also including making a guess about the variance, $\sigma^2$.

I have been reading about possible Bayeisan approaches to estimating the complete set of parameters of a given distribution in light of prior knowledge. This generally makes more sense to me, but I have zero experience with Bayesian analysis. This also seems like a straightforward, relatively simple problem to cut my teeth on.

1) Given my problem, which approach makes the most sense and why? Meta-analysis or a Bayesian approach?

2) If you think the Bayesian approach is best, can you point me to a way to implement this (preferably in R)?

Related question

EDITS:

I have been trying to work this out in what I think is a 'simple' Bayesian manner.

As I stated above, I am not just interested in the estimated mean, $\mu$, but also the variance,$\sigma^2$, in light of prior information, i.e. $P(\mu, \sigma^2|Y)$

Again, I know nothing about Bayeianism in practice, but it didn't take long to find that the posterior of a normal distribution with unknown mean and variance has a closed form solution via conjugacy, with the normal-inverse-gamma distribution.

The problem is reformulated as $P(\mu, \sigma^2|Y) = P(\mu|\sigma^2, Y)P(\sigma^2|Y)$.

$P(\mu|\sigma^2, Y)$ is estimated with a normal distribution; $P(\sigma^2|Y)$ with an inverse-gamma distribution.

It took me a while to get my head around it, but from these links(1, 2) I was able, I think, to sort how to do this in R.

I started with a data frame made up from a row for each of 33 studies/samples, and columns for the mean, variance, and sample size. I used the mean, variance, and sample size from the first study, in row 1, as my prior information. I then updated this with the information from the next study, calculated the relevant parameters, and sampled from the normal-inverse-gamma to get the distribution of $\mu$ and $\sigma^2$. This gets repeated until all 33 studies have been included.

# Loop start values values

i <- 2
k <- 1

# Results go here

muL      <- list()  # mean of the estimated mean distribution
varL     <- list()  # variance of the estimated mean distribution
nL       <- list()  # sample size
eVarL    <- list()  # mean of the estimated variance distribution
distL    <- list()  # sampling 10k times from the mean and variance distributions

# Priors, taken from the study in row 1 of the data frame

muPrior  <- bayesDf[1, 14]    # Starting mean
nPrior   <- bayesDf[1, 10]    # Starting sample size
varPrior <- bayesDf[1, 16]^2  # Starting variance

for (i in 2:nrow(bayesDf)){

# "New" Data, Sufficient Statistics needed for parameter estimation

muSamp    <- bayesDf[i, 14]          # mean
nSamp     <- bayesDf[i, 10]          # sample size
sumSqSamp <- bayesDf[i, 16]^2*(nSamp-1)  # sum of squares (variance * (n-1))

# Posteriors

nPost   <- nPrior + nSamp
muPost  <- (nPrior * muPrior + nSamp * muSamp) / (nPost)
sPost   <- (nPrior * varPrior) +
sumSqSamp +
((nPrior * nSamp) / (nPost)) * ((muSamp - muPrior)^2)
varPost <- sPost/nPost
bPost   <- (nPrior * varPrior) +
sumSqSamp +
(nPrior * nSamp /  (nPost)) * ((muPrior - muSamp)^2)
# Update

muPrior   <- muPost
nPrior    <- nPost
varPrior  <- varPost

# Store

muL[[i]]   <-  muPost
varL[[i]]  <-  varPost
nL[[i]]    <-  nPost
eVarL[[i]] <- (bPost/2) / ((nPost/2) - 1)

# Sample

muDistL  <- list()
varDistL <- list()

for (j in 1:10000){
varDistL[[j]] <- 1/rgamma(1, nPost/2, bPost/2)
v             <- 1/rgamma(1, nPost/2, bPost/2)
muDistL[[j]]  <- rnorm(1, muPost, v/nPost)
}

# Store

varDist    <- do.call(rbind, varDistL)
muDist     <- do.call(rbind, muDistL)
dist       <- as.data.frame(cbind(varDist, muDist))
distL[[k]] <- dist