Windspeed_hourly_mean = xlsread('Windspeed_hourly.xlsx','B2:B8763');
WHM=Windspeed_hourly_mean;
wind60 = WHM .* log(60/0.0001) ./ log(windHeight/0.0001);
windHub = wind60 .* (hubHeight ./ 60) .^ 0.115;
isBelow = (windHub < threshold);
isBelow = [false; isBelow; false];
transitions = diff(isBelow);
starts = find(transitions==1);
ends = find(transitions==-1);
durations = ends - starts;
t_delta = max(durations);
n=histc(durations,1:t_delta);
[y,x]= homemade_ecdf(durations);
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]);
y_fit = f(a_sol, b_sol, x);
Rsquared = my_Rsquared_coeff(y,y_fit);
plot(x, y_fit, '-',x,y, 'r .', 'MarkerSize', 40)
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)
sum_of_squares = sum((data-mean(data)).^2);
sum_of_squares_of_residuals = sum((data-data_fit).^2);
Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
function [v_f,v_x] = homemade_ecdf(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;
v_x = [0; v_unique_data];