[R] a maximazation question

[Jan de Leeuw]

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