Code covered by the BSD License  

Highlights from
pearspdf

5.0

5.0 | 4 ratings Rate this file 26 Downloads (last 30 days) File Size: 4.83 KB File ID: #26516

pearspdf

by

 

28 Jan 2010 (Updated )

Pearson system probability distributions

| Watch this File

File Information
Description

% pearspdf
% [p,type,coefs] = pearspdf(X,mu,sigma,skew,kurt)
%
% Returns the probability distribution denisty of the pearsons distribution
% with mean `mu`, standard deviation `sigma`, skewness `skew` and
% kurtosis `kurt`, evaluated at the values in X.
%
% Some combinations of moments are not valid for any random variable, and in
% particular, the kurtosis must be greater than the square of the skewness
% plus 1. The kurtosis of the normal distribution is defined to be 3.
%
% The seven distribution types in the Pearson system correspond to the
% following distributions:
%
% Type 0: Normal distribution
% Type 1: Four-parameter beta
% Type 2: Symmetric four-parameter beta
% Type 3: Three-parameter gamma
% Type 4: Not related to any standard distribution. Density proportional
% to (1+((x-a)/b)^2)^(-c) * exp(-d*arctan((x-a)/b)).
% Type 5: Inverse gamma location-scale
% Type 6: F location-scale
% Type 7: Student's t location-scale
%
% Examples
%
% See also
% pearspdf pearsrnd mean std skewness kurtosis
%

It's a modification of the pearsrnd function

Required Products Statistics Toolbox
MATLAB release MATLAB 7.9 (R2009b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (25)
22 May 2014 David Holtschlag  
23 May 2012 khoshi

Dear all
I need to apply Pearson type 5 distribution
to number data following code includes all types how can i extract just type 5?

14 Aug 2011 Pierce Brady

Hey guys,
I appreciate all the comments etc.
But i don't have the time to fix these corrections. Its open source so feel free to change it and reuse it yourself.

Carsten Allefeld has made some modifications that fix the normalisation problem. Except in the case of the type 5 distribution.

The sizeOut error has also been fixed.

10 Nov 2010 Daryl

Question:

Occasionally I get an error Undefined function or variable 'sizeOut'. I can't find this value or variable in the code. Can anyone help?

13 Oct 2010 Oleg Nazarevych

See good code in R
http://cran.r-project.org/web/packages/PearsonDS/

02 Oct 2010 Oleg Nazarevych

edit your file:

122: case 1
....
136: p = betapdf(X,m1+1,m2+1)/sigma;
137: p=p/55;

177: case 6
...
187: p = fpdf(X,nu1,nu2)/sigma;
188: p=p/3;

and

195: p = fpdf(X,nu1,nu2)/sigma;
196: p=p/3;

it is not correct, but plot norm PDF hold on Hist (ecdfhist)

30 Sep 2010 Oleg Nazarevych

and not correct for Pearson Type I
coefs = [1.2197 0.30669 -0.07322]

30 Sep 2010 Oleg Nazarevych

for m = {100,10,0.5,3};
not correct

29 Sep 2010 Pierce Brady

add this at line 164

p = p/sigma;

just after
p = pearson4pdf(X,m,nu,a,lambda);

it seems the the pdfs are off by this factor

28 Sep 2010 Oleg Nazarevych

Ok, try your code for m = {0,10,0.5,5};

28 Sep 2010 Pierce Brady

% You need to normalise the histogram properly
[Y,X] = hist(pearsrnd(m{:},100000,1),100)
dX = diff(X)
% Forgot to include this line
dX = dX(1)
bar(X,Y/(sum(Y)*dX))
hold on;
plot(-6:0.01:6,pearspdf(-6:0.01:6,m{:}),'r')

28 Sep 2010 Oleg Nazarevych

Try his code:
%m = {0,1,0,3}; %GOOD
%m = {100,1,1,5}; %GOOD
m = {100,10,1,5}; %BAD. WHY ?! same else scale....
[x,type, coefs] = pearsrnd(m{:},10000,1);
m2 = [mean(x) std(x) skewness(x) kurtosis(x)];
clf;[Y,X] = hist(pearsrnd(m{:},10000,1),100);
%dX = diff(X);
%bar(X,(Y/(sum(Y)*dX)));
[Fxi,Xxi] = ecdf(x);
ecdfhist(Fxi,Xxi,100) %Norm hist
hold on;
xx=min(x):(max(x)-min(x))/1000:max(x);
plot(xx,pearspdf(xx,m{:}),'r');
title(['\bfType (', num2str(type), '), coefs: ',num2str(coefs)])

28 Sep 2010 Oleg Nazarevych

Thax for replay.
Same wrong, arrya size dX[99]<>Y[100]

May be this is true:
%dX = diff(X);
%bar(X,(Y/(sum(Y)*dX)));
[Fxi,Xxi] = ecdf(x);
ecdfhist(Fxi,Xxi,100) %Norm hist
hold on;
plot(-6:0.01:6,pearspdf(-6:0.01:6,m{:}),'r');

28 Sep 2010 Pierce Brady

m and m2 will be different. pearsrnd generates random numbers that follow the distribution.
if you generate more numbers the values will begin to match much better.
i.e.
m = {0,1,0,3};%Norm Distrib
[x,type, coefs] = pearsrnd(m{:},100000,1);
m2 = [mean(x) std(x) skewness(x) kurtosis(x)]

% better

m = {0,1,0,3};%Norm Distrib
[x,type, coefs] = pearsrnd(m{:},100000000000000,1);
m2 = [mean(x) std(x) skewness(x) kurtosis(x)]
% Won't have a enough memory
% m and m2 should be very close if not the same

% As you you generate more and more numbers, the error
% between m and m2 will become smaller and smaller

% You need to normalise the histogram properly
[Y,X] = hist(pearsrnd(m{:},100000,1),100)
dX = diff(X)
bar(X,Y/(sum(Y)*dX))
hold on;
plot(-6:0.01:6,pearspdf(-6:0.01:6,m{:}),'r')

27 Sep 2010 Oleg Nazarevych

1: m = {mean(x),std(x),skewness(x),kurtosis(x)};
2:[r,type,coefs] = pearsrnd(m{:},10000,1);
3:hold on; [Fxi,Xxi] = ecdf(x);
4:ecdfhist(Fxi,Xxi,k);
5:Yxi=pearspdf(sort(x),m{:}); plot(sort(x),Yxi,'-y','LineWidth',1)

Hist genering by pearsrnd and plot of p(x) in same else scale, why ?
Where my bug ?

Q1: How plot p(x) (by pearspdf()) hold on Norm Hist(pearsrnd()) ?

27 Sep 2010 Oleg Nazarevych

m = {0,1,0,3};%Norm Distrib
[x,type, coefs] = pearsrnd(m{:},1000,1);
m2 = [mean(x) std(x) skewness(x) kurtosis(x)]
m2 =
-0.0404 1.0045 0.0485 3.1589

Why m <> m2 ?!

23 Sep 2010 Oleg Nazarevych

*** sigma=std(x)^2;

23 Sep 2010 Oleg Nazarevych

From Matlab help:
[r,type,coefs] = pearsrnd(mu,sigma,skew,kurt,m,n) returns the type of the specified distribution within the Pearson system.
(!) Set m and n to 0 to identify the distribution type without generating any random values.

Example
x=x1,x2,....xn;
mu=mean(x);
sigma=std(x);
skew=skewness(x);
kurt=kurtosis(x);
[r,type,coefs] = pearsrnd(mu,sigma,skew,kurt,0,0)

21 Sep 2010 Oleg Nazarevych

ok. thanx

In string 57: p = NaN(sizeOut,outClass);

if (sigma < 0) || (beta2 <= beta1 + 1)
p = NaN(sizeOut,outClass);
type = NaN;
coefs = NaN(1,3,outClass);
return
end

Not define - sizeOut
May be have - sizeOut=size(X) ?

But type of Pearson Curve define from parameter Kapa:
1) Kapa=skew^2*(kurt+3)^2/(4*(4*kurt-3*skew^2)*(2*kurt-3*skew^2-6));
or
2) Kapa2= - (skew^2*(K+2)^2)/(16*(K+1));

