Re: st: distribution test

[Nick Cox <njcoxstata@gmail.com>]

Also -pexp- (SSC). (It does not do the simulations.)

Nick

On Tue, Aug 30, 2011 at 3:45 PM, Maarten Buis <maartenlbuis@gmail.com> wrote:
> 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 )
>

*
*   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/