# Posterior mean estimator with MCMC (Metropolis Hastings Algorithm) - Concrete example

I have a little project for which I have to estimate parameters on a PSF (Point Spread Function = response of the system to a dirac, i.e a star in my case).

I have the 6 parameters to estimate : $$p=(\theta=[a,b]$$, $$\nu=[r_0,c_0,\alpha,\beta])$$ :

$$\text{PSF}(r,c) = \bigg(1 + \dfrac{r^2 + c^2}{\alpha^2}\bigg)^{-\beta}$$

and the modeling :

$$d(r,c) = a \cdot \text {PSF}_{\alpha,\beta} (r-r_0,c-c_0) + b + \epsilon (r, c)$$

with $$\epsilon$$ being a white gaussian noise.

This is the matrix form :

$$\begin{bmatrix}d(1,1) \\ d(1,2) \\ d(1,3) \\ \vdots \\ d(20,20) \end{bmatrix} = \begin{bmatrix} \text{PSF}_{\alpha,\beta}(1-r_0,1-c_0) & 1 \\ \text{PSF}_{\alpha,\beta}(1-r_0,2-c_0) & 1 \\ \text{PSF}_{\alpha,\beta}(1-r_0,3-c_0) & 1 \\ \vdots & \vdots \\ \text{PSF}_{\alpha,\beta}(20-r_0,20-c_0) & 1 \end{bmatrix} \times \begin{bmatrix}a \\ b \end{bmatrix} + \begin{bmatrix}\epsilon(1,1) \\ \epsilon(1,2) \\ \epsilon(1,3) \\ \vdots \\ \epsilon(20,20) \end{bmatrix}\\$$

So we can write for the 1D data vector "d" :

$$d=H(\nu)\,\theta + \epsilon$$ with $$H$$ the matrix defined above.

Our Teacher told us to start from a posterior distribution $$f(p|d)$$ (from what I have understood, the density function to get $$p$$ parameters knowing the data).

I have to give this expression and its logarithmic form $$\text{log(f(p|d)}$$ as a function of cost function $$J_{\theta,\nu}$$ (I have already implement this cost function but for Maximum Likelihood and after, I used "fminsearch" Matlab function to find best parameters).

Question 1) Someone could help me to find the expression of this posterior ditribution $$f(p|d)$$ in my case ?

Question 2) I have to use the Metropolis Hastings algorithm to find the estimations of the 6 parameters above ( $$p=\theta=[a,b]$$, $$\nu=[r_0,c_0,\alpha,\beta]$$)

I know this algorithm can generate samples following distributions but how to make the link with the estimations of parameters contained into $$d$$ parameters vector ($$p=[a,b,r_0,c_0,\alpha,\beta]$$ ?

Feel free to ask me further informations if you don't understand my concrete example. Technically, it doesn't seem to be very complex but I am a beginner into estimations parameters with MCMC methods.

Any help or suggestion is welcome.

UPDATE 1 :

Question 1) I have found the expression of $$f(p|d)$$ :

$$f(p|d) = \dfrac{f(p)\,f(d|p)}{\int_{p}f(d|p)\,\text{d}p}$$

We can also write :

$$f(p|d) \propto \text{Likelihood(d|p)}\,f(p)\quad(1)$$

with $$f(p)$$ the prior function (that I can take as uniform distribution).

Here's an example of code (R language) from this link with the estiation of a parameter $$\theta$$ :

x <- rexp(100)

N <- 10000
theta <- rep(0,N)
theta[1] <- cur_theta <- 1  # Starting value
for (i in 1:N) {
prop_theta <- runif(1,0,5)  # "Independence" sampler
alpha <- exp(sum(dexp(x,prop_theta,log=TRUE)) - sum(dexp(x,cur_theta,log=TRUE)))
if (runif(1) < alpha) cur_theta <- prop_theta
theta[i] <- cur_theta
}

hist(theta)


and the histogram gives the best estimation of parameter $$\theta$$ :

From what I have found on stackexchange forums : I can take a uniform prior density for $$f(p)$$ and After, I can do sampling following the distribution $$f(d|p)$$.

If I take the product $$f(d|p)\times f(p)$$ and once $$p$$ is given, I plot the histogram for a large number of parameters generated $$p$$ (correspondig to $$f(d|p$$), we can see the peak on histogram which will indicate me the best parameter which fits the data $$d$$ : is the method correct ?

Last question : Why does our teacher asked us to use the cost function used in previous Likelihood method used whith "fminsearch" ?

UPDATE 2: Maybe I did confusions in the expression of posterior probability. Indeed, concerning the Bayes theorem, I saw this :

1) Is $$p(data)$$ a normalization factor ? (is it taken as $$p(data)=1$$ ?)

2) Is $$p(p)=p(\theta)$$ always taken as uniform distribution (maybe to simplify the expression of $$p(p|d)$$) ?

3) How to justify that we can chose a likelihood for $$p(d|p)$$ ? and morever, from what I have seen, it is often chosen as a gaussian likelihood, isn't it ?

UPDATE 3: Are my questions bad formulated ? since I don't manage to get answers/advices. I hope that someone will help me.

Regards