provides the following exercise related to pool-adjacent-violators algorithm.
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
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])
We need to learn first, and it turns out that this is exactly the simply ordered case mentioned in the wiki.
In , there is another exercise about isotonic regression, which is the continuation of the previous section.
