Multidimensional Gaussian in gam

Question

I’ve been reading the documentations and discussions I can find but couldn’t get an answer, so I hope someone can help.

  I want to use gam in the context of an optimization algorithm (EM). In a simpler scenario, if I wanted to fit a 1d Gaussian using gam, I could have done

fit <- gam(y ~ 1)

However, in my case there are 2 complications:

1. My y is a 2d Gaussian; and
2. Not only do I need (mu1, sigma1, mu2, sigma2) from the fitting, but also mu1 and mu2 are some non-linear functions of the observables.

I would like to know whether in this case I can still make use of gam, or do I have to create my own function.

Edit:

To be precise, the problem is in the form

n(x, mu1, sigma1) * n(y, x + mu2(beta), sigma2)

where n() is the Gaussian pdf, x and y are independent variables, and mu2 is a (non-linear) function of an independent variable beta.