for the general case : data is normal-gamma (mean normal and sd is gamma) and I want to estimate the $b$ and $a$ distribution ( they assumed to be normal) in $y=bx+a$ using Bayesian regression. I know I need assumption for my priors but I think I have to be OK with most simplistic case. not very good in reading statistics literature to derive the code from heavy literature.
{trying to be more clear} piece of code is helpful as well as pseudo-code of how to do it. Here I explain the case I did for single parameter as an example. there are 2 ways actually by conditioning on ro/p (precision) or unconditional. Keep these pre-calculation:
Sxx<-sum((x-xbar)^2) , Syy<-sum((y-ybar)^2) , Sxy<-sum((x-xbar)(y-ybar)) , See<-sum( ( y-ybar-b(x-xbar))^2)}
Posterior distribution for $\rho$ is Gamma with shape $(n-2)/2$ and scale $2/S_{ee}$.
Conditional on $\rho$, $\eta$ and $\beta$ are independent and normally distributed, $\eta$ has mean $y$ and precision $n\rho$, $\beta$ has mean $b$ and precision $S_{xx}\rho$.
Uncoditional: Marginal posterior distribution for $\eta$ and $\beta$. $\eta$ has nonstandard t distribution with: Center $\bar y$, spread $(n(n-2)/S_{ee})-1/2$, degrees of freedom $n-2$.
$\beta$ has nonstandard t distribution with: Center $\bar y$, spread $(S_{xx}(n-2)/S_{ee})-1/2$, degrees of freedom $n-2$.
$\eta$ and $\beta$ are uncorrelated but not independent when we marginalize out $\rho$.
I dont know how this could be for MultiVariate Case. I know I can use Winbugs but I need to write it in Matlab.
