Prior estimations for mu and sigma^2 are

mu~ N(ybar, s^2/n), and sigma^2~(n-1)s^2 chi^-2(n-1)

where s^2 is the (n-1) version sample variance and chi^-2 is the inverse chi^2 dist.

random numbers can be generated from the fact that if x~ chi^2 then 1/x is the inverse chi^2.

I have code to calc xbar, and xvar (which is s^2 i think), but am unsure how to change the standard normal dist to a random normal dist. Especially for sigma^2. Any ideas?

I think I can switch to gamma dist maybe?

xbar=mean(x);

s2 = xvar; %? n-1 sample variance?

mu= xbar + sqrt((s2)/n) * randn(m,1);

Thanks