MATLAB Answers

CDF of log pearson

2 views (last 30 days)
Raj Arora
Raj Arora on 10 Jul 2021
Edited: Paul on 11 Jul 2021
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

Accepted Answer

Paul
Paul on 10 Jul 2021
Edited: Paul on 11 Jul 2021
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')
  1 Comment
Raj Arora
Raj Arora on 11 Jul 2021
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!