21 Sep 2010 Pierce Brady

you'll have to write them yourself so

21 Sep 2010 Oleg Nazarevych

toolbox\stats\
can find pdf, cdf and inv-files for all distribution.

We have't its for pears.... ^( Onle pearspdf.m (your local file)

I want have pearcdf.m, pearsinv.m for using mle.m, adding Pearson for Maxumum Likehood estimation (MLE metod).

21 Sep 2010 Oleg Nazarevych

For example
mle(x,'distribution',Pearson)

21 Sep 2010 Pierce Brady

Hey,
I'm not sure where you can find pearsinv
this file is just a modifcation to the pearsrnd function.

but for the cdf
1. caculate the distribution using pearspdf
2. use cumsum on the result to obtain cdf
3. normalise by the total sum, to ensure limit is 0 to 1

a = cumsum(pearspdf(-3:.01:3,0,0,0,3))
b = sum(pearspdf(-3:.01:3,0,0,0,3))
cdf = a/b;

20 Sep 2010 Oleg Nazarevych

Help.
Where can i finf pearscdf.m, pearsinv.m files for Pearson Distribution ?
Mail me - taltek.te@gmail.com

Thanx

11 Jun 2010 John Peterson

Thanks, I've been looking for this. Would it be possible to always have the probability sum to one? Or to create a pearscdf function?

Updates
29 Sep 2010

Added the line at 164
p = p/sigma;

as the type 4 distributions where off by this factor

14 Aug 2011

Many thanks to Carsten Allefeld who has modified the code so that pdfs are normalised correctly.

The type 5 distribution is not normalised correctly, so a warning is thrown to indicate this.

Contact us