Vrftsjtryk

I have the 2 parameters arrays : $theta=[a,b]$, $nu=[r_0,c_0,alpha,beta]$ with the distribution (a point spread function = PSF = response of a Dirac) :

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

I would like to convert this relation under matricial form.

Before, I used the matrix form with only $theta$ and I could write :

$$d=H.theta + epsilonquad(1)$$

i.e the equality :

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

I don't know how to keep the same matricial form but with also $nu=[r_0,c_0,alpha,beta]$ and $theta=[a,b]$ ?

This matricial form will allow me to use Maximum-Likelihood-Estimation on all parameters.

But there, I wonder if I have just to use the same matricial form above $(1)$ but with explicit $H=H(nu)$ ? then we could write :

$ d = H(nu) cdot theta + epsilon quad(1)$

So after, I would generate random values for $[r_{0}, c_{0}, alpha, beta]$ and $[a,b]$, wouldn't it ?

PS : I have posted a first question on a different stack forum but it is more about processing image in astrophysics.

UPDATE 1 : From what I have seen, it seems the Expectation Maximization (EM) is appropriate for my problem.

The likelihood function is expressed as :

$mathcal{L}=prod_{i=1}^{n} text{PSF}_{i}$

So I have to find $nu$ and $theta$ such that $dfrac{partial text{ln}mathcal{L}}{partial nu}=0$ and $dfrac{text{ln}mathcal{L}}{partial theta}=0$

Could anyone help me to implement the EM algorithm to estimate $nu$ and $beta$ arrays of parameters ?