water_level.m

clear all; clc; 
figure(1), clf; 
 
data    =   load('waterstand_HOH.txt');   % loading water levels 
afmeting    =   size(data); 
som    =   afmeting(:,1);  % size of the dataset 
replace    =   true;  % Bootstrap with replacing 
jaar    =   98;  % Number of years over which    
         the dataset is created 
Y    =   sort(data);  
N    =   histc(data,Y); 
Y(N==0)    =   []; 
N(N==0)    =   []; 
sizeN    =   size(N); 
k    =   sizeN(:,1) 
 
for i = 1:k 
    data_orig(i)  =   sum(N(i:k)); 
end 
 
CumHOH_orig  =   data_orig';  
Levels_orig   =   Y; 
CumHOH_aantal_orig =  CumHOH_orig/jaar;  
 
for n = 1:k 
    CumHOH_jaar_orig(n) = 1/CumHOH_aantal_orig(n);  
end 
 
CumHOH_jaar_rot_orig =   CumHOH_jaar_orig';  % Cumulative probabilities of  
        HOH water levels 
for z = 1:1000    % Bootstrap 
    data  =   0; 
    CumHOH_jaar   =   0; 
    CumHOH_jaar_rot   =   0; 
    HOH  =   randsample(data1,som,replace); % Creating dataset out of  
      the data from the HOH  
      dataset  
    Y   =   sort(HOH);  
    N   =   histc(HOH,Y); 
    Y(N==0)   =   []; 
    N(N==0)   =   []; 
    sizeN  =   size(N); 
    k   =   sizeN(:,1); 
 
    for i = 1:k 
        data(i)   =   sum(N(i:k)); 
    end 
    CumHOH   =   data'; 
    Levels   =   Y; 
    CumHOH_aantal   =   CumHOH/jaar; 
 
    for n = 1:k 
        CumHOH_jaar(n)   = 1/CumHOH_aantal(n); 
    end 
 
    CumHOH_jaar_rot   =   CumHOH_jaar'; 
    p   =   polyfit(log(CumHOH_jaar_rot), Levels, 1); 
    A(z)  =   p(:,1); 
    B(z)  =   p(:,2); 
 
    for x = 1:1:10000 
        waterhoogte(x)   =   A(z)*log(x)+B(z); 
    end 
 
    water10000(z)   =   waterhoogte(x); 
    x   =   1:1:10000; 
    line(x,waterhoogte(x)); 
    hold on; 
    grid on; 
         
end 

water10000  =   sort(water10000); 
    water10000  =   water10000'; 
    parmhat   =   lognfit(water10000); 
    mu  =   parmhat(:,1); 
    sigma   =   parmhat(:,2); 
    v   =   300:1:600; 
    p   =   0.001:0.001:1; 
    CDF  =   logncdf(v,mu,sigma); 
    PDF  =   lognpdf(v,mu,sigma); 
    LINV  =   logninv(p,mu,sigma); 
    water_bovengrens   =   LINV(975); 
    water_ondergrens   =   LINV(25); 
 
hold on; 
scatter(CumHOH_jaar_rot_orig,Levels_orig,'.r'); 
plot(10000,water_bovengrens,'or'); 
plot(10000,water_ondergrens,'or');