M
M
Monte-Carlo
Home
ESL-CN
TechNotes
Blog
Search…
Introduction
AIS
Generate R.V.
Gibbs
Markov Chain
MC Approximation
MC Integration
MC Optimization
MCMC
Metropolis-Hastings
PBMCMC
RJMCMC
Diagnosing Convergence
SMCTC
Tempering
Misc
R vs. Julia
References
Powered By
GitBook
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!
Previous
Parallel Tempering
Next
R vs. Julia
Last modified
3yr ago
Copy link
Contents
Pool-adjacent-violators algorithm
Tree ordering