luisfabib
commented
4 years ago Hey, thanks for the inputs. I will try to go through the issues in order to answer or complement them:
Fixed variables and boundaries
Individual parameters can be fixed in DeerLab by defining a new function handle for the model. For example: take a Gauss model dd_gauss(r,[rmean fwhm]) which you would normaslly fit via
paramfit = fitparamodel(V,@dd_gauss,r,K);but now you want to fit rmean with a fixed fwhm, then one needs to use
fwhm = 0.54; %nm
model = @(r,rmean) dd_gauss(r,[rmean fwhm]);
par0 = 5; %nm
paramfit = fitparamodel(V,@dd_gauss,r,K,par0,'lower',1,'upper',20);since the fwhm is defined before the construction of the function handle model, that value will be fixed inside the new model. Then fitparamodel can be called, albeit this time the start values par0 as well as upper/lower boundaries need to be specified.
There is currently no way to fix a parameter of a built-in model while keeping the boundaries/start values information of the other parameters.
Global fit
Similarly, global fitting in DeerLab allows to fit models with fixed parameters, global parameters as well as local parameters.
The following script is an example of fitting two different signal to two different bimodal Gauss distance distributions, whose mean distances are equal for the two signals but their amplitudes differ, i.e. the mean values are global parameters, the amplitudes are local parameters and the widths are fixed to a known value.
This kind of processing requires full scripting as (currently) there are no wrapper functions around this.
clc,clear,clf
rng(1)
%preparations
r = linspace(2,6,300);
t1 = linspace(0,4,200);
t2 = linspace(0,6,100);
%parameters
rmeanA = 3.45;
rmeanB = 5.05;
fwhmA = 0.5;
fwhmB = 0.3;
ampA1 = 0.8;
ampA2 = 0.2;
%two different distributions for two different samples
P1 = dd_gauss2(r,[rmeanA fwhmA ampA1 rmeanB fwhmB]);
P2 = dd_gauss2(r,[rmeanA fwhmA ampA2 rmeanB fwhmB]);
%...and two different signals for two different samples
V1 = dipolarsignal(t1,r,P1,'noiselevel',0.01);
V2 = dipolarsignal(t2,r,P2,'noiselevel',0.01);
%dipolar kernels
K1 = dipolarkernel(t1,r);
K2 = dipolarkernel(t2,r);
%parameters to fit [rmeanA rmeanB ampA_1 ampA_2]
par0 = [2 2 0.5 0.5];
lower = [1 1 0 0];
upper = [20 20 1 1];
%Collect data for global fit
Ks = {K1,K2};
Vs = {V1,V2};
ts = {t1,t2};
%Fit
model = @(t,par,idx) mymodel(t,par,idx,r,Ks);
parfit = fitparamodel(Vs,model,ts,par0,'lower',lower,'upper',upper,'multistart',50);
%Get results
[Vfit1,Pfit1] = mymodel(t1,parfit,1,r,Ks);
[Vfit2,Pfit2] = mymodel(t2,parfit,2,r,Ks);
subplot(221)
plot(r,P1,r,Pfit1)
subplot(222)
plot(t1,V1,t1,Vfit1)
subplot(223)
plot(r,P2,r,Pfit2)
subplot(224)
plot(t2,V2,t2,Vfit2)
function [Vfit,Pfit] = mymodel(t,par,idx,r,Ks)
%Global parameters
rmeanA = par(1);
rmeanB = par(2);
%Fixed parameters
fwhmA = 0.5;
fwhmB = 0.3;
%Local parameters
ampA1 = par(3);
ampA2 = par(4);
if idx == 1
%Use ampA1 for signal Vs{1}
AMP = ampA1;
else
%Use ampA2 for signal Vs{2}
AMP = ampA2;
end
%Generate signal-specific distribution (mix global+local parameters)
Pfit = dd_gauss2(r,[rmeanA fwhmA AMP rmeanB fwhmB]);
%Generate corresponding signal (either V{1} or V{2})
Vfit = Ks{idx}*Pfit;
end
Bilevel optimization
Yes, nested LSQ optimization of full models (we use this term now instead of bilevel optimization) is now way more accessible from the new fitsignal function. fitsignal is a pure wrapper around fitparamodel and fitregmodel so therefore, all workflows that were done before this function are still available if scripted.
The main advantages of fitsignal over a custom-implementation is that this wrapper uses several tricks to largely boost computation times (up to factor 6 in speed for some cases that we have seen) and handles the confidence interval computation for this kind of nested problems which is not trivial.
However, as you said in 0.9.0 there is no way to fix a parameter in fitsignal while fitting the rest. This is one of the many interface issues that are still to be tackled in future releases of DeerLab.
Questions
Am I correct, that bilevel optimization was replaced with fitsignal?
As explained above, it has not been "removed" since it is still scriptable, but fitsignal functions as a wrapper for easier scripting and speed.
Is there a feature for fixed variables or boundaries in your code I missed? Could this I would be willing to provide an example for such analysis, if this is of interest for the repo.
See examples above. If you meant something else please let me know ;)
Is it possible to include the background model, but not the modulation depth in the kernel to use with fitparamodel?
I understand that you want to simulate/fit a system with 100% modulation depth that still has a background function suppressing the modulations? This can be accomplished by just specifying a modulation depth lambda=1, i.e
lam = 1;
B = bg_exp(t,k);
K = dipolarkernel(t,r,lam,B)Do you have better workarounds for this problems/restrictions in mind?
The future releases of DeerLab will primarily focus on the user-program interface to try to simplify different workflows and processing types.
Hopefully this helps, otherwise please let us know ;)
laenan8466
Processing with nested function as in 
stestoll
Here is to note, that the result differs when perfoming the analysis a second time. This could be due to the random MultiStart.
Hey, sadly I havn't followed all changes in the last month. When revisting an 'old' project I came up with some open questions about the design for parametric fitting. I try to group the questions together. I also add some points I noticed, that could be interesting for others and I would like to emphasize because they are genius (tremendous work @luisfabib and all others).
Parametric model fit
Things I noticed/proposed features
Fixed variables and boundaries
With
fitparamodelit is possible to change the inital guess, but it is not possible to fix a parameter individually. This is of interest, if a parameter can be determined with higher accuracy in a single experiment. Possible workarounds could be bilevel fit (see below) or a custom dd-model with a fixed parameter/specialized boundaries.Background for parametric fit
When performing a parametric fit (
fitparamodel) of a mixed model with the example provided it is not possible to include a background correction, if the modulation depth is unknown. The reason for this is, that thedipolarkernelcannot be generated with a background model but unknown modulation depth. See documentation of dipolarkernel. Workaround could be using bilevel fit (see below) orfitsignal(with the restrictions listed below).Global fit (proposal for feature addition)
It could be of interest to not only fit several traces, resulting in one distribution, but fit only some parameters (in a parametric model e.g. the mean position but not the width) and determine others individually.
Bilevel optimization
In a version 0.8 bilevel optimization was existing as an example and was removed in v0.9. I guess the reason that it was removed is that
fitsignalis now doing this job. The way bilevel optimization worked could be used for fixed variables/boundaries or global fit of several (different) samples.Multi-component parametric fit
The possiblities introduced with
fitmultimodelare really great. This could be helpful if having a series of macromolecules with increased size and it is unclear what is included (e.g. polymerization).Questions to @luisfabib
fitsignal?fitparamodel?Maybe it will be helpful for you to devide some feature requests in seperate issues. If you are interested in having s.th. of that I proposed, I can add it in a seperate issue. Just give me a hint. ;-)