I solved the problem. I post the solution here in the hope it will be useful to others.

Again, the problem formulation is:

compute the likelihod of observing sample mean \(\bar\mu\) and sample standard deviation \(\bar\sigma\) in \(n\) samples drawn from the distribution \(N(\mu,\sigma^2)\)

The quantity we are interested in is \(\log P(\bar \mu, \bar\sigma | \mu, \sigma)\) (I use log-likelihood since using log simplifies notation). Since sample mean \(\bar \mu\) and sample variance \(\bar\sigma^2\) of a normally distributed population are two independent random variables, we have that

\(

\log P(\bar \mu, \bar\sigma | \mu, \sigma) = \log P(\bar \mu | \mu, \sigma) + \log P(\bar\sigma | \mu, \sigma)

\)

The random variable \(M=\frac{\bar\mu-\mu}{\sigma/\sqrt{n}}\) is t-distributed with \(n-1\) degrees of freedom. For large \(n\), Student's t-distribution tends to \(N(0,1)\); we assume big n so we can approximate the distribution of \(M\) with \(N(0,1)\). Then (applying the definition of the standard normal distribution's density function),

\(

\log P(\bar \mu | \mu, \sigma) \approx \log\left(\frac{1}{\sqrt{2\pi}}exp(-M^2/2)\right) = - \frac{1}{2}\log{2\pi} - \frac{(\bar\mu-\mu)^2n}{2\sigma^2}

\)

The random variable \(S=\frac{(n-1)\bar\sigma^2}{\sigma^2}\) is chi-distributed with \(n-1\) degrees of freedom. Again, we assume \(n\) to be large. Then, the distribution of the random variable \(Q=\frac{S-n}{\sqrt{2n}}\) tends to \(N(0,1)\) and we have:

\(

\log P(\bar\sigma | \mu, \sigma) \approx \log\left( \frac{1}{\sqrt{2\pi}}exp(-Q^2/2) \right) = - \frac{1}{2}\log{2\pi} - \frac{\left((n-1)\bar\sigma^2-n\sigma^2\right)^2}{4n\sigma^4}

\)

note that \(n\approx n-1\) (n is large), so the above quantity simplifies to

\(

- \frac{1}{2}\log{2\pi} - \frac{n^2\left(\bar\sigma^2-\sigma^2\right)^2}{4n\sigma^4} = - \frac{1}{2}\log{2\pi} - \frac{n\left(\bar\sigma^2-\sigma^2\right)^2}{4\sigma^4}

\)

Putting it all together, we finally obtain

\(

\log P(\bar \mu, \bar\sigma | \mu, \sigma) = - \log{2\pi} - \frac{n}{2\sigma^2}\left( (\bar\mu-\mu)^2 + \frac{(\bar\sigma^2-\sigma^2)^2}{2\sigma^2} \right)

\)