Now I need to doubly apologize - for twice hastily posting something
that doesn't work. Apparently -pearr.ado- only works for Type IV also.
It probably can be fixed without too much trouble so I'll try this and
re-post it when (and if) I can get it to work properly.
Al F.
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Feiveson,
Alan H. (JSC-SK311)
Sent: Friday, November 03, 2006 8:21 AM
To: [email protected]
Subject: st: RE: pearson (correction)
I'm sorry - I think I included the wrong program in my last message -
that one only works for Type IV.
Here is another version that I think works in general. I think the third
argument `x0' is supposed to be an offset. You can probably get around
this by something like
gen x0=0
pearr y x x0
Al F.
==============================================
program define pearr
version 7.0
syntax varlist [aweight] [if] [in]
args y x x0
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)
/*
scalar mu1=-.190783898
scalar mu2=3.238424951
scalar be1=0.829135838
scalar be2=4.862944362
*/
scalar sig=sqrt(scalar(mu2))
local be1 "scalar(be1)"
local be2 "scalar(be2)"
local mu2 "scalar(mu2)"
local mu3 "scalar(mu3)"
local mu4 "scalar(mu4)"
local sig "scalar(sig)"
local skew "scalar(skew)"
`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)
scalar h1=`s'(mu1)-`x0'
local h1 "scalar(h1)"
scalar S=4*`h1'*`mu2'/`mu3'
local S "scalar(S)"
scalar r = 2+`S'
local r "scalar(r)"
scalar asq=(`r'-1)*`mu2'-`h1'*`h1'
scalar a=sqrt(scalar(asq))
local a "scalar(a)"
scalar v=-`r'*`h1'/`a'
scalar m=(`r'+2)/2
local v "scalar(v)"
local a "scalar(a)"
local m "scalar(m)"
scalar L=scalar(r)*scalar(r)+scalar(v)*scalar(v)
local L "scalar(L)"
scalar top4=3*(scalar(a)^4)*`L'*((`r'+6)*`L'-8*`v'*`v')
scalar bot4=(`r'^4)*(`r'-1)*(`r'-2)*(`r'-3)
scalar rhs4=scalar(top4)/scalar(bot4)
di "mu4, rhs ",scalar(mu4)," ",scalar(rhs4)
scalar top3=-4*`L'*`v'*scalar(a)^3
scalar bot3=(`r'^3)*(`r'-1)*(`r'-2)
scalar rhs3=scalar(top3)/scalar(bot3)
di "mu3, rhs ",scalar(mu3)," ",scalar(rhs3)
scalar rhs2=`L'*scalar(a)*scalar(a)/(scalar(r)*scalar(r)*(scalar(r)-1))
di "mu2, rhs ",scalar(mu2)," ",scalar(rhs2)
scalar rhs1=-scalar(a)*`v'/`r'
scalar lhs1=scalar(mu1)-`x0'
di "mu1-x0, rhs ",scalar(lhs1)," ",scalar(rhs1)
if `s'(kap) > 0 { /* Type IV */
tempvar X
gen `X'=`x'-`x0'
qui cap gen zx=. `if' `in'
qui replace zx=`m'*log(1+(`X'/`a')^2)+`v'*atan(`X'/`a') `if' `in'
qui cap gen f0=. `if' `in'
qui replace f0=exp(-zx) `if' `in'
}
qui integ f0 `x' `if' `in'
qui cap gen fh_mom=. `if' `in'
qui replace fh_mom=f0/r(integral) `if' `in'
end
/*
scalar r = 6*(`be2'-`be1'-1)/(2*`be2'-3*`be1'-6)
scalar A0=12*`be2'+12*`be1'+36
scalar A1=10*`be2'+30+12*`be1'
scalar A2=2*`be2'-6+3*`be1'
scalar D=A1*A1-4*A0*A2
scalar r1=(A1-sqrt(D))/(2*A2)
scalar r2=(A1+sqrt(D))/(2*A2)
di "r1, r2 ",scalar(r1)," ",scalar(r2)
scalar r = r1
local r "scalar(r)"
scalar bot=16*(`r'-1)-`be1'*(`r'-2)*(`r'-2)
local bot "scalar(bot)"
scalar v=`r'*(`r'-2)*sqrt(`be1')/sqrt(`bot')
scalar a = sqrt(`mu2'*`bot'/16)
*/
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Feiveson,
Alan H. (JSC-SK311)
Sent: Friday, November 03, 2006 8:08 AM
To: [email protected]
Subject: st: RE: RE: Creating a distribution from moments
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/
*
* 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/
