How to fit the cdf to a function?

16 vues (au cours des 30 derniers jours)
Hozifa
Hozifa le 22 Fév 2023
Modifié(e) : William Rose le 23 Fév 2023
Hi there,
I want my code to give the exponential fit as a function of the x-axis values, in my code, the fit function depends on the generated linespace,
Mnay thanks,
clear
clc
z=[
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
8.2000
8.4000
10.4000
10.4000
11.0000
11.2000
11.4000
11.4000
11.5000
11.5900
11.7000
11.8000
12.0000
12.2000
12.3400
12.4300
12.6000
12.9700
13.2100
13.2900
13.3400
13.3800
14.4000
15.4500
15.8200
15.9900
16.2600
16.4200
16.4400
16.8100
17.6100
17.7000
19.0400
20.8100
21.4000
22.0800
22.4000
22.4300
22.4700
22.5000
22.5000
22.5000
23.4800
23.5500
23.7500
24.1100
24.4800
24.6000
24.7000
26.0900
26.3000
26.4400
26.4500
26.4800
26.5100
27.0000
27.9000
28.3900
28.8500
29.5000
29.8500
30.3000
30.3000
31.0000
31.3000
31.3200
33.1000
33.5700
33.6500
33.9000
34.9300
36.2400
36.6000
37.0200
38.8200
39.6900
40.1600
40.2800
40.3000
40.6200
41.0300
41.1900
42.6800
42.7100
42.7500
43.6500
43.8000
43.9000
46.6000
46.6700
49.4100
49.5400
49.6800
49.8500
50.3000
53.9000
56.6000
60.3000
62.9000
72.9000
72.9000
72.9000
72.9000
72.9000
73.0000
73.0000
73.0000
73.0000
73.3000
73.3000
73.3000
73.3000
73.3000
73.3000
73.3000
73.3000
73.3000
73.3800
73.4000
73.4000
73.4000
73.4000
73.4000
73.4000
73.4000
73.4000
73.4000
73.4200
73.4500
73.6000
73.6000
73.6000
73.6000
73.6000
73.6000
73.6000
73.6000
];
cdfplot(z)
hold on
x = linspace(0,1,length(z));
% x=z;
fitfunc = fittype(@(B,C,x) 0.4*exp(B*x.^C));
x0 = [10 0.5];
f = fit(x',z,fitfunc,'StartPoint',x0)
coeffvals = coeffvalues(f)
plot(f(x),x,'-r',z,x,'b.')
% legend('fitted curve','data')
legend('LOS','Exponential fit')
set(findall(gcf,'type','line'),'LineWidth',1)
set(gca,'fontweight','bold','FontSize',12)
xlabel('Inter cluster delay (ns)','fontweight','bold','fontsize',12)
ylabel ('Probability','fontweight','bold','fontsize',12)

Réponse acceptée

William Rose
William Rose le 22 Fév 2023
Modifié(e) : William Rose le 23 Fév 2023
[edit: correct spelling and grammar]
I put the data in a text file so it is not in your code. It appears that the data are interpreted as being samples of a random variable, which have been sorted into ascending order.
Your question and the comments in the code indicate you want to fit an exponential distribution to the data:
The max. likelihood estimate (MLE) for lambda is 1/mean(x)
Therefore I have computed the CDF with this estimate, and I have added it to the plot. Of course your data shows signs of thresholding and saturation, which violate the assumptions used in the MLE derivation. Thresholding and saturation are types of data censoring. Google, or search matlab answers, for more info on advanced strategies for fitting censored data, if that is of interest.
In the code below, I computed the MLE and added the corresponding CDF to the plot. This method of fitting does not depend on a generated linespace vector, which was your original concern.
z=load('statdata.txt');
cdfplot(z)
hold on
x = linspace(0,1,length(z));
fitfunc = fittype(@(B,C,x) 0.4*exp(B*x.^C));
x0 = [10 0.5];
f = fit(x',z,fitfunc,'StartPoint',x0)
f =
General model: f(x) = 0.4*exp(B*x.^C) Coefficients (with 95% confidence bounds): B = 5.426 (5.396, 5.456) C = 0.582 (0.5531, 0.611)
coeffvals = coeffvalues(f)
coeffvals = 1×2
5.4257 0.5820
lambda=1/mean(z); %MLE for exponential distribution
cdf_mle=1-exp(-lambda*[1:100]); %CDF using MLE for lambda
plot(f(x),x,'-r',z,x,'b.',[1:100],cdf_mle,'-g')
% legend('fitted curve','data')
legend('','Exponential fit','LOS','MLE')
set(findall(gcf,'type','line'),'LineWidth',1)
set(gca,'fontweight','bold','FontSize',12)
xlabel('Inter cluster delay (ns)','fontweight','bold','fontsize',12)
ylabel ('Probability','fontweight','bold','fontsize',12)

Plus de réponses (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by