mypava <- function(y,m)
	{
	## make sure y and m are vectors of counts (no NA's same length
	## with y_i <= m_i and m_i >= 1
	theta <- y/m
	blocks <- 1:length(y)
	repeat
	{
		violators <- ( diff(theta) < 0 )
		if( !any(violators) ) break
		ii <- which.max( violators )
		b1 <- blocks[ii]
		b2 <- blocks[ii+1]
		indices <- (blocks == b1)|(blocks==b2)
		theta[indices] <- sum( y[indices] )/sum( m[indices] )
		blocks[indices] <- b1
	}
	theta
	}