MATLAB Answers

CDF of log pearson

조회 수: 2(최근 30일)
Raj Arora
Raj Arora 2021년 7월 10일
편집: Paul 2021년 7월 11일
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

채택된 답변

Paul
Paul 2021년 7월 10일
편집: Paul 2021년 7월 11일
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
Raj Arora
Raj Arora 2021년 7월 11일
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.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by