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.


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


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

enter image description here

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 :

enter image description here

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.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.