[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] Re: Numerical Gradient?
From: |
Richard Henwood |
Subject: |
[Help-gsl] Re: Numerical Gradient? |
Date: |
Thu, 06 Mar 2008 09:23:20 +0000 |
User-agent: |
KNode/0.10.4 |
David Doria wrote:
> I have a function f() that is not analytic and therefore I need to find a
> numerical gradient. However to initialize a gsl_multimin_function_fdf, I
> still need a function my_df that simply returns the value of the gradient
> at the specified point.
>
> For testing, lets try it with a function that we COULD have taken a normal
> gradient of (this example ships with the package)
>
> double my_f (const gsl_vector *v, void *params)
> {
> double x, y;
> double *dp = (double *)params;
> x = gsl_vector_get(v, 0);
> y = gsl_vector_get(v, 1);
> return 10.0 * (x - dp[0]) * (x - dp[0]) +
> 20.0 * (y - dp[1]) * (y - dp[1]) + 30.0;
> }
>
> and its gradient:
>
> void my_df (const gsl_vector *v, void *params, gsl_vector *df)
> {
> double x, y;
> double *dp = (double *)params;
> x = gsl_vector_get(v, 0);
> y = gsl_vector_get(v, 1);
> gsl_vector_set(df, 0, 20.0 * (x - dp[0]));
> gsl_vector_set(df, 1, 40.0 * (y - dp[1]));
> }
>
> but instead i would like: (and I imagine I should keep the same function
> form (ie. same parameters))
>
> void my_df (const gsl_vector *v, void *params, gsl_vector *df)
> {
> NUMERICAL GRADIENT OF f at v;
> }
>
> Is there a function I can call that will do that?
>
My experience may be applicable to gsl_multimin_function_fdf.
I have implemented fitting (gsl_multifit_fdfsolver) using numerical
differentials. However, I did not find a function which would generically
differentiate a given function with respect to a given parameter. I wrote
the function I wanted to differentiate out multiple times, each time
exposing the parameter I wanted. see below for a terse example of fitting a
log Normal Probability Density Function:
[NOTE: log normal is not a good case as the analytical partial diffentials
exist!]
double myLogNormalPDF_x ( double x, void * p) {...};
double myLogNormalPDF_mu ( double , void * p) {...};
double myLogNormalPDF_sigma ( double , void * p) {...};
and then:
int myLogNormalPDF_f (const gsl_vector * fitParams, void * extraParams,
gsl_vector * fn) {
...
for (i = 0; i < n; i++) {
current_x = xdata[i];
if (current_x > 0.0) {
current_eval = myLogNormalPDF_x (current_x, ¶ms);
gsl_vector_set (fn, i, current_eval - ydata[i]);
}
}
return GSL_SUCCESS;
}
and:
int myLogNormalPDFnum_df ( const gsl_vector * fitParams, void * fitData,
gsl_matrix * J) {
...
for (i = 0; i < n; i++) {
current_x = xdata[i];
if (current_x > 0) {
params.x = current_x;
params.mu = fn_mu;
params.sigma = fn_sigma;
F.function = &myLogNormalPDF_mu;
gsl_diff_central(&F, fn_mu, &df_dmu, &difAbserr);
gsl_matrix_set(J, i, 0, df_dmu);
F.function = &myLogNormalPDF_sigma;
gsl_diff_central(&F, fn_sigma, &df_dsigma, &difAbserr);
gsl_matrix_set(J, i, 1, df_dsigma);
}
else {
gsl_matrix_set(J, i, 0, 0.0);
gsl_matrix_set(J, i, 1, 0.0);
}
}
return GSL_SUCCESS;
}