Re: [R] what is the purpose of an error message in uniroot?
Matt, Some time back I didn't like the uniroot restriction either so I wrote a short function "manyroots" that breaks an interval into many shorter intervals and looks for a single root in each of them. This function is NOT guaranteed to find all roots in an interval even if you specify many subintervals. Graphing is always a good idea. There is always a chance the function dips to or below the axis and back up in an arbitrarily short interval. For example, f(x) = x^2. If one of your subintervals doesn't happen to end at 0, you're out of luck. (This is in contrast to uniroot working on a continuous function that is positive at one end of an interval and negative at the other end: at least one root is guaranteed and uniroot will find it given enough iterations.) Furthermore, if you are working with a polynomial, just use polyroot. Chris (Sorry for the dearth of code comments in the following) manyroots <- function(f, interval, ints=1, maxlen=NULL, lower = min(interval), upper = max(interval), tol = .Machine$double.eps^0.25, maxiter = 1000, ...) { if (!is.numeric(lower) || !is.numeric(upper) || lower >= upper) stop("lower < upper is not fulfilled") if (is.infinite(lower) || is.infinite(upper)) stop("Interval must have finite length") if (!is.null(maxlen)) ints <- ceiling((upper-lower)/maxlen) if (!is.numeric(ints) || length(ints)>1 || floor(ints)!=ints || ints<1) stop("ints must be positive integer") ends <- seq(lower, upper, length=ints+1) fends <- numeric(length(ends)) for (i in seq(along=ends)) fends[i] <- f(ends[i], ...) zeros <- iters <- prec <- rep(NA, ints) for (i in seq(ints)) { cat(i, ends[i], ends[i+1], fends[i], fends[i+1], "\n") if (fends[i] * fends[i+1] > 0) { #cat("f() values at end points not of opposite sign\n") next; } if (fends[i] == 0 & i>1) { #cat("this was found in previous iteration\n") next; } val <- .Internal(zeroin(function(arg) f(arg, ...), ends[i], ends[i+1], tol, as.integer(maxiter))) if (as.integer(val[2]) == maxiter) { warning("Iteration limit (", maxiter, ") reached in interval (", ends[i], ",", ends[i+1], ").") } zeros[i] <- val[1] iters[i] <- val[2] prec[i] <- val[3] } zeros <- as.vector(na.omit(zeros)) fzeros <- numeric(length(zeros)) for (i in seq(along=zeros)) fzeros[i] <- f(zeros[i], ...) list(root = zeros, f.root = fzeros, iter = as.vector(na.omit(iters)), estim.prec = as.vector(na.omit(prec))) } gg <- function(x) x*(x-1)*(x+1) manyroots(gg, c(-4,4), 13, maxiter=200, tol=10^-10) hh <- function(x,x2) x^2-x2 manyroots(hh, c(-10, 10), maxlen=.178, x2=9) manyroots(sin, c(-4,20), maxlen=.01) #but ss <- function(x) sin(x)^2 manyroots(ss, c(-4,20), maxlen=.01) plot(ss, -4,20) abline(h=0) From: [EMAIL PROTECTED] [EMAIL PROTECTED] wrote: > > This is probably a blindingly obvious question: Yes, it is. > > Why does it matter in the uniroot function whether the f() values at > > the end points that you supply are of the same sign? Plot some graphs. Think about the *name* of the function --- *uni*root. Does that ring any bells? And how do you know there *is* a root in the interval in question? Try your ``uniroot2'' on f(x) = 1+x2 and the interval [-5,5]. To belabour the point --- if the f() values are of the same sign, then there are 0, or 2, or 4, or roots in the interval in question. Rolf, Only if f is continuous (of course finding roots of discontinuous functions is a greater challenge) The ***only chance*** you have of there being a unique root is if the f() values are of opposite sign. The algorithm used and the precision estimates returned presumably depend on the change of sign. You can get answers --- sometimes --- if the change of sign is not present, but the results could be seriously misleading. Without the opposite sign requirement the user will often wind up trying to do something impossible or getting results about which he/she is deluded. cheers, Rolf Turner [EMAIL PROTECTED] P. S. If the f() values are of the same sign, uniroot() DOES NOT give a warning! It gives an error. R. T. -- Christopher Andrews, PhD SUNY Buffalo, Department of Biostatistics 242 Farber Hall, [EMAIL PROTECTED], 716 829 2756 __ R-help@stat.math.ethz.ch 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 prov
Re: [R] Switching labels on a factor
Bob, This is I think exactly what one wants to have happen. The first four observations are still women. Both the labels and the underlying integers should change. (If you want to give all the people sex changes, try Relevel in the Epi package. mydata$afterthechange <- Relevel(mydata$gender, list(m="f", f="m")) mydata workshop gender q1 q2 q3 q4 gR afterthechange 11 f 1 1 5 1 f m 22 f 2 1 4 1 f m 31 f 2 2 4 3 f m 42 f 3 1 NA 3 f m 51 m 4 5 2 4 m f 62 m 5 4 5 5 m f 71 m 5 3 4 4 m f 82 m 4 5 5 NA m f unclass(mydata$afterthechange) [1] 1 1 1 1 2 2 2 2 attr(,"levels") [1] "m" "f" Chris Date: Fri, 15 Dec 2006 15:34:15 -0500 From: "Muenchen, Robert A (Bob)" <[EMAIL PROTECTED]> Subject: [R] Switching labels on a factor To: Message-ID: <[EMAIL PROTECTED]> Content-Type: text/plain; charset="US-ASCII" Hi All, I'm perplexed by the way the unclass function displays a factor whose labels have been swapped with the relevel function. I realize it won't affect any results and that the relevel did nothing useful in this particular case. I'm just doing it to learn ways to manipulate factors. The display of unclass leaves me feeling that the relevel had failed. I've checked three books & searched R-help, but found no mention of this particular issue. The program below demonstrates the problem. Is this a bug, or is there a reason for it to work this way? Thanks, Bob mystring<- ("id,workshop,gender,q1,q2,q3,q4 1,1,f,1,1,5,1 2,2,f,2,1,4,1 3,1,f,2,2,4,3 4,2,f,3,1, ,3 5,1,m,4,5,2,4 6,2,m,5,4,5,5 7,1,m,5,3,4,4 8,2,m,4,5,5,9") mydata<-read.table(textConnection(mystring), header=TRUE,sep=",",row.names="id",na.strings="9") mydata # Create a gender Releveled variable, gR. # Now 1=m, 2=f mydata$gR <- relevel(mydata$gender, "m") # Print the data to show that the labels of gR match those of gender. mydata # Show that the underlying codes have indeed reversed. as.numeric(mydata$gender) as.numeric(mydata$gR) # Unclass the two variables to see that print order # implies that both the codes and labels have # flipped, cancelling each other out. For gR, # m appears to be associated with 2, and f with 1 unclass(mydata$gender) unclass(mydata$gR) = Bob Muenchen (pronounced Min'-chen), Manager Statistical Consulting Center U of TN Office of Information Technology 200 Stokely Management Center, Knoxville, TN 37996-0520 Voice: (865) 974-5230 FAX: (865) 974-4810 Email: [EMAIL PROTECTED] Web: http://oit.utk.edu/scc, News: http://listserv.utk.edu/archives/statnews.html -- Christopher Andrews, PhD SUNY Buffalo, Department of Biostatistics 242 Farber Hall, [EMAIL PROTECTED], 716 829 2756 __ R-help@stat.math.ethz.ch 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.
[R] Biobass, SAGx, and Jonckheere-Terpstra test
Horace, Last year I found another Jonckheere-Terpstra function at http://www.biostat.wustl.edu/archives/html/s-news/2000-10/msg00126.html Using that and the version in SAGx and wanting a few other options, I created the version below. Hope it is useful, Chris -- Christopher Andrews, PhD SUNY Buffalo, Department of Biostatistics # # # j.test <- function(x, y, alternative = c("two.sided", "decreasing", "increasing"), asymp=TRUE, correct=FALSE, perm=0, na.action=c("omit", "fail"), permgraph=FALSE, permreps=FALSE, ...) { #Jonckheere-Terpstra test #x = response #y = group membership #alternative hypothesis #asymptotic = asymptotic formula for variance or don't bother #correct = continuity correction #perm = number of repetitions for a permutation test #na.action = what to do with missing data #permgraph = draw a histogram of the permutations? #permreps = return the permutations? #... for graph alternative <- match.arg(alternative) complete <- complete.cases(x,y) nn <- sum(complete) cat(nn, "complete cases.\n") if (!all(complete) && (na.action=="fail")) { cat("option na.action is 'fail' and", sum(!complete), "cases are not complete.\n") return(NULL) } response <- x[complete] groupvec <- y[complete, drop=TRUE] if(!is.ordered(groupvec)) { groupvec <- as.ordered(groupvec) cat("y was not an ordered factor. Redefined to be one.\n") } lev <- levels(groupvec) nlev <- length(lev) if(nlev <= 2) { if (nlev < 2) { cat("Not enough groups (k =", nlev, ") for Jonckheere.\n") return(NULL) } else { cat("Two groups. Could use wilcox.test.\n") } } computestat <- function(response, groupvec, n.groups) { pairwise <- function(x, y) { length(x)*(length(y) + (1+length(x))/2) - sum(rank(c(x,y))[seq(along=x)]) } H<-0 for (i in seq(n.groups-1)) for (j in seq(i+1, n.groups)) H <- H+pairwise(response[groupvec == lev[i]], response[groupvec == lev[j]]) return(H) } retval <- list(statistic=computestat(response, groupvec, nlev), alternative = paste(alternative, paste(levels(groupvec),collapse=switch(alternative, two.sided = ", ", decreasing= " > ", increasing = " < ")), sep=": "), method = "Jonckheere", data.name = paste(deparse(substitute(x)), "by", deparse(substitute(y if (asymp) { ns <- tabulate(as.numeric(groupvec)) retval$EH <- (nn * nn - sum(ns * ns))/4 retval$VH <- if (!any(duplicated(response))) { (nn * nn * (2 * nn + 3) - sum(ns * ns * (2 * ns + 3)))/72 } else { ds <- as.vector(table(response)) ((nn*(nn-1)*(2*nn+5) - sum(ns*(ns-1)*(2*ns+5)) - sum(ds*(ds-1)*(2*ds+5)))/72 + sum(ns*(ns-1)*(ns-2))*sum(ds*(ds-1)*(ds-2))/(36*nn*(nn-1)*(nn-2)) + sum(ns*(ns-1))*sum(ds*(ds-1))/(8*nn*(nn - 1))) } pp <- if (!correct) { pnorm(retval$statistic, retval$EH, sqrt(retval$VH)) } else if (retval$statistic >= retval$EH + 1) { pnorm(retval$statistic - 1, retval$EH, sqrt(retval$VH)) } else if (retval$statistic <= retval$EH - 1) { pnorm(retval$statistic + 1, retval$EH, sqrt(retval$VH)) } else { .5 } retval$p.value <- switch(alternative, two.sided = 2*min(pp,1-pp), decreasing = pp, increasing = 1-pp) } if (perm>0) { reps <- numeric(perm) for (i in seq(perm)) { reps[i] <- computestat(response, sample(groupvec), nlev) } retval$EH.perm <- mean(reps) retval$VH.perm <- var(reps) ppp <- if (!correct) { pnorm(retval$statistic, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic >= retval$EH.perm + 1) { pnorm(retval$statistic - 1, retval$EH.perm, sqrt(retval$VH.perm)) } else if (retval$statistic <= retval$EH.perm - 1) { pnorm(retval$statistic + 1, retval$EH.perm, sqrt(retval$VH.perm)) } else { .5 } retval$p.value.perm <- switch(alternative, two.sided = 2*min(ppp,1-ppp), decreasing = ppp, increasing = 1-ppp) rr <- rank(c(retval$statistic, reps))[1]/(perm+2) retval$p.value.rank <- switch(alternative, two.sided = 2*min(rr,1-rr), decreasing = rr, increasing = 1-rr) if (permreps) retval$reps <- reps if (permgraph) { hist(reps, xlab="Jonckheere Statistic", prob=TRUE, main="Histogram of Permutations", sub=paste(perm, "permutations")) abline(v=retval$statistic, col="red") points((-3:3)*sqrt(retval$VH.perm)+retval$EH.perm, rep(0,7), col="blue", pch=as.character(c(3:1,0:3))) } } class(retval) <- "htest" names(retval$statistic) <- "J" return(retval) #returns list (of class htest) with the following components #Always: #statistic = value of J #alternative = same as input #method = "Jonckheere" #data.name = source of data based on call #Most of the time (i.e., when asymp=TRUE): #EH = expected value based on sample
Re: [R] Finding all possible partitions of N units into k classes
nkpartitions <- function(n, k, exact=FALSE, print=FALSE) { # n objects # k subgroups # exactly k or at most k? # print results as they are found? if (n != floor(n) | n<=0) stop("n must be positive integer") if (k != floor(k) | k<=0) stop("k must be positive integer") if (print) { printnkp <- function(a) { for (j in seq(max(a))) cat("{", seq(along=a)[a==j], "} "); cat("\n") } } # How many? Stirling2nd <- function(n, k) { sum((-1)^seq(0,k-1) * choose(k, seq(0,k-1)) * (k-seq(0,k-1))^n) / factorial(k) } rows <- Stirling2nd(n,k) if (!exact & k>1) { for (i in seq(k-1,1)) { rows <- rows + Stirling2nd(n,i) } } if (print) cat("rows =",rows,"\n") # Allocate space theparts <- matrix(NA, nrow=rows, ncol=n) # begin counting howmany <- 0 # all in one group a <- rep(1,n) # does this count? if (!exact | (k==1)) { # increase count, store, and print howmany <- howmany + 1 theparts[howmany,] <- a if (print) printnkp(a) } # search for others repeat { # start at high end last <- n repeat { # increment it if possible if ((a[last] <= max(a[1:(last-1)])) & (a[last] < k)) { a[last] <- a[last]+1 # does this count? if (!exact | max(a)==k) { # increase count, store, and print howmany <- howmany + 1 theparts[howmany,] <- a if (print) printnkp(a) } # start again at high end. break } # otherwise set to 1 and move to a different object a[last] <- 1 if (last>2) { last <- last-1 next } # report the partitions return(theparts) } } } nkpartitions(5,3) nkpartitions(5,3,T) -- Dr Christopher Andrews University at Buffalo, Department of Biostatistics 242 Farber Hall [EMAIL PROTECTED] 716 829 2756 __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html