Demonstration of two-parameter Rouse regression

twoparam_nonlinear_rouse.m

% load some example data
sampleDepth = [0.9500; 0.9500; 0.9500; 0.8500; 0.8500; 0.8500; 0.7500; 0.7500; 0.7500; 0.5000; 0.5000; 0.5000; 0.1000; 0.1000; 0.1000];
sampleZ = [0.0850; 0.0850; 0.0850; 0.2550; 0.2550; 0.2550; 0.4250; 0.4250; 0.4250; 0.8500; 0.8500; 0.8500; 1.5300; 1.5300; 1.5300];
conc = [1.5897; 1.8884; 2.3833; 1.0620; 0.8321; 0.8760; 0.4087; 0.5752; 0.4819; 0.4908; 0.6069; 0.5313; 0.2321; 0.2569; 0.2492];


% organize into a model structure
modelParams.flowDepth = 1.7;
modelParams.b = 0.0850;
modelParams.Hbb = (modelParams.flowDepth - modelParams.b) / modelParams.b; % this is ((H-b)/b), a constant


% datatable for regression
samplesTable = array2table([sampleDepth, sampleZ, conc], 'VariableNames', {'sampleDepth', 'sampleZ', 'conc'});
  

% regression formula has four vars (conc, cb, sampleZ, and Rou. others are hardcoded into formula)
modelParams.formula = ['conc ~ cb * ( ((', num2str(modelParams.flowDepth), '-sampleZ)/sampleZ) / ', num2str(modelParams.Hbb), ' )^Rou'];


% samplesTable supplies the conc and sampleZ, so guess is needed for other params
modelParams.guess = [nanmean(samplesTable.conc(samplesTable.sampleDepth > 0.9)) 0.5]; % guesses for cb and Rou
modelParams.model = fitnlm(samplesTable, modelParams.formula, modelParams.guess);


% pull regressed variables outs
modelParams.Rou = modelParams.model.Coefficients.Estimate(1);
modelParams.cb = modelParams.model.Coefficients.Estimate(2);


% sanity check
if modelParams.Rou > 1; error('Rouse > 1, error'); end % is this fair?


% evaluate the regression
nEvalPts = 50;
modelEvalZs = linspace(modelParams.flowDepth*0.05, modelParams.flowDepth, nEvalPts)';
modelEvalCs = modelParams.cb .* ( ((modelParams.flowDepth-modelEvalZs)./modelEvalZs) ./ modelParams.Hbb ) .^ modelParams.Rou;


% plot it
figure(); hold on;
plot(samplesTable.conc, samplesTable.sampleZ, 'ko')
plot(modelEvalCs, modelEvalZs, 'r-')