Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

# Re: st: distribution test

 From Maarten Buis <[email protected]> To [email protected] Subject Re: st: distribution test Date Tue, 30 Aug 2011 16:45:50 +0200

```On Tue, Aug 30, 2011 at 11:13 AM, Nick Cox  wrote:
> In addition to Maarten's good advice, and setting aside the question
> of covariates, not mentioned in the original, at least as quoted here:
>
> There is a direct method to check for fit to an exponential
> distribution: a quantile-quantile plot.

Another alternative is the probability-probability plot. The QQ-plot
is good at detecting deviations in the tails of the distribution,
while the PP-plot is good at detecting deviations at the center of the
distribution. Just like the QQ-plot, I like to add a number of
simulations assuming the model is correct to help assess how much
deviation from the 45 degree line is to be expected. Below is an
example how I would create such a graph in Stata:

*------------------------- begin example ----------------------
// create some expnential data
local lambda = 2
drop _all
set obs 500
gen y = -1/`lambda'*ln(1-runiform())
label var y "observed"

// estimate parameter
sum y, meanonly
local lambdahat = 1/r(mean)
di as txt "ML estimate of lambda is: " ///
as result `lambdahat'

// create variables for pp-plot
gen byte touse = !missing(y)
sort y
gen F = 1-exp(-`lambdahat'*y)
label var F "observed"
gen Psubi = sum(touse)
replace Psubi = cond(F>=.,.,Psubi/(Psubi[_N]+1))

// create 20 random variables assuming the model is correct
forvalues i = 1/20 {
gen y`i' = -1/`lambdahat'*ln(1-runiform()) if touse
gen F`i' = 1-exp(-`lambdahat'*y`i')
sort y`i'
drop y`i'
gen Psubi`i' = sum(touse)
replace Psubi`i' = cond(F>=.,.,Psubi`i'/(Psubi`i'[_N]+1))
#delim ;
local gr `"`gr' || line F`i' Psubi`i',
sort lstyle(solid) lcolor(gs12)"' ;
#delim cr
}

// create the graph
twoway `gr'                               || ///
scatter F Psubi,                       ///
aspect(1) msymbol(oh)               || ///
function reference = x,                ///
lstyle(solid)                          ///
legend(order( 1 "simulations" 21 22 )) ///
xtitle("Empirical Pr(Y {&le} y) = i/(N+1)") ///
ytitle("theoretical Pr(Y {&le} y)")
*-------------------------- end example -----------------------
(For more on examples I sent to the Statalist see:
http://www.maartenbuis.nl/example_faq )

Hope this helps,
Maarten

--------------------------
Maarten L. Buis
Institut fuer Soziologie
Universitaet Tuebingen
Wilhelmstrasse 36
72074 Tuebingen
Germany

http://www.maartenbuis.nl
--------------------------
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```