Multidimensional Gaussian integral with constraints

foobster

1. The problem statement, all variables and given/known data
The larger context is that I'm looking at the scenario of fitting a polynomial to points with Gaussian errors using chi squared minimization. The point of this is to describe the likelihood of measuring a given parameter set from the fit. I'm taking N equally spaced x values and saying that the probability of measuring y at each x value is described by a normal distribution centered around the value of some parent function. I then compute the value of chi squared leaving the y values as variables and assuming a fit function, for example a + b*x + c*x*x. Taking the partial derivatives with respect to the fit parameters gives me a set of constraint equations.

Now I'm trying to integrate over all of the possible y configurations that satisfy my constraints so I can get a likelihood depending only on the fit parameters and the parameters of the original function that generated the y values. All of the y terms can be easily substituted out using the constraint equations except for the Ʃy^2 term in the exponential. If somebody could explain how to do this even with the simplest case of a constraint like [itex]a_{1}y_{1}+a_{2}y_{2}+a_{3}y_{3}+...=A[/itex] then I would really appreciate it. In this case I tried solving for [itex]y_{1}[/itex], substituting it into the exponential, then integrating over all the other ys but I got confused about how to deal with the correlated terms in a general way. Either specific help on this or a pointer to more information on this type of problem would be great.

On a side question, does anybody know of any software that can help deal with symbolic algebra and calculus in arbitrary dimensions? I frequently find myself trying to work with equations and integrals that have some undefined N terms or dimensions and haven't been able to figure out how to do this in Mathematica.