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