-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of
Ted.Harding at nessie.mcc.ac.uk
Sent: Sunday, May 08, 2005 5:14 PM
To: Scott Briggs
Cc: r-help at stat.math.ethz.ch
Subject: RE: [R] Monotonic regression
On 08-May-05 Scott Briggs wrote:
Hi, I'm trying to find an implementation of monotonic
and I haven't been able to find anything that's really related to
this. isoMDS in the MASS package uses monotonic
I was wondering if there is any standalone function for monotonic
regression?
Basically what I'm trying to do is implement monotonic regression
where I can see not just the results of each iteration but also be
able to tweak the input in order to test for or "kick" the
out of a local minimum so that I can make sure I have the global
minimum.
Any help would be much appreciated. Thanks!
Scott
You may probably find PAVA ("Pool Adjacent Violators Algorithm")
useful. Below is code for a simple version which I have been using
for a few years. I forget where I found it!
An R site search comes up with code for a version with more
complex functionality at
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/9807.html
(contributed to r-help by Jan de Leeuw on 01 Jul 2004). I have
not tested this code.
pava<-function(x,wt=rep(1,length(x)))
{
n<-length(x)
if(n<=1) return(x)
if(any(is.na(x)) || any(is.na(wt))) {
stop("Missing values in 'x' or 'wt' not allowed")
}
lvlsets<-(1:n)
repeat {
viol<-(as.vector(diff(x))<0)
if(!(any(viol))) break
i<-min( (1:(n-1))[viol])
lvl1<-lvlsets[i]
lvl2<-lvlsets[i+1]
ilvl<-(lvlsets==lvl1 | lvlsets==lvl2)
x[ilvl]<-sum(x[ilvl]*wt[ilvl])/sum(wt[ilvl])
lvlsets[ilvl]<-lvl1
}
x
}
# Examples:
# > x<-c(1,0,1,0,0,1,0,1,1,0,1,0)
# > x
# [1] 1 0 1 0 0 1 0 1 1 0 1 0
# > pava(x)
# [1] 0.4 0.4 0.4 0.4 0.4 0.5 0.5 0.6 0.6 0.6 0.6 0.6
# >
# > pava(c(0,0,2/4,1/5,2/4,1/2,4/5,5/8,7/11,10/11),
# c(5,4,4,5,4,2,5,8,11,11))
# [1] 0.0000000 0.0000000 0.3333333 0.3333333 0.5000000
# [6] 0.5000000 0.6666667 0.6666667 0.6666667 0.9090909
# (example where data are {ri,ni} so x={ri/ni} and w={ni})
Best wishes,
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 08-May-05 Time: 20:56:28
------------------------------ XFMail ------------------------------