How do I create a Weibull fit over an empirical cumulative distribution function (ecdf)?

clear all
close all
clc

%Data import
Windspeed_hourly_mean = xlsread('Windspeed_hourly.xlsx','B2:B8763');
WHM=Windspeed_hourly_mean;

%Parameters
threshold = 3;
windHeight = 50;
hubHeight = 150;

%Windspeed calculation at hub height
wind60 = WHM .* log(60/0.0001) ./ log(windHeight/0.0001);
windHub = wind60 .* (hubHeight ./ 60) .^ 0.115;

isBelow = (windHub < threshold); % gives us a vector of logical values
isBelow = [false; isBelow; false]; % make sure the vector starts and ends with false
transitions = diff(isBelow); % gives 1 where false->true and -1 where true->false 
starts = find(transitions==1); % indexes of where each period starts
ends = find(transitions==-1); % indexes of where each period ends
durations = ends - starts;
max(durations); % largest duration

t_delta = max(durations);

plot(windHub,'b')

n=histc(durations,1:t_delta);

[f,x]= ecdf(durations)

ecdf(durations)

Answer 1

[Mathieu NOE]

hello

this is the suggestion from the guy who has no toolboxes (almost !)

• alpha is a_sol in my code
• beta is b_sol
• M = 1 obviously

clear all
close all
clc
%Data import
Windspeed_hourly_mean = xlsread('Windspeed_hourly.xlsx','B2:B8763');
WHM=Windspeed_hourly_mean;
%Parameters
threshold = 3;
windHeight = 50;
hubHeight = 150;
%Windspeed calculation at hub height
wind60 = WHM .* log(60/0.0001) ./ log(windHeight/0.0001);
windHub = wind60 .* (hubHeight ./ 60) .^ 0.115;
isBelow = (windHub < threshold); % gives us a vector of logical values
isBelow = [false; isBelow; false]; % make sure the vector starts and ends with false
transitions = diff(isBelow); % gives 1 where false->true and -1 where true->false 
starts = find(transitions==1); % indexes of where each period starts
ends = find(transitions==-1); % indexes of where each period ends
durations = ends - starts;
max(durations); % largest duration
t_delta = max(durations);
n=histc(durations,1:t_delta);

[y,x]= homemade_ecdf(durations); %https://fr.mathworks.com/matlabcentral/fileexchange/32831-homemade-ecdf

% curve fit using fminsearch
f = @(a,b,x) (1- exp(-a*x.^b));
obj_fun = @(params) norm(f(params(1), params(2),x)-y);
sol = fminsearch(obj_fun, [0.5,1]); % "good" initial guess => will give good fit 
a_sol = sol(1);
b_sol = sol(2);
y_fit = f(a_sol, b_sol, x);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

plot(x, y_fit, '-',x,y, 'r .', 'MarkerSize', 40)
ylim([0 1.1]);
title(['Weibull Fit to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end


function [v_f,v_x] = homemade_ecdf(v_data)
    nb_data = numel(v_data);
    v_sorted_data = sort(v_data);
    v_unique_data = unique(v_data);
    nb_unique_data = numel(v_unique_data);
    v_data_ecdf = zeros(nb_unique_data,1);
    for index = 1:nb_unique_data
        current_data = v_unique_data(index);
        v_data_ecdf(index) = sum(v_sorted_data <= current_data)/nb_data;
    end
%     v_x = [v_unique_data(1); v_unique_data];
    v_x = [0; v_unique_data];
    v_f = [0; v_data_ecdf];
end