Hi Brian, I use the weights both in r[i] and Jacobians. What I realized after hours of playing with my code is that the coefficient values are quite different than the known solution from a different software (on the same input data). So, there is probably some other problem in my code and re-weighting only makes it more important and visible. Now I have a different set of input data for which re-weighting does not trigger the solver error. In this case, re-weighting makes final values of coefficients slightly different with a better R-squared value than with fixed weights. But still it is far away from coefficients computed by another software.
I really tried find the problem in my code, but with no success yet. The only thing, which I think is really different than in the example in the gsl documentation (with static weights, of course), it seems to be the "data" structure, which is used as a second parameter of func_f and func_df functions. I hope thiese structure data are used only in the code of these functions and not by any other internal solver procedure, so this structure can have all the members I need (not only n, y and sigma) in any order. Is it right ? If yes, then I really do not know what could be wrong, but it would mean the problem is somewhere in my head :-) Thanks again. Sincerely, David Dne 10.11.2010 6:44, Brian Hawkins napsal(a): > Hi David, > > I don't think there's anything wrong with the approach of re-computing the > weights. The solver has no knowledge of the weights, only the functions you > provide. Keep in mind that you should be supplying the residual function, > e.g. > > r[i] = (y[i] - f(p, x[i])) * weight[i] > > and similarly the rows of the Jacobian are scaled by the weights. While I > think varying the weights between iterations is probably okay, you should > verify that they are well-behaved. As a first cut, none are going to 0 or > inf relative to the other weights and/or machine precision. I don't know to > what extent the GSL solvers are robust to poorly-scaled systems. Also if > your weights vary strongly between iterations, that would suggest you're > taking too-large a step or maybe the weight model is wrong. > > The part of the GSL lmsder algorithm that is somewhat mysterious to me is > the "trust region" aspect. In solvers I've written in the past, I > essentially hand-tuned the step size when necessary. You might check to see > that you're taking reasonably sized steps (s->dx). I don't know if there's > any interface for manual control of the step size (didn't see it in the > book). > > Of course, try to make sure the parameters are observable in your data set. > By that I mean you can see a pattern in the data that suggests you're > fitting an appropriate model. I'm not sure offhand what your model or data > look like. Maybe you're only testing with a subset of your data, and it's > not enough. > > It's also possible that your model is just really tough to fit. For > example, your performance could be very sensitive to the initial guess. > What happens when you make your initial guess is very close to a known > solution (on known/simulated data)? > > Hope that helps, > Brian > > Message: 1 >> Date: Sun, 07 Nov 2010 22:02:16 +0100 >> From: David Komanek <[email protected]> >> Subject: Re: [Help-gsl] Re: iteratively re-weighted least squares >> fitting >> To: [email protected] >> Message-ID: <[email protected]> >> Content-Type: text/plain; charset=ISO-8859-2 >> >> Dear Brian, >> >> here is the main loop: >> >> computeWeightsByYi(); // computes weights for the first round >> do { >> iter++; >> status = gsl_multifit_fdfsolver_iterate(s); >> if (status) break; >> status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); >> if (movingWeights) { // boolean to switch re-weighted and >> "normal" regression >> gsl_vector *params = gsl_vector_alloc(s->x->size); >> gsl_vector_memcpy(params, s->x); >> setParameters(params); // copies current model parameters >> into the enclosing object member data field >> computeWeightsByFx(); // computes new weights from the new >> set of parameters for the next iteration >> } >> } >> while ((status == GSL_CONTINUE) && (iter < 10000)); >> >> The solver is gsl_multifit_fdfsolver_lmsder. >> >> The function is f(x) = a * exp(b * (x+0.5)) / (1 + a * exp(b * (x+0.5))) >> >> With nearly certainty, the problem is on my side, but I wanted to be >> sure that re-weighting is o.k. with GSL or that I will need another >> approach. I am not a matematician, but I am pretty sure the partial >> derivatives etc. are o.k. in my case - my colleague uses another >> computing method with a success, but he has it as a set of VBA macros >> in Excel and we need to do it in C/C++. For some reason I think GSL will >> be better choice than blindly rewrite code lines from VBA to C++. What I >> like to know is if I am using the GSL right for this type of problem. >> >> Thank you. Kind regards, >> >> David >> >> . >> Dne 6.11.2010 18:46, Brian Hawkins napsal(a): >>> David, >>> >>> What's your convergence criterion? Is your system full-rank? Have you >> had >>> success with this problem using a different solver? I'm having a hard >> time >>> understanding why the GSL solver in particular would be giving you >> trouble. >>> Regards, >>> Brian >>> >>> On Sat, Nov 6, 2010 at 9:01 AM, <[email protected]> wrote: >>>> Message: 1 >>>> Date: Fri, 05 Nov 2010 23:09:54 +0100 >>>> From: David Komanek <[email protected]> >>>> Subject: [Help-gsl] iteratively re-weighted least squares fitting >>>> To: [email protected] >>>> Message-ID: <[email protected]> >>>> Content-Type: text/plain; charset=ISO-8859-2 >>>> >>>> Dear all, >>>> >>>> I ran into problems using weighted least squares fitting in GSL. For >>>> some reason I need to use IRLS modification of this method, so the >>>> weights are recomputed in every iteration. In the case the weights are >>>> computed at the beginning and being constant throug all the iterations, >>>> the procedure works fine. But when I adjust the weights in every >>>> iteration, this usually leads to an error: >>>> >>>> 27 iteration is not making progress towards solution >>>> >>>> I think it is because there are some internally tested conditions and >>>> some of them are not satisfied in this case. For example, in SAS, there >>>> is a special parameter to relax those conditions: >>>> >>>> >>>> >> http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_nlin_sect034.htm >>>> Is something like this possible with GSL ? >>>> >>>> Thank you in advance. >>>> >>>> David >>>> > _______________________________________________ > Help-gsl mailing list > [email protected] > http://lists.gnu.org/mailman/listinfo/help-gsl _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
