CDF of log pearson

Can somebody explain why when I use this code, I get CDF as negative and decreasing function (magnitude)? It should be non decreasing function

where Q include the data given below

23.81

33.98

62.01

140.45

184.91

257.97

325.64

410.59

543.68

I have integrated the pdf to calculate cdf of log pearson type III distribution

MATLAB CODE

Q = load('F:\DISCHARGE data.txt')

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%%

syms e

fun = (1/(e*gamma(alpha_p).*abs(beta_b))).*(((log(e)-y_a)/beta_b).^(alpha_p - 1)).*(exp(-((log(e)-y_a)/beta_b)));

wer = vpa(simplify(int(fun,e)))

CDF = round(double(subs(wer,Q)),4)

OUTPUT

-0.9674

-0.9175

-0.7612

-0.4678

-0.3738

-0.2743

-0.2158

-0.1669

-0.1196

Answer 1

[Paul]

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))

ans = -6.5190e-09 - 5.3012e-08i

double(fun(3))

ans = 7.2787e-16

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')

[Raj Arora]

Thanks paul for your suggestion and yes you are absolutely correct it is incomplete gamma function and the distribution above is log pearson type III distribution.