[See at end] On 05-Apr-2012 00:34:30 Peter Ehlers wrote: > Hi Ted, > > On 2012-04-04 15:06, Ted Harding wrote: >> Greetings! >> I want to have the coefficients that R uses in shapiro.test() >> for the Shapiro-Wilk test for a prticular sample size, i.e. >> the a[i] in >> >> W = Sum(a[i]*x[i])/(Sum(x[i] - mean(x))^2) >> >> (where the x[i] are sorted). Two questions: >> >> Q1: >> Is there a readymade R function from which I can extract these? >> >> Q2: >> I was wondering if I might be able to modify the code for the >> function shapiro.test() so as to obtain these. When I enter >> >> shapiro.test >> >> I get: >> >> function (x) >> { >> DNAME<- deparse(substitute(x)) >> x<- sort(x[complete.cases(x)]) >> stopifnot(is.numeric(x)) >> n<- length(x) >> if (n< 3 || n> 5000) >> stop("sample size must be between 3 and 5000") >> rng<- x[n] - x[1L] >> if (rng == 0) >> stop("all 'x' values are identical") >> if (rng< 1e-10) >> x<- x/rng >> n2<- n%/%2L >> sw<- .C(R_swilk, init = FALSE, as.single(x), n, n1 = n, >> n2, a = single(n2), w = double(1), pw = double(1), ifault = >> integer(1L)) >> if (sw$ifault&& sw$ifault != 7) >> stop(gettextf("ifault=%d. This should not happen", sw$ifault), >> domain = NA) >> RVAL<- list(statistic = c(W = sw$w), p.value = sw$pw, method = >> "Shapiro-Wilk normality test", >> data.name = DNAME) >> class(RVAL)<- "htest" >> return(RVAL) >> } >> <environment: namespace:stats> >> >> >> So, on the off-chance that the variable 'sw' computed therein might >> contain something useful, I changed "return(RVAL)" to "return(sw)", >> just in case the coefficients might be lurking as a component of sw, >> and used this to define a function SW_ted(). I then ran >> >> SW_ted(rnorm(30)) >> # Error in SW_ted(rnorm(30)) : object 'R_swilk' not found >> >> Since shapiro.test(rnorm(30)) works perfectly, and since the >> "stats:" namespace is already present, I am wondering why >> "object 'R_swilk' not found" when it clearly can be found by >> shapiro.test(). >> >> So why is it that in the ".C" call: >> >> sw<- .C(R_swilk, ... ) >> >> my modifiction of shapiro.test() doesn't find it? >> >> (No doubt there is some dumb oversight behind this, but I'd be >> grateful to be told what it is)! > > No dumb oversight that I can see, but: > Q1: I have no idea, but I doubt that there is an existing function. > > Q2: I think that the environment of the newly created function > SW_ted needs to be set appropriately; try this: > > environment(SW_ted) <- environment(shapiro.test) > > More information should be available from the constants used in > swilk.c (in the R sources) which is a translation of the Fortran > code at http://lib.stat.cmu.edu/apstat/R94. > > Peter Ehlers > >> >> With thanks, >> Ted.
Thanks, Peter! Your suggestion environment(SW_ted) <- environment(shapiro.test) got it working, and now I can see that the variable 'sw' has six components, of which shapiro.test() returns only two: statistic = c(W = sw$w), p.value = sw$pw With X <- rnorm(30), SW <- SW_ted(X), one component, SW$init, is itself a list which includes a vector of 30 numbers which are the same as sort(X). Another is SW$a, a vector of 15 numbers: $a [1] 0.41668701 0.30027449 0.25057057 0.21649644 0.18856336 0.16440846 [7] 0.14279744 0.12299577 0.10452516 0.08705221 0.07033122 0.05417194 [13] 0.03842017 0.02294516 0.00763082 and I strongly suspect that these, in effect, will give the coefficients. Now to check it out! Thanks again, and best wishes, Ted. ------------------------------------------------- E-Mail: (Ted Harding) <ted.hard...@wlandres.net> Date: 05-Apr-2012 Time: 09:11:51 This message was sent by XFMail ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.