Re: [Help-gsl] Non-linear Multi-parameter Least-Squares Example

[John Gehman]

G'day Gordan,

The Jacobian matrix J is the differential of your objective equation  
with respect to each parameter (corresponding  to the second index)  
at each data point (corresponding to the first index).

The relevant equation is  [a * sin(b / t + c) + d - yi ] / sigma[i],  
i.e. you're trying to minimize the difference between the model and  
the data (yi), where some points may be more reliable than others  
(given sigma[i]),  by optimizing the amplitude, period, phase, and  
offset of your sine wave. The Jacobian provides the analytic  
equations to inform the system of the sensitivity of the residuals in  
the evolving fit to changes in each of the parameters.

Taking the partial derivative of your equation with respect to each  
of the floating parameters, and setting them appropriately into the  
matrix J, assuming your vector of floating parameters elsewhere is  
ordered (a,b,c,d), and retaining your s = sigma[i], I get:


gsl_matrix_set (J, i, 0, sin(b / t + c) / s);	        # derivative  
with respect to a
gsl_matrix_set (J, i, 1, cos(b / t + c) * a/(t*s) ); # derivative  
with respect to b
gsl_matrix_set (J, i, 2, cos(b / t + c) * a/s);       # derivative  
with respect to c
gsl_matrix_set (J, i, 3, 1/s);                                #  
derivative with respect to d

The analytic derivatives are usually my problem, however, so please  
confirm them!

Regards,
john

---------------------------------------------------------

Dr John Gehman
Research Fellow
School of Chemistry
University of Melbourne (Australia)


On 17/07/2007, at 9:46 PM, Gordan Bobic wrote:

> Hi,
>
> I have a question regarding the example 37.9 in the manual.
>
> In the df function, there is the following code:
>
> for (i = 0; i < n; i++)
> {
>         /* Jacobian matrix J(i,j) = dfi / dxj, */
>         /* where fi = (Yi - yi)/sigma[i],      */
>         /*       Yi = A * exp(-lambda * i) + b  */
>         /* and the xj are the parameters (A,lambda,b) */
>
>         double t = i;
>         double s = sigma[i];
>         double e = exp(-lambda * t);
>
>         gsl_matrix_set (J, i, 0, e/s);
>         gsl_matrix_set (J, i, 1, -t * A * e/s);
>         gsl_matrix_set (J, i, 2, 1/s);
> }
>
> I am trying to adapt this to fit a generic sine curve. So, I am  
> trying to
> do something like this:
>
> for (i = 0; i < n; i++)
> {
>         /* Jacobian matrix J(i,j) = dfi / dxj, */
>         /* where fi = (Yi - yi)/sigma[i],      */
>         /*       Yi = a * sin(b / t + c) + d   */
>         /* and the xj are the parameters (a,b,c,d) */
>
>         double t = i;
>         double s = sigma[i];
>         double e = sin(b / t + c);
>
>         gsl_matrix_set (J, i, 0, e/s);
>         gsl_matrix_set (J, i, 1, cos (b / t + c) * e/s);
>         gsl_matrix_set (J, i, 2, 1/s);
>         gsl_matrix_set (J, i, 3, 1/s);
> }
>
> Is this a correct adaptation, or am I misunderstanding what this is
> supposed to do? Is the 4th parameter in this case the
> differential/derivative of the Yi equation with respect to a, b, c  
> and d
> respectively?
>
> Thanks.
>
> Gordan
>
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl