Here is a version of pooled-adjacent-violators in R, which
does weighted mean, weighted median, and weighted p-fractile.

===============================================================

pava<-function(x,w=rep(1,length(x)),block=weighted.mean){
nblock<-n<-length(x); blocklist<-array(1:n,c(n,2)); blockvalues<-x; active<-1
repeat{
if (!is.up.satisfied(blockvalues,active)) {
blockmerge<-merge.block.up(blocklist,blockvalues,x,w,active,block)
blockvalues<-blockmerge$v; blocklist<-blockmerge$l
nblock<-nblock-1
while (!is.down.satisfied(blockvalues,active)) {
blockmerge<-merge.block.up(blocklist,blockvalues,x,w,active-1,block)
blockvalues<-blockmerge$v; blocklist<-blockmerge$l;
nblock<-nblock-1; active<-active-1;
}
}
else if (active == nblock) break() else active<-active+1
}
put.back(n,blocklist,blockvalues)
}

merge.block.up<-function(blocklist,blockvalues,x,w,i,block){
n<-length(blockvalues); nn<-1:n; ii<-which(i+1!=nn)
blocklist[i,]<-c(blocklist[i,1],blocklist[i+1,2])
indi<-blocklist[i,1]:blocklist[i+1,2]
blockvalues[i]<-block(x[indi],w[indi])
blocklist<-blocklist[ii,]
if (length(ii) == 1) dim(blocklist)<-c(1,2)
blockvalues<-blockvalues[ii]
list(v=blockvalues,l=blocklist)
}

put.back<-function(n,blocklist,blockvalues){
x<-rep(0,n);nb<-length(blockvalues)
for (i in 1:nb) {
x[blocklist[i,1]:blocklist[i,2]]<-blockvalues[i]}
return(x)
}

is.up.satisfied<-function(x,i) (i == length(x))||(x[i]<=x[i+1])

is.down.satisfied<-function(x,i) (i == 1)||(x[i-1]<=x[i])

weighted.median<-function(x,w=rep(1,length(x))){
ox<-order(x);x<-x[ox];w<-w[ox];k<-1
low<-cumsum(c(0,w)); up<-sum(w)-low; df<-low-up
repeat{
if (df[k] < 0) k<-k+1
else if (df[k] == 0) return((w[k]*x[k]+w[k-1]*x[k-1])/(w[k]+w[k-1]))
else return(x[k-1])
}
}

weighted.pth.fractile<-function(x,w=rep(1,length(x)),a=1,b=1){
ox<-order(x);x<-x[ox];w<-w[ox];k<-1
low<-cumsum(c(0,w)); up<-sum(w)-low; df<-a*low-b*up
repeat{
if (df[k] < 0) k<-k+1
else if (df[k] == 0) return((w[k]*x[k]+w[k-1]*x[k-1])/(w[k]+w[k-1]))
else return(x[k-1])
}
}

On Monday, December 23, 2002, at 08:15 AM, Thomas Lumley wrote:

On Sat, 21 Dec 2002, Shuangge Ma wrote:

Dear Sir/Madam:
this is shuangge Ma, graduate student in UW-Madison statistics department.
I have a computational question.
I have a function f(x,y). I want to find the y(x) that maximize f(x,y)
under the constraint y(x) is a non-decreasing step function.
Is there any R package or algorithm I can use for this purpose?
thanks a lot for your time and help,

Often for a problem like this a finite set of possible locations for the
steps are known (for a not-necessarily-unique solution) -- eg the observed
values of x. In that case the answer can be parametrised by the step
heights at each of these x values, with the constraint that these are
non-negative. The L-BFGS-B method of optim() will probably work.

If you don't know where the steps are, it's likely to be much harder.

One important special case that's worth mentioning is the isotonic
regression problem. If you have data (x_i,y_i) and are trying to fit an
increasing function by least squares or maximum likelihood the (or at
least a) solution is usually the isotonic regression, given by the
Pool-Adjacent-Violators Algorithm.


-thomas

______________________________________________
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help


===
Jan de Leeuw; Professor and Chair, UCLA Department of Statistics;
Editor: Journal of Multivariate Analysis, Journal of Statistical Software
US mail: 9432 Boelter Hall, Box 951554, Los Angeles, CA 90095-1554
phone (310)-825-9550; fax (310)-206-5658; email: [EMAIL PROTECTED]
homepage: http://gifi.stat.ucla.edu
------------------------------------------------------------------------ -------------------------
No matter where you go, there you are. --- Buckaroo Banzai
http://gifi.stat.ucla.edu/sounds/nomatter.au
------------------------------------------------------------------------ -------------------------

______________________________________________
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help


Reply via email to