3
$\begingroup$

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$ :

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.

Regards

$\endgroup$

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.