Curve fitting is an integral part of quantitative analysis. Usually the experimental data is fitted to a predictive analytical model. Here an example code has been explained to fit CPMG relaxation dispersion NMR data using Richard Craver model. The codes are written in both matlab as well as C language. In this file please set the parameterMatlab: In the attached downloadable matlab zip file the following files are present: 1) Curve_fit.m param_search=1; % If exhaustive search for parameters is needed; takes long time but gives excellent fit or param_search=0; % for quick fit,sometimes fitting may fail criteria_res_norm=1e-2; % Minimum squared sum of residuals, lower the value more stringent the fitting will be and yields % excellent results function z=model_func(a,x); function z=model_func(a,x); [row,col]=size(x); for i=1:row z(i,1)=a(1)*x(i)+a(2); end note that in the first definition of z=a(1).*x+a(2); we have used .* for element wise multiplication of vector elements. In the attachment file the "Richard Craver" model was defined. 3) Data: This file contains the x and y values as first and second column. Example content of Data file is: 0 21.8783973163 1 32.6781264109 2 51.4573515575 3 90.2250441136 4 115.6549521303 5 188.5256850821 6 218.6186036507 7 268.4985192822 8 294.7800862444 9 314.6647993941 10 287.9013078371 11 277.8269737847 12 327.831676702 13 297.2861193764 14 312.7020443178 The first column is X value, the second column is Y value. 4) Limits: This file contains the lower and upper limit of the fitting parameter. The number of rows corresponds to the number of fit parameters in the model function. For the above example of linear line fit, we have defined the model function as z=a(1).*x+a(2). The first and the second parameter are a(1) and a(2) which corresponds to 'm' and 'c' in the analytical expression. So the content of the limit files are as follows, -10000 10000 0 5 -20000 20000 0 4 The first row -10000 10000 0 5 corresponds to paramater 'm'. The first column -10000 is the lower limit, 10000 is the upper limit within which the algorithm will search for the best value of 'm'. the third column '0' specifies that the fit value is either positive or negative. '1' is for only poisitive, '-1' is for only negative value. The fourth column 5 specifies that the algorithm should split the interval -10000 to 10000 in to 5 geometric progression based partition (say, -10000, -100, 0, 100, 10000). lesser the magnitude of this value the faster the program completes, but may fail to converge to best solution. Larger the value for this entry longer the program runs but ensures good results. 5) Results_out: This is the name of the results of the fitting with standard error and the back calculated value of Y, along with the confidence interval. Example fit result for 4 CPMG parameters are as follows: The fitted parameters with standard error and confidence interval at 95% level (alpha=5%) __________________________________________________________________________ -10.000353 +/- 0.011333 [-10.023791 to -9.973288] -2591.00502 +/- 1.836048 [-2595.399443 to -2587.217504] 0.502156 +/- 0.005422 [0.487923 to 0.512086] 12.999901 +/- 0.010593 [12.974459 to 13.021662] The fitted data and its confidence interval at 95% level (alpha=5%) _________________________________________________________ x_expt y_exp y_pred [y_ci_lower y_ci_upper] ------------------------------------------------------------------------------ 20.000000 7.870000 7.871026 [7.857646 to 7.884407] 32.800000 7.800000 7.800339 [7.786799 to 7.813879] 65.600000 7.780000 7.778339 [7.763084 to 7.793593] 131.200000 7.010000 7.011901 [7.005356 to 7.018447] 262.400000 6.740000 6.736953 [6.731555 to 6.742352] 328.000000 5.670000 5.671648 [5.667968 to 5.675328] 393.600000 4.970000 4.969932 [4.967188 to 4.972676] 459.200000 4.500000 4.500416 [4.497963 to 4.502869] 524.800000 4.180000 4.175870 [4.173345 to 4.178394] 590.400000 3.940000 3.944011 [3.941288 to 3.946734] 656.000000 3.770000 3.773382 [3.770445 to 3.776319] 721.600000 3.640000 3.644522 [3.641391 to 3.647654] 787.200000 3.550000 3.545007 [3.541709 to 3.548306] 852.800000 3.470000 3.466652 [3.463213 to 3.470090] C Language: If you want to perform Curve fit using a c program, you can try Curve_fit_C.zip (see the bottom of this page). This program requires GSL scientific library (>=1.16) for its compilation. Please look in to the README file for how to compile this program. The attached file contains the source file: LVM_primer.c and LVM_algorithm_modified.c Please change the parameters: #define Trial_max 100 // Increase for deeper search but will take more time #define Geometric_search 2 // 0) Quick random search uses CHI_SQUARE_LIMIT & Trial_max // 1) Complete random search tries Trial_max times and chooses the best of lowest chi-square // 2) Completely organized; chooses initial parameters based on geometric progression and outputs best chi-squared values. #define CHI_SQUARE_LIMIT 0.00001 // Lower the chi-square the better the fitting will be but may take long time. If fitting // does'nt occur give a higher value like 0.1 #define Iter_max 10000 // Higher this value better the fit will be but it may take long time. Modify the model function func_calc int func_calc(double *param, void *calc_data) { int n = ((Data_nongsl *)calc_data)->n; int p = ((Data_nongsl *)calc_data)->p; double *x_data = ((Data_nongsl *)calc_data)->x_data; int i,j; for(j=0;j<p;j++){if(param[j]==0){param[j]=pow(10,-10);}} for(i=0;i<n;i++) { ((Data_nongsl *)calc_data)->y_calc_data[i]= param[0]+param[1]*exp(-1*param[2]*x_data[i]); } } Change only the line highlighted in blue. In the above case the model function is an exponential: Y=A+B*exp(-1*C*X); for a linear fit the above line would be: ((Data_nongsl *)calc_data)->y_calc_data[i]= param[0]+param[1]*x_data[i]; note that param[0] corresponds to 'c' and param[1] corresponds to 'm'. If the model includes many fit parameters, use param[0] to param[n], where n is the index for the last parameter. Accordingly, the number of rows in the Limit file should match with the number of parameters present in the model (see below for the preparation of input files). Compile the program LVM_primer.c as follows (no need to compile LVM_algorithm_modified.c as this file will be compiled automatically when LVM_primer.c is compiled), gcc -o LVM_primer LVM_primer.c -lm /home/Software/GSL/lib/libgsl.a /home/Software/GSL/lib/libgslcblas.a Run the program as ./LVM_primer This program requires two input files: Data, Limit. Prepare the input files are explained above for matlab curve fit program or see the example files present in the attached file. The model considered here is Richard Craver model for CPMG relaxation dispersion. Fit_results and Raw_results are the output file. Raw_results contains all the fits including the ones that are bad. Fit_results is the best of all the fits present in Raw_results. (Please click the download arrow to download the files correctly  ) |
