tsEva/tsGetReturnPeriodOfLevel.m at master · menta78/tsEva · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
function [ myRetPeriod, myRetPeriodCISup, myRetPeriodCIInf ] = tsGetReturnPeriodOfLevel( retPeriod, retLevel, retLevError, myLevel, varargin )
  % given a list of retPeriod and corresponding retLevel with error,
  % estimates the return period for a level myLevel, and the related confidence interval.

  args.logExtrap = true;
  args.linExtrap = false;
  args.cpP = .68; % by default estimating 1-sigma confidence interval
  args = tsEasyParseNamedArgs(varargin, args);
  cpP = args.cpP;
  linExtrap = args.linExtrap;
  logExtrap = args.logExtrap && ~linExtrap;

  [rlLow, rlHigh] = tsEstimateConfidenceIntervalOfRL(retLevel, retLevError, cpP);

  myLevNonNan = ~isnan(myLevel) & ~isinf(myLevel);

  if logExtrap
    myRetPeriod_ = tsInterp1Extrap(retLevel, retPeriod, myLevel(myLevNonNan), logExtrap)';
    try
      myRetPeriodCISup_ = tsInterp1Extrap(rlLow, retPeriod, myLevel(myLevNonNan), logExtrap)';
    catch
      myRetPeriodCISup_ = interp1(rlLow, retPeriod, myLevel(myLevNonNan));
      myRetPeriodCISup_(isnan(myRetPeriodCISup_)) = inf;
    end
    try
      myRetPeriodCIInf_ = tsInterp1Extrap(rlHigh, retPeriod, myLevel(myLevNonNan), logExtrap)';
    catch
      myRetPeriodCIInf_ = interp1(rlHigh, retPeriod, myLevel(myLevNonNan));
    end
  elseif linExtrap
    myRetPeriod_ = interp1(retLevel, retPeriod, myLevel(myLevNonNan), 'linear', 'extrap');
    myRetPeriodCISup_ = interp1(rlLow, retPeriod, myLevel(myLevNonNan), 'linear', 'extrap');
    myRetPeriodCIInf_ = interp1(rlHigh, retPeriod, myLevel(myLevNonNan), 'linear', 'extrap');
  else
    myRetPeriod_ = interp1(retLevel, retPeriod, myLevel(myLevNonNan));
    myRetPeriodCISup_ = interp1(rlLow, retPeriod, myLevel(myLevNonNan));
    myRetPeriodCIInf_ = interp1(rlHigh, retPeriod, myLevel(myLevNonNan));
  end

  myRetPeriod = ones(size(myLevel))*nan;
  myRetPeriod(myLevNonNan) = myRetPeriod_;

  myRetPeriodCISup = ones(size(myLevel))*nan;
  myRetPeriodCISup(myLevNonNan) = myRetPeriodCISup_;

  myRetPeriodCIInf = ones(size(myLevel))*nan;
  myRetPeriodCIInf(myLevNonNan) = myRetPeriodCIInf_;

end