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
