Re: [Gretl-devel] Goodness of Fit test to normal distribution and QQ plot
Allin Cottrell schrieb:
On Mon, 30 Nov 2009, Sven Schreiber wrote:
> Costia schrieb:
>> 2. I haven't found a Q-Q plot test in GRETL, (like in SAS e.g.). I think
>> it would be a good thing to add it.
> I agree. It's not difficult to use a quick-and-dirty script to get it
> done, but properly handling missing values etc. needs more care.
I haven't messed much with Q-Q plots. Can anyone confirm if the
following script is correct? If so, then I guess we could do it
as a built-in at some point.
sorry for not answering earlier...
<script>
nulldata 500
scalar n = $nobs
# artificial series to test
genr y = sort(normal())
series Qemp, Qnorm
scalar p, N, nl, nh
loop i=1..n -q
p = i/(n+1)
N = (n + 1) * p - 1;
given your def. of p, wouldn't that be equiv. to N = i-1 ?
nl = floor(N);
accordingly redundant? (nl = i-1)
nh = ceil(N);
ditto?
if (nh == 0 || nh == n)
Qemp[i] = NA
don't understand the reason for this: here it seems the first item in
the empirical quantiles will always be set to NA?
else
Qemp[i] = quantile(y, p)
endif
Qnorm[i] = critical(N, 1-p)
or rather Qnorm[i] = 1-critical(N,p)? (could make a difference for
discrete distributions, not for continuous as here)
endloop
and the 45-degree line is needed for comparison
altogether I think something like this also would do the trick (bugs
probably included, even if it's only two lines...):
matrix mp = transp(seq(1,n))/n
# (col vector of cumulative relative frequencies for n obs)
matrix mQtheory = invcdf(N,mp)
# BTW, the help about invcdf() says P(X<x), but shouldn't it be P(X<=x)?
and your y (the sorted input series) can already be plotted against
mQtheory to give the QQ plot (if I'm not missing something right now
which may well be the case since I'm a bit tired). You might argue that
this doesn't really deal with ties, but from the various possible
definitions I've seen of quantiles per se and how to deal with ties I
think it's actually legitimate.
Why not just copy the convention and maybe the code from R?
good night,
sven