File indexing completed on 2024-05-12 15:27:07
0001 /*************************************************************************** 0002 File : nsl_fit.c 0003 Project : LabPlot 0004 Description : NSL (non)linear fit functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2016-2018 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 ***************************************************************************/ 0008 0009 /*************************************************************************** 0010 * * 0011 * This program is free software; you can redistribute it and/or modify * 0012 * it under the terms of the GNU General Public License as published by * 0013 * the Free Software Foundation; either version 2 of the License, or * 0014 * (at your option) any later version. * 0015 * * 0016 * This program is distributed in the hope that it will be useful, * 0017 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 0019 * GNU General Public License for more details. * 0020 * * 0021 * You should have received a copy of the GNU General Public License * 0022 * along with this program; if not, write to the Free Software * 0023 * Foundation, Inc., 51 Franklin Street, Fifth Floor, * 0024 * Boston, MA 02110-1301 USA * 0025 * * 0026 ***************************************************************************/ 0027 0028 #include "nsl_fit.h" 0029 #include "nsl_common.h" 0030 #include "nsl_sf_basic.h" 0031 #include <gsl/gsl_math.h> 0032 #include <gsl/gsl_sf_erf.h> 0033 #include <gsl/gsl_sf_gamma.h> 0034 #include <gsl/gsl_sf_psi.h> 0035 #include <gsl/gsl_randist.h> 0036 0037 const char* nsl_fit_model_category_name[] = {i18n("Basic functions"), i18n("Peak functions"), i18n("Growth (sigmoidal)"), i18n("Statistics (distributions)"), 0038 i18n("Custom")}; 0039 0040 const char* nsl_fit_model_basic_name[] = {i18n("Polynomial"), i18n("Power"), i18n("Exponential"), i18n("Inverse exponential"), i18n("Fourier")}; 0041 const char* nsl_fit_model_basic_equation[] = {"c0 + c1*x", "a*x^b", "a*exp(b*x)", "a*(1-exp(b*x)) + c", "a0 + (a1*cos(w*x) + b1*sin(w*x))"}; 0042 const char* nsl_fit_model_basic_pic_name[] = {"polynom", "power", "exponential", "inv_exponential", "fourier"}; 0043 0044 const char* nsl_fit_model_peak_name[] = {i18n("Gaussian (normal)"), i18n("Cauchy-Lorentz"), i18n("Hyperbolic secant (sech)"), i18n("Logistic (sech squared)"), 0045 i18n("Voigt profile"), i18n("Pseudo-Voigt (same width)")}; 0046 const char* nsl_fit_model_peak_equation[] = {"a/sqrt(2*pi)/s * exp(-((x-mu)/s)^2/2)", "a/pi * g/(g^2+(x-mu)^2)", "a/pi/s * sech((x-mu)/s)", 0047 "a/4/s * sech((x-mu)/2/s)**2", "a*voigt(x - mu, s, g)", "a*pseudovoigt1(x - mu, et, w)"}; // eta is already used as function 0048 const char* nsl_fit_model_peak_pic_name[] = {"gaussian", "cauchy_lorentz", "sech", "logistic", "voigt", "pseudovoigt1"}; 0049 0050 const char* nsl_fit_model_growth_name[] = {i18n("Inverse tangent"), i18n("Hyperbolic tangent"), i18n("Algebraic sigmoid"), i18n("Logistic function"), 0051 i18n("Error function (erf)"), i18n("Hill"), i18n("Gompertz"), i18n("Gudermann (gd)")}; 0052 const char* nsl_fit_model_growth_equation[] = {"a * atan((x-mu)/s)", "a * tanh((x-mu)/s)", "a * (x-mu)/s/sqrt(1+((x-mu)/s)^2)", "a/(1+exp(-k*(x-mu)))", 0053 "a/2 * erf((x-mu)/s/sqrt(2))", "a * x^n/(s^n + x^n)", "a*exp(-b*exp(-c*x))", "a * asin(tanh((x-mu)/s))"}; 0054 const char* nsl_fit_model_growth_pic_name[] = {"atan", "tanh", "alg_sigmoid", "logistic_function", "erf", "hill", "gompertz", "gd"}; 0055 0056 const char* nsl_fit_weight_type_name[] = {"No", "Instrumental (1/col^2)", "Direct (col)", "Inverse (1/col)", "Statistical (1/data)", "Statistical (Fit)", "Relative (1/data^2)", "Relative (Fit)"}; 0057 0058 /* 0059 see https://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf 0060 and https://lmfit.github.io/lmfit-py/bounds.html 0061 */ 0062 double nsl_fit_map_bound(double x, double min, double max) { 0063 if (max <= min) { 0064 printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); 0065 return DBL_MAX; 0066 } 0067 0068 /* not bounded */ 0069 if (min == -DBL_MAX && max == DBL_MAX) 0070 return x; 0071 0072 /* open bounds */ 0073 if (min == -DBL_MAX) 0074 return max + 1. - sqrt(x*x + 1.); 0075 if (max == DBL_MAX) 0076 return min - 1. + sqrt(x*x + 1.); 0077 0078 return min + (sin(x) + 1.) * (max - min)/2.; 0079 0080 /* alternative transformation for closed bounds 0081 return min + (max - min)/(1. + exp(-x)); 0082 */ 0083 } 0084 0085 /* 0086 see https://seal.web.cern.ch/seal/documents/minuit/mnusersguide.pdf 0087 and https://lmfit.github.io/lmfit-py/bounds.html 0088 */ 0089 double nsl_fit_map_unbound(double x, double min, double max) { 0090 if (max <= min) { 0091 printf("given bounds must fulfill max > min (min = %g, max = %g)! Giving up.\n", min, max); 0092 return DBL_MAX; 0093 } 0094 if (x < min || x > max) { 0095 printf("given value must be within bounds! Giving up.\n"); 0096 return -DBL_MAX; 0097 } 0098 0099 /* not bounded */ 0100 if (min == -DBL_MAX && max == DBL_MAX) 0101 return x; 0102 0103 /* open bounds */ 0104 if (min == -DBL_MAX) 0105 return sqrt(gsl_pow_2(max - x + 1.) - 1.); 0106 if (max == DBL_MAX) 0107 return sqrt(gsl_pow_2(x - min + 1.) - 1.); 0108 0109 return asin(2. * (x - min)/(max - min) - 1.); 0110 0111 /* alternative transformation for closed bounds 0112 return -log((max - x)/(x - min)); 0113 */ 0114 } 0115 0116 /********************** parameter derivatives ******************/ 0117 0118 /* basic */ 0119 double nsl_fit_model_polynomial_param_deriv(double x, int j, double weight) { 0120 return sqrt(weight)*pow(x, j); 0121 } 0122 double nsl_fit_model_power1_param_deriv(unsigned int param, double x, double a, double b, double weight) { 0123 if (param == 0) 0124 return sqrt(weight)*pow(x, b); 0125 if (param == 1) 0126 return sqrt(weight)*a*pow(x, b)*log(x); 0127 0128 return 0; 0129 } 0130 double nsl_fit_model_power2_param_deriv(unsigned int param, double x, double b, double c, double weight) { 0131 if (param == 0) 0132 return sqrt(weight); 0133 if (param == 1) 0134 return sqrt(weight)*pow(x, c); 0135 if (param == 2) 0136 return sqrt(weight)*b*pow(x, c)*log(x); 0137 0138 return 0; 0139 } 0140 double nsl_fit_model_exponentialn_param_deriv(unsigned int param, double x, double *p, double weight) { 0141 if (param % 2 == 0) 0142 return sqrt(weight)*exp(p[param+1]*x); 0143 else 0144 return sqrt(weight)*p[param-1]*x*exp(p[param]*x); 0145 } 0146 double nsl_fit_model_inverse_exponential_param_deriv(unsigned int param, double x, double a, double b, double weight) { 0147 if (param == 0) 0148 return sqrt(weight)*(1. - exp(b*x)); 0149 if (param == 1) 0150 return -sqrt(weight)*a*x*exp(b*x); 0151 if (param == 2) 0152 return sqrt(weight); 0153 0154 return 0; 0155 } 0156 double nsl_fit_model_fourier_param_deriv(unsigned int param, unsigned int degree, double x, double w, double weight) { 0157 if (param == 0) 0158 return sqrt(weight)*cos(degree*w*x); 0159 if (param == 1) 0160 return sqrt(weight)*sin(degree*w*x); 0161 0162 return 0; 0163 } 0164 0165 /* peak */ 0166 double nsl_fit_model_gaussian_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0167 double s2 = s*s, norm = sqrt(weight)/M_SQRT2/M_SQRTPI/s, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); 0168 0169 if (param == 0) 0170 return norm * efactor; 0171 if (param == 1) 0172 return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; 0173 if (param == 2) 0174 return A * norm/s2 * (x-mu) * efactor; 0175 0176 return 0; 0177 } 0178 double nsl_fit_model_lorentz_param_deriv(unsigned int param, double x, double A, double s, double t, double weight) { 0179 double norm = sqrt(weight)/M_PI, denom = s*s+(x-t)*(x-t); 0180 0181 if (param == 0) 0182 return norm * s/denom; 0183 if (param == 1) 0184 return A * norm * ((x-t)*(x-t) - s*s)/(denom*denom); 0185 if (param == 2) 0186 return A * norm * 2.*s*(x-t)/(denom*denom); 0187 0188 return 0; 0189 } 0190 double nsl_fit_model_sech_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0191 double y = (x-mu)/s, norm = sqrt(weight)/M_PI/s; 0192 0193 if (param == 0) 0194 return norm/cosh(y); 0195 if (param == 1) 0196 return A/s * norm * (y*tanh(y)-1.)/cosh(y); 0197 if (param == 2) 0198 return A/s * norm * tanh(y)/cosh(y); 0199 0200 return 0; 0201 } 0202 double nsl_fit_model_logistic_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0203 double y = (x-mu)/2./s, norm = sqrt(weight)/4./s; 0204 0205 if (param == 0) 0206 return norm/cosh(y)/cosh(y); 0207 if (param == 1) 0208 return A/s * norm * (2.*y*tanh(y)-1.)/cosh(y); 0209 if (param == 2) 0210 return A/s * norm * tanh(y)/cosh(y)/cosh(y); 0211 0212 return 0; 0213 } 0214 0215 double nsl_fit_model_voigt_param_deriv(unsigned int param, double x, double a, double mu, double s, double g, double weight) { 0216 #if !defined(_MSC_VER) 0217 if (s <= 0 || g < 0) 0218 return 0; 0219 0220 double y = x - mu, norm = a * sqrt(weight/2./M_PI)/(s*s*s); 0221 0222 double v = nsl_sf_voigt(y, s, g), im_w = nsl_sf_im_w_of_z(y); 0223 if (param == 0) 0224 return sqrt(weight) * v; 0225 if (param == 1) 0226 return a*sqrt(weight)*y/(s*s) * v - norm * g * im_w; 0227 if (param == 2) 0228 // return a*sqrt(weight)*g/M_PI/(s*s*s) + norm*sqrt(2.*M_PI)*v * (y*y+mu*mu-s*s) - norm/s*im_w*2.*g*y; 0229 return a/(s*s*s)*sqrt(weight)*(g/M_PI + v*(y*y -g*g -s*s) + im_w*2.*g*y/s); 0230 if (param == 3) 0231 return -a*sqrt(weight)/M_PI/(s*s) + norm*M_SQRT2*M_SQRTPI*s*v*g + im_w; 0232 #endif 0233 0234 return 0; 0235 } 0236 0237 double nsl_fit_model_pseudovoigt1_param_deriv(unsigned int param, double x, double a, double eta, double w, double mu, double weight) { 0238 double y = x - mu, norm = sqrt(weight); 0239 double sigma = w/sqrt(2.*M_LN2); 0240 0241 if (param == 0) 0242 return norm * nsl_sf_pseudovoigt1(y, eta, w); 0243 if (param == 1) 0244 return a * norm * (gsl_ran_cauchy_pdf(y, w) - gsl_ran_gaussian_pdf(y, sigma)); 0245 if (param == 2) 0246 return a/w * norm * (eta*(1.-2.*w*w)*gsl_ran_cauchy_pdf(y, w) + (eta-1.)*gsl_ran_gaussian_pdf(y, sigma)*(1.-2.*M_LN2*y*y/w/w)); 0247 if (param == 3) 0248 return 2.*a*y/w/w * norm * (eta*M_PI*w*gsl_pow_2(gsl_ran_cauchy_pdf(y, w)) + (1.-eta)*M_LN2*gsl_ran_gaussian_pdf(y, sigma)); 0249 0250 return 0; 0251 } 0252 0253 /* growth */ 0254 double nsl_fit_model_atan_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { 0255 double norm = sqrt(weight), y = (x-mu)/s; 0256 if (param == 0) 0257 return norm * atan(y); 0258 if (param == 1) 0259 return -A/s * norm * 1./(1+y*y); 0260 if (param == 2) 0261 return -A/s * norm * y/(1.+y*y); 0262 0263 return 0; 0264 } 0265 double nsl_fit_model_tanh_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { 0266 double norm = sqrt(weight), y = (x-mu)/s; 0267 if (param == 0) 0268 return norm * tanh(y); 0269 if (param == 1) 0270 return -A/s * norm * 1./cosh(y)/cosh(y); 0271 if (param == 2) 0272 return -A/s * norm * y/cosh(y)/cosh(y); 0273 0274 return 0; 0275 } 0276 double nsl_fit_model_algebraic_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { 0277 double norm = sqrt(weight), y = (x-mu)/s, y2 = y*y; 0278 if (param == 0) 0279 return norm * y/sqrt(1.+y2); 0280 if (param == 1) 0281 return -A/s * norm * 1./pow(1.+y2, 1.5); 0282 if (param == 2) 0283 return -A/s * norm * y/pow(1.+y2, 1.5); 0284 0285 return 0; 0286 } 0287 double nsl_fit_model_sigmoid_param_deriv(unsigned int param, double x, double A, double mu, double k, double weight) { 0288 double norm = sqrt(weight), y = k*(x-mu); 0289 if (param == 0) 0290 return norm/(1. + exp(-y)); 0291 if (param == 1) 0292 return -A*k * norm * exp(-y)/gsl_pow_2(1. + exp(-y)); 0293 if (param == 2) 0294 return A/k * norm * y*exp(-y)/gsl_pow_2(1. + exp(-y)); 0295 0296 return 0; 0297 } 0298 double nsl_fit_model_erf_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { 0299 double norm = sqrt(weight), y = (x-mu)/(sqrt(2.)*s); 0300 if (param == 0) 0301 return norm/2. * gsl_sf_erf(y); 0302 if (param == 1) 0303 return -A/M_SQRT2/M_SQRTPI/s * norm * exp(-y*y); 0304 if (param == 2) 0305 return -A/M_SQRTPI/s * norm * y*exp(-y*y); 0306 0307 return 0; 0308 } 0309 double nsl_fit_model_hill_param_deriv(unsigned int param, double x, double A, double n, double s, double weight) { 0310 double norm = sqrt(weight), y = x/s; 0311 if (param == 0) 0312 return norm * pow(y, n)/(1.+pow(y, n)); 0313 if (param == 1) 0314 return A * norm * log(y)*pow(y, n)/gsl_pow_2(1.+pow(y, n)); 0315 if (param == 2) 0316 return -A*n/s * norm * pow(y, n)/gsl_pow_2(1.+pow(y, n)); 0317 0318 return 0; 0319 } 0320 double nsl_fit_model_gompertz_param_deriv(unsigned int param, double x, double a, double b, double c, double weight) { 0321 if (param == 0) 0322 return sqrt(weight)*exp(-b*exp(-c*x)); 0323 if (param == 1) 0324 return -sqrt(weight)*a*exp(-c*x-b*exp(-c*x)); 0325 if (param == 2) 0326 return sqrt(weight)*a*b*x*exp(-c*x-b*exp(-c*x)); 0327 0328 return 0; 0329 } 0330 double nsl_fit_model_gudermann_param_deriv(unsigned int param, double x, double A, double mu, double s, double weight) { 0331 double norm = sqrt(weight), y = (x-mu)/s; 0332 if (param == 0) 0333 return -asin(tanh(y)); 0334 if (param == 1) 0335 return -A/s * norm * 1./cosh(y); 0336 if (param == 2) 0337 return -A/s * norm * y/cosh(y); 0338 0339 return 0; 0340 } 0341 0342 /* distributions */ 0343 double nsl_fit_model_gaussian_tail_param_deriv(unsigned int param, double x, double A, double s, double a, double mu, double weight) { 0344 if (x < a) 0345 return 0; 0346 double s2 = s*s, N = erfc(a/s/M_SQRT2)/2., norm = sqrt(weight)/M_SQRT2/M_SQRTPI/s/N, efactor = exp(-(x-mu)*(x-mu)/(2.*s2)); 0347 0348 if (param == 0) 0349 return norm * efactor; 0350 if (param == 1) 0351 return A * norm/(s*s2) * ((x-mu)*(x-mu) - s2) * efactor; 0352 if (param == 2) 0353 return A/norm/norm * efactor * exp(-a*a/(2.*s2)); 0354 if (param == 3) 0355 return A * norm/s2 * (x-mu) * efactor; 0356 0357 return 0; 0358 } 0359 double nsl_fit_model_exponential_param_deriv(unsigned int param, double x, double A, double l, double mu, double weight) { 0360 if (x < mu) 0361 return 0; 0362 double y = l*(x-mu), efactor = exp(-y); 0363 0364 if (param == 0) 0365 return sqrt(weight) * l * efactor; 0366 if (param == 1) 0367 return sqrt(weight) * A * (1. - y) * efactor; 0368 if (param == 2) 0369 return sqrt(weight) * A * gsl_pow_2(l) * efactor; 0370 0371 return 0; 0372 } 0373 double nsl_fit_model_laplace_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0374 double norm = sqrt(weight)/(2.*s), y = fabs((x-mu)/s), efactor = exp(-y); 0375 0376 if (param == 0) 0377 return norm * efactor; 0378 if (param == 1) 0379 return A/s*norm * (y-1.) * efactor; 0380 if (param == 2) 0381 return A/(s*s)*norm * (x-mu)/y * efactor; 0382 0383 return 0; 0384 } 0385 double nsl_fit_model_exp_pow_param_deriv(unsigned int param, double x, double a, double s, double b, double mu, double weight) { 0386 double norm = sqrt(weight)/2./s/gsl_sf_gamma(1.+1./b), y = (x-mu)/s, efactor = exp(-pow(fabs(y), b)); 0387 0388 if (param == 0) 0389 return norm * efactor; 0390 if (param == 1) 0391 return norm * a/s * efactor * (b * y * pow(fabs(1./y), 1.-b) * GSL_SIGN(y) - 1.); 0392 if (param == 2) 0393 return norm * a/b * gsl_sf_gamma(1.+1./b)/gsl_sf_gamma(1./b) * efactor * (gsl_sf_psi(1.+1./b) - gsl_pow_2(b) * pow(fabs(y), b) * log(fabs(y))); 0394 if (param == 3) 0395 return norm * a*b/s * efactor * pow(fabs(y), b-1.) * GSL_SIGN(y); 0396 0397 return 0; 0398 } 0399 double nsl_fit_model_maxwell_param_deriv(unsigned int param, double x, double a, double s, double weight) { 0400 double s2 = s*s, s3 = s*s2, norm = sqrt(weight)*M_SQRT2/M_SQRTPI/s3, x2 = x*x, efactor = exp(-x2/2./s2); 0401 0402 if (param == 0) 0403 return norm * x2 * efactor; 0404 if (param == 1) 0405 return a * norm * x2*(x2-3.*s2)/s3 * efactor; 0406 0407 return 0; 0408 } 0409 double nsl_fit_model_poisson_param_deriv(unsigned int param, double x, double A, double l, double weight) { 0410 double norm = sqrt(weight)*pow(l, x)/gsl_sf_gamma(x+1.); 0411 0412 if (param == 0) 0413 return norm * exp(-l); 0414 if (param == 1) 0415 return A/l * norm *(x-l)*exp(-l); 0416 0417 return 0; 0418 } 0419 double nsl_fit_model_lognormal_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0420 double norm = sqrt(weight)/M_SQRT2/M_SQRTPI/(x*s), y = log(x)-mu, efactor = exp(-(y/s)*(y/s)/2.); 0421 0422 if (param == 0) 0423 return norm * efactor; 0424 if (param == 1) 0425 return A * norm * (y*y - s*s) * efactor; 0426 if (param == 2) 0427 return A * norm * y/(s*s) * efactor; 0428 0429 return 0; 0430 } 0431 double nsl_fit_model_gamma_param_deriv(unsigned int param, double x, double A, double k, double t, double weight) { 0432 double factor = sqrt(weight)*pow(x, k-1.)/pow(t, k)/gsl_sf_gamma(k), efactor = exp(-x/t); 0433 0434 if (param == 0) 0435 return factor * efactor; 0436 if (param == 1) 0437 return A * factor * (log(x/t) - gsl_sf_psi(k)) * efactor; 0438 if (param == 2) 0439 return A * factor/t * (x/t-k) * efactor; 0440 0441 return 0; 0442 } 0443 double nsl_fit_model_flat_param_deriv(unsigned int param, double x, double A, double b, double a, double weight) { 0444 if (x < a || x > b) 0445 return 0; 0446 0447 if (param == 0) 0448 return sqrt(weight)/(b-a); 0449 if (param == 1) 0450 return - sqrt(weight) * A/gsl_pow_2(a-b); 0451 if (param == 2) 0452 return sqrt(weight) * A/gsl_pow_2(a-b); 0453 0454 return 0; 0455 } 0456 double nsl_fit_model_rayleigh_param_deriv(unsigned int param, double x, double A, double s, double weight) { 0457 double y=x/s, norm = sqrt(weight)*y/s, efactor = exp(-y*y/2.); 0458 0459 if (param == 0) 0460 return norm * efactor; 0461 if (param == 1) 0462 return A*y/(s*s) * (y*y-2.)*efactor; 0463 0464 return 0; 0465 } 0466 double nsl_fit_model_rayleigh_tail_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0467 double norm = sqrt(weight)*x/(s*s), y = (mu*mu - x*x)/2./(s*s); 0468 0469 if (param == 0) 0470 return norm * exp(y); 0471 if (param == 1) 0472 return -2. * A * norm/s * (1. + y) * exp(y); 0473 if (param == 2) 0474 return A * mu * norm/(s*s) * exp(y); 0475 0476 return 0; 0477 } 0478 double nsl_fit_model_levy_param_deriv(unsigned int param, double x, double A, double g, double mu, double weight) { 0479 double y=x-mu, norm = sqrt(weight)*sqrt(g/(2.*M_PI))/pow(y, 1.5), efactor = exp(-g/2./y); 0480 0481 if (param == 0) 0482 return norm * efactor; 0483 if (param == 1) 0484 return A/2.*norm/g/y * (y - g) * efactor; 0485 if (param == 2) 0486 return A/2.*norm/y/y * (3.*y - g) * efactor; 0487 0488 return 0; 0489 } 0490 double nsl_fit_model_landau_param_deriv(unsigned int param, double x, double weight) { 0491 if (param == 0) 0492 return sqrt(weight) * gsl_ran_landau_pdf(x); 0493 0494 return 0; 0495 } 0496 double nsl_fit_model_chi_square_param_deriv(unsigned int param, double x, double A, double n, double weight) { 0497 double y=n/2., norm = sqrt(weight)*pow(x, y-1.)/pow(2., y)/gsl_sf_gamma(y), efactor = exp(-x/2.); 0498 0499 if (param == 0) 0500 return norm * efactor; 0501 if (param == 1) 0502 return A/2. * norm * (log(x/2.) - gsl_sf_psi(y)) * efactor; 0503 0504 return 0; 0505 } 0506 double nsl_fit_model_students_t_param_deriv(unsigned int param, double x, double A, double n, double weight) { 0507 if (param == 0) 0508 return sqrt(weight) * gsl_ran_tdist_pdf(x, n); 0509 if (param == 1) 0510 return sqrt(weight) * A * gsl_sf_gamma((n+1.)/2.)/2./pow(n, 1.5)/M_SQRTPI/gsl_sf_gamma(n/2.) * pow(1.+x*x/n, - (n+3.)/2.) 0511 * (x*x - 1. - (n+x*x)*log1p(x*x/n) + (n+x*x)*(gsl_sf_psi((n+1.)/2.) - gsl_sf_psi(n/2.)) ) ; 0512 0513 return 0; 0514 } 0515 double nsl_fit_model_fdist_param_deriv(unsigned int param, double x, double A, double n1, double n2, double weight) { 0516 double norm = sqrt(weight) * gsl_sf_gamma((n1+n2)/2.)/gsl_sf_gamma(n1/2.)/gsl_sf_gamma(n2/2.) * pow(n1, n1/2.) * pow(n2, n2/2.) * pow(x, n1/2.-1.); 0517 double y = n2+n1*x; 0518 0519 if (param == 0) 0520 return sqrt(weight) * gsl_ran_fdist_pdf(x, n1, n2); 0521 if (param == 1) 0522 return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n2*(1.-x) + y*(log(n1) + log(x) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n1/2.))); 0523 if (param == 2) 0524 return A/2. * norm * pow(y, -(n1+n2+2.)/2.) * (n1*(x-1.) + y*(log(n2) - log(y) + gsl_sf_psi((n1+n2)/2.) - gsl_sf_psi(n2/2.))); 0525 0526 return 0; 0527 } 0528 double nsl_fit_model_beta_param_deriv(unsigned int param, double x, double A, double a, double b, double weight) { 0529 double norm = sqrt(weight) * A * gsl_sf_gamma(a+b)/gsl_sf_gamma(a)/gsl_sf_gamma(b) * pow(x, a-1.) * pow(1.-x, b-1.); 0530 0531 if (param == 0) 0532 return sqrt(weight) * gsl_ran_beta_pdf(x, a, b); 0533 if (param == 1) 0534 return norm * (log(x) - gsl_sf_psi(a) + gsl_sf_psi(a+b)); 0535 if (param == 2) 0536 return norm * (log(1.-x) - gsl_sf_psi(b) + gsl_sf_psi(a+b)); 0537 0538 return 0; 0539 } 0540 double nsl_fit_model_pareto_param_deriv(unsigned int param, double x, double A, double a, double b, double weight) { 0541 if (x < b) 0542 return 0; 0543 0544 double norm = sqrt(weight) * A; 0545 if (param == 0) 0546 return sqrt(weight) * gsl_ran_pareto_pdf(x, a, b); 0547 if (param == 1) 0548 return norm * pow(b/x, a) * (1. + a * log(b/x))/x; 0549 if (param == 2) 0550 return norm * a*a * pow(b/x, a-1.)/x/x; 0551 0552 return 0; 0553 } 0554 double nsl_fit_model_weibull_param_deriv(unsigned int param, double x, double A, double k, double l, double mu, double weight) { 0555 double y = (x-mu)/l, z = pow(y, k), efactor = exp(-z); 0556 0557 if (param == 0) 0558 return sqrt(weight) * k/l * z/y * efactor; 0559 if (param == 1) 0560 return sqrt(weight) * A/l * z/y*(k*log(y)*(1.-z) + 1.) * efactor; 0561 if (param == 2) 0562 return sqrt(weight) * A*k*k/l/l * z/y*(z-1.) * efactor; 0563 if (param == 3) 0564 return sqrt(weight) * A*k/l/l * z/y/y*(k*z + 1. - k) * efactor; 0565 0566 return 0; 0567 } 0568 double nsl_fit_model_frechet_param_deriv(unsigned int param, double x, double A, double g, double s, double mu, double weight) { 0569 double y = (x-mu)/s, efactor = exp(-pow(y, -g)); 0570 0571 if (param == 0) 0572 return g * sqrt(weight)/s * pow(y, -g-1.) * efactor; 0573 if (param == 1) 0574 return sqrt(weight) * A/s * pow(y, -2.*g-1.) * (g*log(y)*(1.-pow(y, g))+pow(y, g)) * efactor; 0575 if (param == 2) 0576 return A * sqrt(weight) * gsl_pow_2(g/s)*pow(y, -2.*g-1.) * (pow(y, g)-1.) * efactor; 0577 if (param == 3) 0578 return A * sqrt(weight) * g/(s*s)*pow(y, -g-2.) * (g+1.-g*pow(y, -g)) * efactor; 0579 0580 return 0; 0581 } 0582 double nsl_fit_model_gumbel1_param_deriv(unsigned int param, double x, double A, double s, double mu, double b, double weight) { 0583 double norm = sqrt(weight)/s, y = (x-mu)/s, efactor = exp(-y - b*exp(-y)); 0584 0585 if (param == 0) 0586 return norm * efactor; 0587 if (param == 1) 0588 return A/s * norm * (y - 1. - b*exp(-y)) * efactor; 0589 if (param == 2) 0590 return A/s * norm * (1. - b*exp(-y)) * efactor; 0591 if (param == 3) 0592 return -A * norm * exp(-y) * efactor; 0593 0594 return 0; 0595 } 0596 double nsl_fit_model_gumbel2_param_deriv(unsigned int param, double x, double A, double a, double b, double mu, double weight) { 0597 double y = x - mu, norm = A * sqrt(weight) * exp(-b * pow(y, -a)); 0598 0599 if (param == 0) 0600 return sqrt(weight) * gsl_ran_gumbel2_pdf(y, a, b); 0601 if (param == 1) 0602 return norm * b * pow(y, -1. -2.*a) * (pow(y, a) -a*(pow(y, a)-b)*log(y)); 0603 if (param == 2) 0604 return norm * a * pow(y, -1. -2.*a) * (pow(y, a) - b); 0605 if (param == 3) 0606 return norm * a * b * pow(y, -2.*(a + 1.)) * ((1. + a)*pow(y, a) - a*b); 0607 0608 return 0; 0609 } 0610 double nsl_fit_model_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { 0611 if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) 0612 return 0; 0613 k = round(k); 0614 n = round(n); 0615 0616 double norm = sqrt(weight) * gsl_sf_fact((unsigned int)n)/gsl_sf_fact((unsigned int)(n-k))/gsl_sf_fact((unsigned int)(k)); 0617 if (param == 0) 0618 return sqrt(weight) * gsl_ran_binomial_pdf((unsigned int)k, p, (unsigned int)n); 0619 if (param == 1) 0620 return A * norm * pow(p, k-1.) * pow(1.-p, n-k-1.) * (k-n*p); 0621 if (param == 2) 0622 return A * norm * pow(p, k) * pow(1.-p, n-k) * (log(1.-p) + gsl_sf_psi(n+1.) - gsl_sf_psi(n-k+1.)); 0623 0624 return 0; 0625 } 0626 double nsl_fit_model_negative_binomial_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { 0627 if (k < 0 || k > n || n < 0 || p < 0 || p > 1.) 0628 return 0; 0629 0630 double norm = A * sqrt(weight) * gsl_sf_gamma(n+k)/gsl_sf_gamma(k+1.)/gsl_sf_gamma(n); 0631 if (param == 0) 0632 return sqrt(weight) * gsl_ran_negative_binomial_pdf((unsigned int)k, p, n); 0633 if (param == 1) 0634 return - norm * pow(p, n-1.) * pow(1.-p, k-1.) * (n*(p-1.) + k*p); 0635 if (param == 2) 0636 return norm * pow(p, n) * pow(1.-p, k) * (log(p) - gsl_sf_psi(n) + gsl_sf_psi(n+k)); 0637 0638 return 0; 0639 } 0640 double nsl_fit_model_pascal_param_deriv(unsigned int param, double k, double A, double p, double n, double weight) { 0641 return nsl_fit_model_negative_binomial_param_deriv(param, k, A, p, round(n), weight); 0642 } 0643 double nsl_fit_model_geometric_param_deriv(unsigned int param, double k, double A, double p, double weight) { 0644 if (param == 0) 0645 return sqrt(weight) * gsl_ran_geometric_pdf((unsigned int)k, p); 0646 if (param == 1) 0647 return A * sqrt(weight) * pow(1.-p, k-2.) * (1.-k*p); 0648 0649 return 0; 0650 } 0651 double nsl_fit_model_hypergeometric_param_deriv(unsigned int param, double k, double A, double n1, double n2, double t, double weight) { 0652 if (t > n1 + n2) 0653 return 0; 0654 0655 double norm = sqrt(weight) * gsl_ran_hypergeometric_pdf((unsigned int)k, (unsigned int)n1, (unsigned int)n2, (unsigned int)t); 0656 if (param == 0) 0657 return norm; 0658 if (param == 1) 0659 return A * norm * (gsl_sf_psi(n1+1.) - gsl_sf_psi(n1-k+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); 0660 if (param == 2) 0661 return A * norm * (gsl_sf_psi(n2+1.) - gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2+1.) + gsl_sf_psi(n1+n2-t+1.)); 0662 if (param == 3) 0663 return A * norm * (gsl_sf_psi(n2+k-t+1.) - gsl_sf_psi(n1+n2-t+1.) - gsl_sf_psi(t-k+1.) + gsl_sf_psi(t+1.)); 0664 0665 return 0; 0666 } 0667 double nsl_fit_model_logarithmic_param_deriv(unsigned int param, double k, double A, double p, double weight) { 0668 if (param == 0) 0669 return sqrt(weight) * gsl_ran_logarithmic_pdf((unsigned int)k, p); 0670 if (param == 1) 0671 return A * sqrt(weight) * pow(1.-p, k-2.) * (1.-k*p); 0672 0673 return 0; 0674 } 0675 double nsl_fit_model_sech_dist_param_deriv(unsigned int param, double x, double A, double s, double mu, double weight) { 0676 double norm = sqrt(weight)/2./s, y = M_PI_2*(x-mu)/s; 0677 0678 if (param == 0) 0679 return norm * 1./cosh(y); 0680 if (param == 1) 0681 return -A/s * norm * (y*tanh(y)+1.)/cosh(y); 0682 if (param == 2) 0683 return A*M_PI_2/s * norm * tanh(y)/cosh(y); 0684 0685 return 0; 0686 }