model parameter estimation from RMSE between modeled outputs and observations

Hello,

I am a hydrologist, and trying to estimate soil parameters by using RMSE between modeled outputs (streamflow) and observed ones.

For example, a watershed model is applied from 1980 to 2000 with an initial set of parameters. RMSE can be calcuated between modeled streamflow and observed ones.

With the RMSE values, is there any matlab function to search a new set of parameters to try to reduce RMSE?

The curve for RMSE over an iteration would be ideally converged in a global minimization point.

Thanks!

Answer 1

[Jeff Miller]

I assume the model is too complex for regression, etc. In that case, you might be able to do this with fminsearch, if there are not too many parameters. You write a function to compute RMSE for any given parameter vector, and fminsearch will try to find the parameter values that will minimize that function.

[DK]

Thanks.

Step1) Simulated with the paramters and observed outputs give multiple RMSE values.

Step2) With given RMSE values, I am trying to search the best set of the parameters to minimize the RMSE values

Is there any example for doing this? I am not trying to explain it physically at this stage.

[Jeff Miller]

This forum seems to have a number of examples using fminsearch in various ways. Not sure if any of them minimize RMSE specifically, but the general principles are the same regardless of what error function you compute.

[DK]

Dear Jeff,

Below is my code to conduct a sensitivity test with varying two parameters (dsnowi, grfactor) of the function (msnow). RMSE (YEH, and Tb) is calculated at the end of the loop, but it would be great, preceeding guesses of the parameters help to search the optimized ones with the minimized RMSE.

function rresult = calibgrain(DBn_2,dateob,Tb)

dss = linspace(0.01,3,100);

ggf = linspace(0.5,5,100);

rresult = [];

for ii = 1:length(dss)

for jj = 1:length(ggf)

dsnowi = dss(ii);grfactor = ggf(jj);

[dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);

rmsecomp = calcrmse(YEH,dateob,Tb);

rresult = [rresult ;dsnowi grfactor rmsecomp];

jj

end

ii

%jj

end

Answer 2

[Jeff Miller]

Something like this. If it is too slow, note that you can use 'optimset' to pass fminsearch options that make it finish faster (at the loss of a little precision in your parameter estimates).

function rresult = calibgrain(DBn_2,dateob,Tb)
 
dss = linspace(0.01,3,100);
ggf = linspace(0.5,5,100);
 
rresult = [];
for ii = 1:length(dss)
    for jj = 1:length(ggf)
        dsnowi = dss(ii);grfactor = ggf(jj);
        startparms = [dsnowi,grfactor];
        bestparms = fminsearch(@RMSE1,startparms);
        dsnowi = bestparms(1);
        grfactor = bestparms(2);
        [dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);
        rmsecomp = calcrmse(YEH,dateob,Tb);
        rresult = [rresult ;dsnowi grfactor rmsecomp];
        jj
    end
    ii
    %jj

    function thisrmse = RMSE1(parms)
    % nest this function inside calibgrain so that it has
    % access to all of calibgrain's variables.
    dsnowi = parms(1);
    grfactor = parms(2);
    [dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);
    thisrmse = calcrmse(YEH,dateob,Tb);
    end

end

[DK]

Error: File: calibgrain.m Line: 21 Column: 5

Function definition is misplaced or improperly nested.

Maybe, I was not educated with 'function eval', but this stil has the error associated with the sub-function with RMSE1.

[DK]

Dear Jeff,

I modified the code, but RMSE1 should have fixed inputs, DBn_2, dateob, and Tb.

How can I put them into the funsciton RMSE1, which is also evaluated in the main function, calibgrain?

function rresult = calibgrain(DBn_2,dateob,Tb)

dss = linspace(0.01,3,100);

ggf = linspace(0.5,5,100);

rresult = [];

for ii = 1:length(dss)

for jj = 1:length(ggf)

dsnowi = dss(ii);grfactor = ggf(jj);

startparms = [dsnowi,grfactor];

bestparms = fminsearch(@RMSE1,startparms);

dsnowi = bestparms(1);

grfactor = bestparms(2);

[dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);

rmsecomp = calcrmse(YEH,dateob,Tb);

rresult = [rresult ;dsnowi grfactor rmsecomp];

%jj

end

%ii

%jj

end

end

function thisrmse = RMSE1(parms)

% nest this function inside calibgrain so that it has

% access to all of calibgrain's variables.

dsnowi = parms(1);

grfactor = parms(2);

[dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);

thisrmse = calcrmse(YEH,dateob,Tb);

end

[Jeff Miller]

Sorry, I messed up the alignment of your end's. It should look like this:

function rresult = calibgrain(DBn_2,dateob,Tb)
    
    dss = linspace(0.01,3,100);
    ggf = linspace(0.5,5,100);
    
    rresult = [];
    for ii = 1:length(dss)
        for jj = 1:length(ggf)
            dsnowi = dss(ii);grfactor = ggf(jj);
            startparms = [dsnowi,grfactor];
            bestparms = fminsearch(@RMSE1,startparms);
            dsnowi = bestparms(1);
            grfactor = bestparms(2);
            [dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);
            rmsecomp = calcrmse(YEH,dateob,Tb);
            rresult = [rresult ;dsnowi grfactor rmsecomp];
        end    % end of jj loop
    end % end of ii loop
    
    function thisrmse = RMSE1(parms)
        % nest this function inside calibgrain so that it has
        % access to all of calibgrain's variables.
        dsnowi = parms(1);
        grfactor = parms(2);
        [dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH, YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);
        thisrmse = calcrmse(YEH,dateob,Tb);
    end  % end of function RMSE1
    
end % end of function calibgrain

RMSE1 has access to all of calibgrain's parameters because it is 'nested inside' the overall calibgrain function. But RMSE1 can't be inside one of the 'for' loops as I had it originally.

[DK]

Hi Jeff,

You are very helpful! By the way, I got the error message at 200th line of fminsearch

Not enough input arguments.

Error in calibgrain/RMSE1 (line 25)

[dsns_2,lws,btotals,phis,dsnows_y2_r1,tsnows,sds_0,kss,tsss,gds,outmelts,YEH,YEV,YEVH,

YEHV,result,compaction,ritisi,stor,ro] = msnow(1,6552,DBn_2,dsnowi,grfactor);

Error in fminsearch (line 200)

fv(:,1) = funfcn(x,varargin{:});

Error in calibgrain (line 11)

bestparms = fminsearch(@RMSE1,startparms);

Any solution for this?

Thank you very much!

[Jeff Miller]

Sorry, I don't see what is causing that. Is it possible that your original call to calibgrain (outside of all this code) did not pass its required 3 parameters?