Not familiar with the disribution, but can offer the following about this code:
% data
Q=[23.81
33.98
62.01
140.45
184.91
257.97
325.64
410.59
543.68];
% parameters
n = length(Q);
Z = log(Q);
mu = mean(Z);%%mean
Sigma = std(Z);%%Standard deviattion%
Skew = sum(n*((Z-mu).^3))/((Sigma.^3)*(n-1)*(n-2));%%Skewness%%
kurt = (((n*(n+1))/((n-1)*(n-2)*(n-3)))*(sum(((Z-mu)/Sigma).^4)))-((3*((n-1).^2))/((n-2)*(n-3)));
alpha_p = 4/(Skew^2);%%Shape_alpha_b%%
beta_b = (Sigma)/(alpha_p^0.5); %%scale_beta_1bya%%
y_a = mu-(Sigma*(alpha_p^0.5));%%location_y_m%%
% pdf
syms e positive
fun(e) = (1/(e*gamma(alpha_p).*abs(beta_b))).*(((log(e)-y_a)/beta_b).^(alpha_p - 1)).*(exp(-((log(e)-y_a)/beta_b)));
The support of the pdf is not given. Note that for e less than 3 fun(e) is complex, and f(3) ~ 0
double(fun(2))
double(fun(3))
So we'll assume the support of fun is 3 < e < inf.
The CDF is the integral of the pdf. Note that we must integrate from 3 to q to get P(3 < Q < q). Compare this to the original code
syms q positive
F(q) = simplify(int(fun(e),e,3,q));
Plot the CDF:
fplot(F(q),[3 1000])
That looks like a CDF (note that F(q) = 0 for q < 3). No idea if it's the CDF you're expecting.
Looks like it's very close to, if not exactly, the incomplete gamma function with support x >= exp(y_a)
hold on;
x = 3:10:1000;
plot(x,gammainc((log(x)-y_a)/beta_b,alpha_p),'ro')
