a maximazation question

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

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 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

______________________________________________
R-help at stat.math.ethz.ch 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: deleeuw at stat.ucla.edu
homepage: http://gifi.stat.ucla.edu
   
------------------------------------------------------------------------ 
-------------------------
           No matter where you go, there you are. --- Buckaroo Banzai
                    http://gifi.stat.ucla.edu/sounds/nomatter.au
   
------------------------------------------------------------------------ 
-------------------------