Misc
Pool-adjacent-violators algorithm
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.
1
pooladv <- function(f, w)
2
{
3
n = length(f)
4
lag = diff(f)
5
if (sum(lag < 0) == 0) # f is isotonic
6
return(f)
7
while (TRUE) {
8
idx = which(lag < 0)[1] + 1
9
newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
10
f[idx] = newf
11
f[idx-1] = newf
12
lag = diff(f)
13
if (sum(lag < 0) == 0)
14
return(f)
15
}
16
}
17
18
## example
19
pooladv(c(23, 27, 25, 28), rep(1, 4))
20
# ans: 23, 26, 26, 28
Copied!
We can also use another scientific programming language, Julia.
1
function pooladv(f::Array, w::Array)
2
n = size(f, 1)
3
lag = diff(f)
4
if all(i -> i >= 0, lag)
5
return(f)
6
end
7
while true
8
for i = 1:(n-1)
9
global idx
10
idx = i+1
11
lag[i] < 0 && break
12
end
13
newf = (w[idx]*f[idx] + w[idx-1]*f[idx-1])/(w[idx]+w[idx-1])
14
f[idx] = newf
15
f[idx-1] = newf
16
lag = diff(f)
17
if all(i -> i>=0, lag)
18
return(f)
19
end
20
end
21
end
22
23
## example
24
g = pooladv([23, 27, 25, 28], ones(4))
25
for i in eachindex(g)
26
println(g[i])
27
end
Copied!
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.
1
function treeordering(f::Array, w::Array)
2
n = size(f, 1)
3
lag = diff(f)
4
# if f is isotonic
5
if all(i -> i >= 0, lag)
6
return(f)
7
end
8
for j = 1:n
9
Aj = sum(f[1:j].*w[1:j])/sum(w[1:j])
10
if Aj < f[j+1]
11
return([Aj*ones(j); f[(j+1):end]])
12
end
13
end
14
end
15
16
## example
17
treeordering([18,17,12,21,16], [1,1,3,3,1])
Copied!
Last modified 3yr ago
Copy link