Robert and Casella (2005) provides the following exercise related to pool-adjacent-violators algorithm.
We need to learn isotonic regression first, and it turns out that this is exactly the simply ordered case mentioned in the wiki.
Use R program to solve this problem, and the code is as follows.
pooladv <- function(f, w)
{
n = length(f)
lag = diff(f)
if (sum(lag < 0) == 0) # f is isotonic
return(f)
while (TRUE) {
idx = which(lag < 0)[1] + 1
newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
f[idx] = newf
f[idx-1] = newf
lag = diff(f)
if (sum(lag < 0) == 0)
return(f)
}
}
## example
pooladv(c(23, 27, 25, 28), rep(1, 4))
# ans: 23, 26, 26, 28
We can also use another scientific programming language, Julia.
function pooladv(f::Array, w::Array)
n = size(f, 1)
lag = diff(f)
if all(i -> i >= 0, lag)
return(f)
end
while true
for i = 1:(n-1)
global idx
idx = i+1
lag[i] < 0 && break
end
newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
f[idx] = newf
f[idx-1] = newf
lag = diff(f)
if all(i -> i>=0, lag)
return(f)
end
end
end
## example
g = pooladv([23, 27, 25, 28], ones(4))
for i in eachindex(g)
println(g[i])
end
Tree ordering
In Robert and Casella (2005), there is another exercise about isotonic regression, which is the continuation of the previous section.
The following Julia program can be used to solve this exercise.
function treeordering(f::Array, w::Array)
n = size(f, 1)
lag = diff(f)
# if f is isotonic
if all(i -> i >= 0, lag)
return(f)
end
for j = 1:n
Aj = sum(f[1:j].*w[1:j])/sum(w[1:j])
if Aj < f[j+1]
return([Aj*ones(j); f[(j+1):end]])
end
end
end
## example
treeordering([18,17,12,21,16], [1,1,3,3,1])
