Reza - If you really want to fit a Pearson distribution to your moments,
I have an old Stata program to do it, but I don't remember exactly what
the output is. Basically if y is your data and x is a x-axis plotting
variable (can be the same as y or cover the same range as y butt say,
equally spaced) then typing
pearson y x
produces a Pearson density fy0. If you plot fy0 vs x and superimpose a
histomgram of y you should get a good approximation. The parameters that
produce fy0 are buried in a list of scalars produced by the program.
Maybe from a description of the Pearson fitting process you can figure
out what they are and what "Type" your distribution is. It's been so
long since I wrote this, I don't remember.
Al Feiveson
= ================================
program define pearson
version 7.0
syntax varlist [aweight] [if] [in]
args y x
qui summ `y' [`weight' `exp'] `if' `in' ,det
scalar mu1=r(mean)
scalar mu2=r(Var)
scalar skew=r(skewness)
scalar kurt=r(kurtosis)
scalar sig=r(sd)
local s "scalar"
`s' mu3=`s'(skew)*`s'(sig)^3
`s' mu4=`s'(kurt)*`s'(mu2)^2
`s' be1=`s'(skew)^2
`s' be2=`s'(kurt)
local be1 "scalar(be1)"
local be2 "scalar(be2)"
local mu2 "scalar(mu2)"
local sig "scalar(sig)"
local skew "scalar(skew)"
/* calculate Pearson diff-eq parameters for mean=0 */
tempvar xp
gen `xp'=`x'-scalar(mu1)
scalar Ap=10*`be2'-18-12*`be1'
scalar b0 = -(4*`be2'-3*`be1')*`mu2'/scalar(Ap)
scalar b1=-(`be2'+3)*`sig'*`skew'/scalar(Ap)
scalar b2=-(2*`be2'-3*`be1'-6)/scalar(Ap)
scalar a=scalar(b1)
local a "scalar(a)"
local b0 "scalar(b0)"
local b1 "scalar(b1)"
local b2 "scalar(b2)"
scalar B0=`b0'+`a'*`b1'+`a'*`a'*`b2'
scalar B1=`b1'+2*`a'*`b2'
scalar B2=`b2'
scalar D=B1*B1-4*B0*B2
di "D = ",`s'(D)
if D < 0 { /* complex roots - Type IV */
scalar del=sqrt(-D/(4*B2*B2))
scalar ga=B1/(2*B2)
tempvar X
gen `X'=`xp'-`s'(`a')+`s'(ga)
cap gen zx1=.
replace zx1=log(`X'*`X'+`s'(del)*`s'(del))/(2*B2)
cap gen zx2=.
replace zx2=`s'(ga)*atan(`X'/`s'(del))/(B2*`s'(del))
cap gen f0=.
replace f0=exp(zx1-zx2)
}
integ f0 `xp'
replace f0=f0/r(integral)
`s' top = `be1'*(`be2'+3)^2
`s' bot = 4*(2*`be2'-3*`be1'-6)*(4*`be2'-3*`be1')
`s' kap = `s'(top)/`s'(bot)
di "kappa = ",`s'(kap)
end
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Nick Cox
Sent: Friday, November 03, 2006 6:09 AM
To: [email protected]
Subject: st: RE: Creating a distribution from moments
This is called the method of moments
and Karl Pearson thought, more or less, that it was the best way to fit
a distribution, so long as you know what kind of distribution you are
fitting. But it has never really recovered from criticisms made by
Ronald A. Fisher and others several decades ago, except that we still
use it routinely in problems like plugging the mean and the standard
deviation into a Gaussian (so-called normal) formula.
(That is usually very close to the maximum likelihood solution, but most
summary programs, like -summarize-, use n - 1 not n as divisor for the
variance estimator.)
Yes, you can have a go at this, but you need to look at books like the
multivolume reference by N.L. Johnson, S. Kotz and friends to see the
ways to do it. Usually there is some indirectness:
for example, I can't recall any named distribution for which one of the
parameters _is_ the skewness or the kurtosis: rather if you have k
parameters, you typically end up with k simultaneous equations in those
parameters and the first k moments will be useful in solving those.
(Strictly, there is a terminology problem here:
the skewness and kurtosis, at least as
named in Stata, are ratios derived from moments, not moments themselves,
but we know what you mean.)
More important, this is only the method of choice if these summaries are
_all_ you have to go on. If you have the data, use all the data.
Nick
[email protected]
Reza C Daniels
> I have summary statistics for four moments of the density of y:
> mean=877; std. deviation=611; skewness=0.658; kurtosis=2.278. Is it
> possible for me to use this information to generate a hypothetical
> density whose four moments approximate these values?
>
> The analogy here would be creating, for example, a Gaussian
> distribution
> by:
>
> set obs 1000
> gen z1=invnorm(uniform())
>
> My question relates to when one does not want a standard density but
> possibly some parameterization of a Beta or Gamma that fits the four
> moments.
>
> I have searched through the probability distributions and density
> functions in Stata (I'm using version 8.2SE), but it is not
> immediately obvious to me that I can do this.
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
