fitting is not working
Ältere Kommentare anzeigen
%% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess = min(R2);
tau_guess = 0.963;
p0 = [R2m_guess, tau_guess];
%% Bounds (lower, upper)
lb = [0, 0, 1e-4];
ub = [1, 1000, 500];
%% Fit using lsqcurvefit (Optimization Toolbox)
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);
R2m_fit = params_fit(1);
tau_fit = params_fit(2);
fprintf('\n--- Fit results ---\n');
fprintf('R2m = %.5f\n', R2m_fit);
fprintf('tau = %.5f\n', tau_fit);
fprintf('Residual norm = %.5f\n', resnorm);
%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;
1 Kommentar
Walter Roberson
vor etwa eine Stunde
Your upper and lower bounds are three elements each, but you are only searching over 2 parameters.
Antworten (3)
As far as I can see, the function f(x) = p*(1-2*p./x.*tanh(x/(2*p))) is monotonically increasing (the first parameter params(1) is not important in the analysis since it's only a translation of the curve in the upward or downward direction). Thus you will never get the bellshape of your experimental data for a fit result. Try plotting the function from above for different values of p, and you will see what I mean.
Thus in short: the model function is not suited to reflect your experimental data.
x = 1:14;
p = [0.003;0.03;0.3;3;30;300];
f = p.*(1-2*p./x.*tanh(x./(2*p)));
plot(x,f)
9 Kommentare
John D'Errico
vor 10 Minuten
Bearbeitet: John D'Errico
vor 9 Minuten
It appears you (both) may be using the wrong model form. That is, when I look online for that model, I see a different form.
In there (aswell in all other places I saw that model), the form appears to have an extra 1/x term in it. Something more like this (again, I've dropped extraneous constants so we can see the fundamental shape.)
f = @(x,p) 1./x.*(1-2*p./x.*tanh(x./(2*p)));
figure
fplot(@(x) f(x,1),[0,50],'r')
hold on
fplot(@(x) f(x,2),[0,50],'g')
fplot(@(x) f(x,3),[0,50],'b')
Anyway, this now has a skewed bell shape to it. At the same time, the shape is still not terribly consistent with the data provided, so I doubt it will help greatly.
Ehtisham
vor 34 Minuten
Replace
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
by the modified model function.
Ehtisham
vor etwa 9 Stunden
Bearbeitet: Walter Roberson
vor etwa 8 Stunden
% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess = min(R2);
tau_guess = 0.963;
p0 = [R2m_guess, tau_guess];
%% Bounds (lower, upper)
lb = [0, 0, 1e-4];
ub = [1, 100, 50];
%% Fit using lsqcurvefit
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);
R2m_fit = params_fit(1);
tau_fit = params_fit(2);
fprintf('\n--- Fit results ---\n');
fprintf('R2m = %.5f\n', R2m_fit);
fprintf('tau = %.5f\n', tau_fit);
fprintf('Residual norm = %.5f\n', resnorm);
%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;
i replace that but giving the same plot
Torsten
vor 7 Minuten
I don't see changes in
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
compared to your previous code.
And you still use three-element vectors for lower and upper bounds on the parameters instead of two-element vectors.
Ehtisham
vor etwa 8 Stunden
Bearbeitet: Walter Roberson
vor etwa 8 Stunden
See i change the funcztion but still getting the same result
% Luz-Meiboom fit of R2 vs pH
clear; clc; close all;
%% Data
pH = [1.88, 2.45, 3.39, 3.87, 4.53, 4.68, 5.12, 5.25, 5.55, 5.82, ...
6.01, 6.27, 6.66, 6.98, 7.75, 8.25, 8.96, 9.23];
R2 = [0.01137, 0.00973, 0.00803, 0.01824, 0.0547, 0.0712, 0.10246, ...
0.12065, 0.12992, 0.22537, 0.18405, 0.05077, 0.0203, 0.01554, ...
0.02122, 0.0371, 0.03736, 0.00611];
pH = pH(:);
R2 = R2(:);
%% Luz-Meiboom model function
% params(1) = R2m
% params(2) = tau
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
%% Initial parameter guesses
R2m_guess = min(R2);
tau_guess = 1;
p0 = [R2m_guess, tau_guess];
%% Bounds (lower, upper)
lb = [0, 1e-4];
ub = [1, 50];
%% Fit using lsqcurvefit
opts = optimoptions('lsqcurvefit', 'Display', 'iter');
[params_fit, resnorm, residual, exitflag] = lsqcurvefit(luzMeiboomFun, p0, pH, R2, lb, ub, opts);
R2m_fit = params_fit(1);
tau_fit = params_fit(2);
fprintf('\n--- Fit results ---\n');
fprintf('R2m = %.5f\n', R2m_fit);
fprintf('tau = %.5f\n', tau_fit);
fprintf('Residual norm = %.5f\n', resnorm);
%% Generate smooth fitted curve for plotting
pH_fine = linspace(min(pH), max(pH), 300);
R2_fit_curve = luzMeiboomFun(params_fit, pH_fine);
%% Plot: pH vs R2
figure('Color','w');
plot(pH, R2, 'ko', 'MarkerFaceColor','r', 'MarkerSize',7, 'DisplayName','Data');
hold on;
plot(pH_fine, R2_fit_curve, 'b-', 'LineWidth', 1.8, 'DisplayName','Luz-Meiboom fit');
xlabel('pH', 'FontSize', 12);
ylabel('R_2 (s^{-1})', 'FontSize', 12);
title('R_2 vs pH — Luz-Meiboom Fit', 'FontSize', 13);
legend('Location','best');
grid on;
box on;
hold off;
Walter Roberson
vor etwa 8 Stunden
Your code
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
does not contain a 1./x term.
Ehtisham
vor etwa 7 Stunden
I didn't say to use
luzMeiboomFun = @(params, x) params(1) + params(2).* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
- this is what you implemented right from the beginning and which will for sure be inadequate for your experimental data. I said to replace this function by the correct model function (which we don't know how it looks). As far as I understand, @John D'Errico made the suggestion to use
luzMeiboomFun = @(params, x) params(1) + params(2)./x.* ...
(1 - (2.*params(2)./x).*tanh(x./(2.*params(2))));
instead, but his last comment makes me guess that it won't improve things greatly.
So search for the adequate Luz Meiboom function for your application.
Hi @Ehtisham,
Just want to close the loop on this one. Walter was right that the bounds had an extra entry that wasn't even being used (3 bounds for 2 parameters), and Torsten's point was the real issue: the original Luz-Meiboom function only goes up as x increases, so it could never match a curve that rises then falls the way this R2 vs pH data does — no amount of retuning the bounds or starting guesses was ever going to fix that, because the shape of the equation itself was wrong for this job. The fix was to stop plugging pH directly into that equation and instead let pH control the exchange rate (kex), using the standard acid/base-catalyzed exchange formula, so pH shapes the curve indirectly instead of being treated like a CPMG pulse spacing. After fitting the parameters in log-scale (they ranged from single digits to billions, which was making the optimizer wander uselessly) and running it from a bunch of different starting points to rule out getting stuck, the model landed on a real bell-shaped curve peaking right around pH 5.7-5.8, matching the data with R^2 = 0.85 and physically sensible numbers (the fitted CPMG spacing came out to a normal 60-170 microseconds). The one spot it still doesn't nail is the single sharp spike at pH 5.82, which a smooth physical curve like this can't fully reproduce — worth checking whether that point (and pH 6.01 next to it) were repeated measurements or a one-off, since that's a data question rather than a fitting one at this point. Attaching a PDF with the final fitted plot, the script output, and the full script for reference.
Umar
vor etwa 3 Stunden
0 Stimmen
P.S. — meant to add this context earlier. Zooming out from the fitting mechanics: that bell shape in R2 vs pH (low at both ends, spiking around pH 5.8) is a classic sign of a proton exchange process passing through "coalescence" as pH changes. Some proton in the sample — could be on a side chain like histidine, or an amide/hydroxyl proton — is exchanging with the solvent, and how fast that exchange happens depends on pH, since both H+ and OH- catalyze it (faster exchange at low pH from acid catalysis, faster again at high pH from base catalysis, slower in between). R2 only gets strongly broadened by this kind of exchange when the exchange rate lines up closely with the frequency difference between the two states the proton is hopping between — too fast or too slow relative to that, and the broadening drops off. That's exactly why R2 climbs, peaks, then falls instead of just rising or falling monotonically with pH.
So what's really being pulled out of this data isn't just "a fit" — it's the underlying exchange rate constants (acid-catalyzed and base-catalyzed rate constants, or equivalently the exchange rate at any given pH) and, likely, something like a pKa: the pH where exchange is fastest, which usually reflects the protonation behavior of whatever residue or group is being studied. That's a fairly standard target for NMR relaxation data like this. The original code just had the wrong equation for it — built for how R2 depends on CPMG pulse spacing, not on pH — so it could never reproduce the bell shape no matter how the fit was tuned. That's the real motivation behind switching to the acid/base-catalyzed exchange-rate model above.
vor etwa 3 Stunden
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
