頻度論とベイズをつなぐ統計的信頼度
下平英寿
大阪大学 基礎工学研究科
理化学研究所 革新知能統合研究センター
２０１６年１１月１６日
第19回情報論的学習理論ワークショップ (IBIS2016)1

最近の様子@下平研2

本日のトーク： 地味にやってる話• Multiple comparisons of log-likelihoods with applications tophylogenetic inference. H Shimodaira, M Hasegawa. Molecularbiology and evolution 16, 1114-1116, 1999
• An approximately unbiased test of phylogenetic tree selection. H
Shimodaira. Systematic biology 51 (3), 492-508, 2002
• Approximately unbiased tests of regions using multistep-multiscalebootstrap resampling. H Shimodaira. The Annals of Statistics 32,
2004• Testing regions with nonsmooth boundaries via multiscalebootstrap. H Shimodaira. Journal of Statistical Planning andInference 138, 1227-1241, 2008
• Higher-order accuracy of multiscale-double bootstrap for testingregions. H Shimodaira. Journal of Multivariate Analysis 130, 208223, 2014
マルチスケール・ブートストラップ法による信頼度計算の研究
3

領域の検定と信頼度

𝜇𝑌|𝜇 データの分布𝐻𝑦confidence 𝜇 ∈ 𝐻 𝑦
信頼度が小さい時，仮説を棄却する

4

多変量正規モデルで考える𝜇𝐻

𝑌

𝑌|𝜇~𝑁 𝜇, 𝐼
density ∝ exp(− y − µ / 2)
2

5

ベイズ事後確率は容易に計算
𝑌|𝜇~𝑁 𝜇, 𝐼

µ ~ uniformµHydensity ∝ exp(− y − µ / 2)2

µ | y ~ N(y, I )muの事後分布 ＝ ブートストラップ分布

confidence23456 𝜇 ∈ 𝐻 𝑦 = 𝑃 𝜇 ∈ 𝐻|𝑦
Hの事後確率＝ ブートストラップ確率6頻度論のp-値が満たすべき性質
𝑌𝜇Confidence = 0.05

𝜕𝐻

𝐻
Type-I errorを有意水準(=0.05)で抑えたい．しかし検出力を上げるには等号に近い方が良い．

P confidence:;5< 𝜇 ∈ 𝐻 𝑌 < 0.05 ≤ 0.05 under 𝜇 ∈ 𝐻
不偏検定では，仮説境界上でType-I errorを有意水準と同じにしたい

P confidence:;5< 𝜇 ∈ 𝐻 𝑌 < 0.05 = 0.05 under 𝜇 ∈ 𝜕𝐻
7

符号付き距離を使う

H

µ̂

d(y)
y

符号付き距離(signed distance)

Y | y ~ N( µ̂, I )confidence:;5< 𝜇 ∈ 𝐻 𝑦 = 𝑃 𝑑 𝑌 > 𝑑 𝑦 |𝜇̂L
KM+ 𝑂H (𝑛 )8もしパラメータにアクセスすれば簡単
ブートストラップ分布をひっくり返す

データからブートストラップ分布

H : q £ q0

BP(1)
qˆ*

q0

qˆ

帰無仮説からシミュレーション
p-value
qˆ*

q0

qˆ
9

ベイズ事後確率はp-値を近似する
多変量正規モデル，平坦な事前分布，仮説領域の境界が滑らかであることを仮定すると
O
KM
confidence:;5< 𝜇 ∈ 𝐻 𝑌 = confidence23456 𝜇 ∈ 𝐻 𝑌 + 𝑂H (𝑛 )たしかに近似してるけど，誤差が大きい．誤差の主な原因は，仮説境界の「平均曲率」．
頻度論の意味でp-値を追求してみる．
ブートストラップ法で計算したい．カウント値だけをつかう．
多重比較との関連．

10たとえば、モデル選択• 多項式回帰モデルの次数選択• 適切なモデルをAIC最小化法で選ぶ

y = β0 + β1x + β 2 x + ε211各領域がモデルに相当
多項式回帰の次数選択： 次数＝０，１，２から選ぶ

y = β0 + ε

y = β0 + β1x + β 2 x 2 + ε

y = β0 + β1x + ε( β̂1 , β̂ 2 )Efron and Tibshirani (1998)“The problem of regions”

12たとえば、進化系統樹の推定• DNA配列データ• 系統樹上のマルコフ過程• 最尤法でモデル選択する

13

各領域が系統樹に相当
1

1

2

3

4

3

2

4

5

5
1

2

3

4

5

1

Dn

2

3

4

5

D¥
14

Pvclust: Suzuki and Shimodaira (2006)
PvclustCluster dendrogram with AU/BP values (%)

fetal_lung

X165.96_Adeno
X178.96_Adeno
X185.96_Adeno

X11.00_Adeno
X181.96_Adeno
X156.96_Adeno
X198.96_Adeno
X180.96_Adeno
X187.96_Adeno
X68.96_Adeno
X137.96_Adeno
X12.00_Adeno
X199.97_Adeno_c
X199.97_Adeno_p
X204.97_Adeno
X257.97_Adeno
X132.95_Adeno
X319.00PT_Adeno
X320.00_Adeno_p
X320.00_Adeno_c
X207.97_SCLC
X314.99_SCLC
X230.97_SCLC
X315.99_SCLC
X315.99_node
X161.96_Adeno
X299.99_Adeno
X75.95_combined
X69.96_Adeno
X157.96_SCC
X219.97_SCC
X166.96_SCC
X3_SCC
X58.95_SCC
X232.97_SCC
X232.97_node
X220.97_SCC
X220.97_node
X42.95_SCC
X239.97_SCC
X245.97_SCC
X245.97_node
X246.97_SCC_p
X246.97_SCC_c
X237.97_Adeno
X147.96_Adeno
X6.00_LCLC
X256.97_LCLC
X248.97_LCLC
X191.96_Adeno
X139.97_LCLC
X59.96_SCC
X234.97_Adeno
X184.96_Adeno
X184.96_node
X218.97_Adeno
X319.00MT1_Adeno
X319.00MT2_Adeno
X265.98_Adeno
X80.96_Adeno
X223.97_Adeno

85 26
74 29
71
77 30
70
83 27
69
70 4
68
79 16
67
69 3
66
65 72 16
84 58
99 96
73 59
68 21
64
62 22
71 44 62
61 93 6860
100 9988 5676 6263 59
58
57
56
71 37
55
71 39
54
95 95
80 17
66 45
61 2753
91 77
52
75 21
51
99 97
50
67
31
49
47
59
100 100
4887 584793 7946
45
100 100
63 53
44
43
42
41 99 93
40
39 91 79
98 94
38
96 97
92 86
37
99 96
100 100
36
97 89
35
34
33
87755930
32
31
83
60
30
75 47
96 93
96 92
2927
84 49
92
95
28
100
100
26
25
24
100 10021
96 90
23
22
10061 10042
100 10020
19
1817
100 100
99 99
16
100 100
100 100
15
14
100 100
100 100 100 100
13
12
11
10
10079 100966
100 100
100 100
100 100
86
5 7 100 100
99 100 4
3
2
100 100

X306.99_normal
X219.97_normal
X222.97_normal
X314.99_normal
X315.99_normal
X222.97_Adeno
X226.97_Adeno
X306.99_Adeno
X306.99_node
X313.99PT_Adeno
X313.99MT_Adeno

0.6
0.0

0.2

0.4

Height

0.8

1.0

au bp
edge #

1

Distance: correlation
Cluster method: averageFig. 1. Hierarchical clustering of 73 lung tumors. The data are expression pattern of 916 genes of Garber et al. (2001). Values at branches are AU p-values (left),
BP values (right), and cluster labels (bottom). Clusters with AU % 95 are indicated by the rectangles. The fourth rectangle from the right is a cluster labeled 62 with
AU ¼ 0.99 and BP ¼ 0.96.

Function seplot provides a graphical interface for examining
standard errors, while print gives more detailed information15

1,2,3

4

1,2

cag-GFP

doi:10.1038/nature12969

5

7

6

6

3

cag-GFPWe recently discovered an unexpected phenomenon of somatic cell
reprogramming into pluripotent cells by exposure to sublethal stimwhich cells
we callinvolves
stimulus-triggeredof pluripotency
derivation of STAP stem cells fromuli,STAP
the lossacquisitionof(STAP)1. This reprogramming does not require nuclear transfer2,3
competence to differentiate into placental
lineages.
Consistent
withor genetic manipulation4. Here we report that reprogrammed STAP
this idea, STAP stem cells show little cells,
expressionof trophoblast
markerunlike embryonicstem (ES) cells,
can contribute to both embryonic and placental tissues, as seen in a blastocyst injection assay.
genes (Fig. 1f).
Mouse STAP cells lose the ability to contribute to the placenta as
We next examined whether an alterationin culturemarker
conditions
could
well as trophoblast
expression
on converting into ES-likeby treatment
with adrenocorticotropic
hormone (ACTH)
induce in vitro conversion of STAP stem
cellscells
into
cells similar
to trophoand leukaemia inhibitory factor (LIF). In contrast, when cultured
fromblastocysts
duringproblast stem cells8,9, which can be derived
with Fgf4, STAP cells give rise to proliferative stem cells with enhanced
longed adhesion culture in the presenceof Fgf4.
When wecultured
trophoblastic
characteristics.
Notably,
unlike conventional trophoblast stem
cells,
theoneFgf4-induced
stemwellcells from STAP cells conSTAP cell clusters under similar conditions
(Fig.
2a;
cluster per
tribute to both embryonic and placental tissues in vivo and transformin a 96-well plate), flat cell colonies grewout bycells
days
7–10
(Fig.
2b,
left;
into ES-likewhencultured
with
LIF-containing
medium. Taken

4

2

2,8

cag-GFP

together, the developmental potential of STAP cells, shown by chimaera formation and in vitro cell conversion, indicates that they
represent a unique state of pluripotency.
We recently discovered an intriguing phenomenon of cellular fateconversion: somatic cells regain pluripotency after experiencing sublethal stimuli such as a low-pH exposure1. When splenic CD451 lymphocytes are exposed to pH 5.7 for 30 min and subsequently cultured
in the presence of LIF, a substantial portion of surviving cells start toexpress the pluripotent cell marker Oct4 (also called Pou5f1) at day 2.
1.2 cell clusters form with a bona fide pluripotency
By day 7, pluripotent
marker profile and acquire the competence for three-germ-layer differentiation as shown
1.0by teratoma formation. These STAP cells can also
efficiently contribute to chimaeric mice and undergo germline transmission using 0.8
a blastocyst injection assay1. Although these characteristics resemble those of ES cells, STAP cells seem to differ from ES
cells in their limited
0.6capacity for self-renewal (typically, for only a few

12
10
8
6
4
2
0

FI
-S
C
FI
-S 1
C
FI
2
G -SC
FP 3
G ES
FP 1
G -ES
FP 2
-

LETTER

Percentage ofed efficiently in the presence of Fgf4 by passaging for moreassages with trypsin digestion every third day. Hereafter,
iferative cells induced from STAP cells by Fgf4 treatment
Bright-field
DAPI
d to as Fgf4-induced stem cells. This type of derivation into
st-stem-like cells is not common with ES cells (unless genetBidirectional developmental potential innipulated)13 or STAP stem cells.
blastocyst injection assay, unlike
STAP stem cells, thepla-with acquired pluripotency
reprogrammed
cells
tribution of Fgf4-induced stem
cells
(cag-GFP-labelled)
was, Mitsutaka Kadota , Munazah Andrabi , Nozomu Takata , Mikiko Tokoro ,HarukoObokata, Yoshiki Sasai , Hitoshi Niwa
, Shigenobu Yonemura , Charles A. Vacanti & Teruhiko Wakayama
with 53% of embryos (Fig. 2f,Yukarig; nTerashita
5LETTER
60).
InRESEARCH
the chimaeric

1

2

3

4

5

Embryo-specific
Embryo+placenta+yolk sac

0.8
0.6
0.4
0.2

0

ST

2

3

ES STAP CD45 TS

Eomes

Elf5

1.6
1.4
1.2
1.0
0.8
0.6
0.4
0.2

0

TS

STAP CD45 ES

1.4
1.2
1.0

Cdx2

0.8

Eomes

0.6

Elf5

0.4
0.2

ES

TS

0

ST
A

C

3

1

2

TS

C

-S

AP

ST

ST

1

Cdx2

Relative expression levels(TS=1.0)

Relative expression levels(ES=1.0)

1.0

Relative expression levels(TS=1.0)

20
18
16
14
12
10
8
6
4
2
0

AP

1

1.2

PSC

Long exposureC

cag-GFP

Rex1

1.4

ES
STAP
n = 50 n = 10-SBright-field

Nanog

Oct4

1.6

Placental contribution

Percentage of GFP+ cells

STAP
STAPchimaera
chimera

100
90
80
70
60
50
40
30
20
10
0

-S

Long exposureAP

cag-GFP

Percentage of embryos witheach contribution pattern

Bright-field

4

Laboratory for Cellular Reprogramming, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan. 2Laboratory for Genomic Reprogramming, RIKEN Center for Developmental Biology, Kobe 6500047, Japan. 3Laboratory for Tissue Engineering and Regenerative Medicine, Brigham and Women’s Hospital, Harvard Medical School, Boston, Massachusetts 02115, USA. 4Laboratory for Organogenesis
and Neurogenesis, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan. 5Laboratory for Pluripotent Stem Cell Studies, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan.
6
Genome Resource and Analysis Unit, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan. 7Electron Microscopy Laboratory, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan.
8
Faculty of Life and Environmental Sciences, University of Yamanashi, Yamanashi 400-8510, Japan.

Relative expression levels(TS=1.0)

Relative expression levels(ES=1.0)

0.10
0.05

Contribution pattern

ES chimaera

0.00

Height0.15

Figure 2 | Fgf4 treatment induces someaucharacter in STAP cells.
i trophoblast-lineagejka, Schematic of Fgf4 treatment to induce
Fgf4-induced
stem cells from STAP cells. b, Fgf4
ng
CD45
***
1.2
NS
ytreatment promoted the generation of flat cellme1.0
clusters that expressed Oct4-GFP at moderateay
0.8levels(right).
Top
andmiddle:
days
1
and
7
of
NS
100
lls
he
0.6
culture with Fgf4, respectively. Bottom: culture
ly, typically in ,30%
of wells).
The Fgf4-induced
cells bar,strongly
expressed
after
the first
passage.
Scale
50
mm.
0.4
9–12
0.4
STAP
he the trophoblast marker proteinsintegrin100
a7 (Itga7) and eomesoderc,d,Immunostainingof
Fgf4-induced
cells
with
0.2
in min (Eomes) (Fig. 2c, d) and marker genes (for example, Cdx2; Fig. 2e).
0.2100integrin
the trophoblast
stem
cell
markersa
These Fgf4-induced
cells with trophoblast
markerexpression
could a7
0
0 d
c
***
lls be expanded efficiently
the presence of
Fgf4
by passaging
for50
more
(c) andineomesodermin(d).
Scale bar,
mm.
TS
TS
JAKi+ JAK
JAKi– JAKi+ JAKi– JAKi+
100
cy than 30 passages
with trypsin
digestion
every thirdday. Hereafter,e, qPCR
analysisof marker
expression.
FI-SC
1d these proliferative cells induced from STAP cells FI-SC
by Fgf4 treatment
ESFI-SCf, g, Placental contribution of
Fgf4-induced stem
re, are referred to as Fgf4-induced stem cells. This type of derivation into
ES
STAP-SC
cells (FI-SCs)labelledwithconstitutive
hat trophoblast-stem-likecells is not (geneticallycommon with ES
cells (unlessgenetEomes
Nanog
Cdx2
Rex1
Oct4
Elf5
13
lls ically manipulated)
or
STAP
stem
cells.
GFP expression).b Scale bars: 5.0 mm (f (left panel)ef
In the blastocyst
assay,
stemIncells,
the pla-toandinjectiong); 50 mm(f,unlikerightSTAP
panel).
additions
to cental contribution of Fgf4-induced stem cells (cag-GFP-labelled) was
placental contribution, Fgf4-induced stem cells
he observed with 53% of embryos (Fig. 2f, g; n 5 60). In the chimaeric
contributed to the embryonic portion at a
Extended Data Figure
6 | Characterizationof ES-like cells
converted from
Oct4-gfp Fgf4-induced stem cells in trophoblast stem-cell medium containing
moderatelevel (g). h, Quantificationof placental
Extended
Data
Figure 6 | Characterizationof colonies
ES-likecellswasconvertedfrom from
Oct4-gfp Fgf4-inFigure2|
Fgf4
treatment
inducessome
Macmillan
Publishers
Limited.
Allrightsres©2014
observed
Fgf4-induced stem
cells and comparisonof STAP
cells withearly embryos. Fgf4.
No substantial
formationofOct4-GFP1
contribution
by FACScharacter
analysis.
UnlikeFgf4-induced stem cells and comparison of STAP cells with early embryos. Fgf4. No substan
trophoblast-lineage
in STAP
cells.Fgf4dissociateda,Fgf4-induced
stem ofcells
in cells
MEK-inhibitor-containinga, Immunohistochemistryof ES-likecells
fornottrophoblastanddissociated Fgf4ImmunohistochemistryES-likefor trophoblast and pluripotencymedium.
Schematic
Fgf4
treatment
tocontributeinduceinduceda, cells,
ESofcells
did
topluripotency
Data
Figure6|
CharacterizationofES-likecells
fromFgf4-indu
Scalebar, 100 mm
markers.
ES-likecellsforconverted
convertedfrom Fgf4-induced
stem cellsOct4-gfp
noclustering
longer
Fgf4-induced
stem
cells
STAP cells.
b, embryonic
Fgf4
analysisof FACS-sorted
Oct4-GFP-strong
STAP
cells
pluripotent
markerFigure1Fgf4-induced| STAPfrom
cells
contributeto bothand longer
placental tissues
Scale
bar,
100
mm.
d,
Clustertree
diagramfromhierarchicalof
globalmarkers. ES-likeplacental
cellsExtended
converted
from
stem
cells
no
tissues
atvivo.
a detectablei, Cluster
ina, theb, E12.5
embryoslevel.
frominjected withtree
ES cells (a) andgenes (left) and trophoblast
marker genes(right).
Values are shownas ratio
toexpression profi
expressedthetrophoblast
marker(integrinalpha 7), but they did express thetreatment
promoted
generationof blastocysts
flat cell
expression
profiles.
Red,
AU
P
values.
expressed
the trophoblast
marker(integrin
alpha
but
they
did
STAP
cells (b). Bothcells are 7),
geneticallylabelled
with
GFPexpress
driven by a the the expression
level in ES cells.
Error bars represent
s.d. e, Contribution
to As this analysis included morula andgf4
medium
Fgf4-induced
stem
cells
and
comparisonof STAPcells
with(Oct4,
earlyembryos. Fgf4. No substantia
blastocyst embry
pluripotency
diagramfromthathierarchicalclustering
clusters
expressed
Oct4-GFP
atmoderate
constitutive promoter.
Progenyof
STAP ofcells global
also contributed to placental
placental tissues. Unlike parental
STAP cellsmarkers
and trophoblast
stemNanog(TS) cells,and SSEA-1). Scale bar, 100 mm.
ontribution
embryos
from
which
only
small fromamountsof RNApluripotency markers (Oct4,
Nanog
and
SSEA-1).
Scalebar,
100 mm.
we used pre-amp
b, Pluripotency
ES-likecells
converted
Fgf4-induced
stemcould
cells asbe obtained,tissuesandfetalmembranes(b),whereasES-cell-derived
cells were not found STAP blastocyst
stem cells (STAP-SCs)
did not retain
the ofabilityfor placental
levels(right).
Top
andmiddle:
days
1
and
7
of
expression
profiles.
approximately
unbiased
Eomes
in theseRed,tissues (a).
Scale bar, 5.0 mm. c, PercentagesofES-likefetuses in which injected
contributions.
Three
independent
testedandall showed
substantial
dissociatedFgf4-ina, Immunohistochemistryfortrophoblastandpluripotency
Sequencing(Closhownlines
by wereteratoma
formation.
Thosecells successfully
formedteratomas
we usedpre-amplification
withthe
SMARTer
Ultra LowRNA
kit for Illuminab, Pluripotency of ES-likecellswithconvertedfrom
Fgf4-induced
cells ascells
cultureFgf4,
respectively.
Bottom:cultureofstemcells contributed only to the embryonic portion (red) or also to placental
contributions to the embryonic portions. f, qPCR analysis of trophoblast
P values.
j, theqPCRanalysisofcells
profile of STAP
tissuesfrom
all threelayers: neuroepithelium (left, arrow
after
first
passage.
bar,Fgf4-induced
50
mm.
and
yolkcellssac Scaletissues(blue).
***P, 0.001with Fisher’steratomas
exact test. d, qPCR
markerSequencing
gene expression incontaining
STAP
stem cells.
Error
bars
represent
s.d. germ(ClontechLaboratories).
e,f,Volcanoplotoftheexpression
shown by teratoma formation.
Those
successfully
formedScale bar, 100
mm.
markers.
ES-likecells
converted
Fgf4-induced
stem
cells
no longer
showing
greaterd
indicates),muscletissue (middle,
arrow indicates) and bronchial-likec, d,under
Immunostainingof Fgf4-induced
cells with(culturedfeeder-free
conditions)with(left,or from
profile
of
STAP
cells
comparedtothemorula(e)andblastocyst
(f).
Genes
containing tissues
from allthreegermlayers:neuroepithelium
arrow
red
and
are
cons
epithelium(right).
Scale
bar,
100
mm.
c,
MEK
inhibitor
treatment
assay
for
the trophoblast stem cell markers integrin a7
16profiles
26
without
JAK
inhibitor(JAKi)
treatment
forexpressed
the indicates)
trophoblast
marker(integrinshowing
alphagreater7), but
express
the1.0 3 10expression
are highlighted
inthan they
10-folddid
changeand P value
indicates), muscle
tissue(middle,arrow
andbronchial-like(c)and eomesodermin(d). Scale bar,
50 mm.
pluripotent
markerqPCR(Oct4,
analysisof assay fore,bar,
qPCR
analysisexpression.
red and are considered
up-100(or down-)
cells.
epithelium (right).
Scale
100
mm.ofgenes.
c,marker
MEKk,
inhibitor
treatment
blastocyst
embryos
pluripotency
markers
Nanog
and SSEA-1).
Scale bar,
mm. regulated in the STAP

en
ta
hat

Percentage of GFP+ cell

20
18
16
14
12
10
8
6
4
2
0

LETTER RESEARCH
DAPI

mediummes

Figure 2 | Fgf4 treatment induces sometrophoblast-lineage character in STAP cells.
a, Schematic of Fgf4 treatment to induce
Fgf4-induced stem cells from STAP cells. b, Fgf4
treatment promoted the generation of flat cellclusters that expressed Oct4-GFP at moderatelevels (right). Top and middle: days 1 and 7 ofculture with Fgf4, respectively. Bottom: culture
after the first passage. Scale bar, 50 mm.
c, d, Immunostaining of Fgf4-induced cells withthe trophoblast stem cell markers integrin a7(c) and eomesodermin (d). Scale bar, 50 mm.
e, qPCR analysis of marker expression.
f, g, Placental contribution of Fgf4-induced stem
cells (FI-SCs) (genetically labelled with constitutive
GFP expression). Scale bars: 5.0 mm (f (left panel)and g); 50 mm (f, right panel). In addition toplacental contribution, Fgf4-induced stem cells
contributed to the embryonic portion at a

cag-GFP

1.2

***

1.0
0.8
0.6
0.4
0.2
0

NS

JAKi– JAKi+ JAKi– JAKi+
ES
Oct4

FI-SC
Nanog

Rex1

k

Relative expression levels(TS=1.0)

1(FI-SC, ES, STAP-SC)

Relative expression levels(ES=1.0)

0.10
0.00

0.05

Height0.15

derivation of STAP stem cells from STAP cells involves the loss ofcompetence to differentiate into placental lineages. Consistent withthis idea, STAP stem cells show little expression of trophoblast marker
genes (Fig. 1f).
We next examined whether an alteration in culture conditions could
induce in vitro conversion of STAP cells into cells similar to trophocag-GFP
cag-GFP
blast stem cells8,9, which can be derived from blastocysts during prolonged adhesion culture in the presence of Fgf4. When we cultured
auconditions (Fig. 2a; one cluster per well
STAP cell clusters under similarijin a 96-well plate), flat cell colonies grew out by days 7–10 (Fig. 2b, left;
typically in ,30% of wells).
CD45The Fgf4-induced cells strongly expressed
the trophoblast marker proteins9–12 integrin a7 (Itga7) and eomesodermin (Eomes) (Fig. 2c, d) and marker genes (for example, Cdx2; Fig. 2e).
These Fgf4-induced cells with trophoblast
marker expression could
100
be expanded efficiently in the presence of Fgf4 by passaging for morethan 30 passages with trypsin digestion every third day. Hereafter,
STAP
these proliferative cells induced
from STAP
100cells by Fgf4 treatment
are referred to as Fgf4-induced stem cells. This type100
of derivation into
trophoblast-stem-like cells is not common with ES cells (unless genetTS
ically manipulated)13 or STAP stem cells.
100
In the blastocyst injection assay, unlike STAPFI-SC
stem cells, the placental contribution of Fgf4-induced stem cells (cag-GFP-labelled) was
ES STAP-SC
observed with 53% of embryos (Fig. 2f, g; n 5 60). In the chimaeric

FI
-S
C
FI
-S 1
C
FI
2
S
G
C
FP 3
G ES
FP 1
G -ES
FP 2
-E
S
3

各領域がクラスタに相当

Bright-field

2

3

4

5

6

(TS, ES, STAP-SC)
NS

1.2
1.0
0.8
0.6
0.4
0.2
0

TS

JAKi+

P values(culture
without
pluripot
FI-SCs w
treatmenshown a(j) or tro
NS, notgroups w
n 5 3. S
three plu
trophobsignifica

JAKi–

FI-SC
Cdx2

Eomes

Elf5

Itga7

30 JANUARY©2014 Macmillan Publishers Limited. All rights reserved(ES, STAP-SC)がクラスタになる全てのデンドログラムに相当する領域

Oct4
Cdx2
Eomes

17

系統樹の信頼度は？
Dn

計算結果

human
cow
rabbitmouse
opposum

計算
サンプルサイズ

n

D¥

計算

？

18

「計算」はブラックボックスとみなす
計算

Dn
human
cow
rabbitmouse
opposumℓ1 (θˆ1 )

計算結果ℓ 2 (θˆ2 )maxℓ15 (θˆ15 )19

反復計算してカウントだけ使う
Dn

計算

Yes

Dn*1

計算

Yes

Dn*2

計算

No

*10000
n

D

計算

Yes

#{Yes}BP =≈ P(Yes | Data)
10000
Bootstrap Probability
ブートストラップ確率20

ランダムに要素をコピーする
Dn

Dm*1
*2
m

D

1

2

3

n

1

2

3

m

Dm*10000

ブートストラップ・リサンプリング

Efron (1979)
Felsenstein (1985)
21

Bootstrap: n=10,242 columns
human
cow
rabbitmouse
opposum

Assume i.i.d. for columns

Dn

*1
m

*100
m

D ,..., D

We assume m = O(n) , typically m=10242, 5000, 20000, say
In the “m-out-of-n bootstrap” 0 < m  n , typically m=30, say
22

BPはベイズ事後確率
Efron and Tibshirani (1998)
*1
n

*2
n

#{Tree}BP = P(Tree | data) ≈
100

*3
n

D , D , D ,...

1

1

Dn

D¥

2

3

4

3

2

4

5

5

1

2

3

4

5

1

2

3

4

5

∗
のバラツキ（分散）はmに反比例する
Shimodiara (2002)のマルチスケール・ブートストラップ：𝐷Q

mを変化させたBPから，各仮説領域までの符号付き距離や，仮説境界の平均曲率がわかる！
23

たとえば，多重比較
𝑿~𝑁(𝜼, UOΣ)
𝐻O: 𝜂O ≥ 𝜂M,…, 𝜂O ≥ 𝜂H

3.5
𝑿 = 4.0
0.5
3.5
𝜼 = 3.0
0.5𝑿𝜼𝜂O ≥ 𝜂M

𝜂O ≥ 𝜂L
𝐻O 凸多面錐24

多重比較の場合
𝑿~𝑁(𝜼, UOΣ)𝒀 = 𝑛Σ KO/M𝑿
𝝁 = 𝑛Σ KO/M𝜼

𝑑O(𝑌) = max 𝑐M 𝑋M − 𝑋O , … , 𝑐H 𝑋H − 𝑋O

𝒀~𝑁 𝝁, 𝐼

confidenceQb 𝜇 ∈ 𝐻O 𝑦 = 𝑃(𝑑O 𝑌 > 𝑑O (𝑦)|0)𝐻M𝑦Least favorable configuration

𝜂O = 𝜂M = ⋯ = 𝜂H

𝐻L𝑑O (𝑦)

仮説領域（凸錘）の頂点

𝐻O25多重比較は保守的
各領域のType-I error
𝑃 confidenceQb 𝜇 ∈ 𝐻j 𝑦 < 0.05 ≤ 0.05 under 𝐻j
仮説境界上で等号なら嬉しいが，等号になるのはl.f.c.のみ

仮説の信頼集合𝑆 𝑦 = 𝐻j : confidenceQb 𝜇 ∈ 𝐻j 𝑦 ≥ 0.05
𝑃 𝐻j ∉ 𝑆 𝑦≤ 0.05 under 𝐻j26

Region and a data-pointv

AU = PV + O(n−3/2 )
!
1
u2
v=−
+
64
3
!
"
1
u2 #
∂H = (u, v) : v = −
+
64
3
!
"
1
u2 #
H = (u, v) : v ≤ −+64
3
σ2 = τ 2 = 1u27

regionH = (u, v) | v ≤ −h(u), u ∈ Rq

nd Tibshirani (1998)
considered the local
coordinates (u, v) ∈ R!
"q
H = (u, v) is| vdenoted≤ −h(u),u∈R
as R(h). The boundary
surface
−2 −2

(8)

=+
PV
+ O(n
)
DAUDAU= PV
O(n) ModelMultivariate
Normal. , u ) ∈ R , v ∈ R for representing H in a neighbourhoodof (0,!
" 0qq

s R(h). The boundary surface

+1 ,

q+1
y=
(u,
v)
y! = (u,v)∈
R
is denoted as B(h).

∂H = (u, v) | v = −h(u), u ∈ Rq

= dim
dimuu
"qq =

q For a continuous functionn 1 (Region
boundary
surface).
∂H and= (u,v) | v = −h(u),u ∈ Ridentity
matrix Iq+1 ,es region
B(h).

Y ∼ Nq+1 (µ, Iq+1 ).

H"= R(h)!
Higher
order
terms
of
doublebootstrapgher
order terms
of ofdoubleprobabn arbitrary-shaped
regionpositive bootstrap
volume.
Givinganpro
obsqu ∈ R (1)
∼ Nq+1 (µ, Iq+1H).= (u, v) | v ≤ −h(u),Let H ⊂ Rq+1 be an arbitra
erested in testing
theH µ̂(H|y)
nullhypothesisµ∈
H
againstthealternboundary
surface:=(0,−h(0))∈∂H
region:∂H=
B(h)= R(h)µ̂(H|y)= (0, −h(0)) ∈ ∂H
as
R(h).
The
boundary
surfaceed region of positive volume. Giving
an observed
value
y, we are interested in
By ∂Htaking ofthe H
origina point on ∂H
the axes
The boundary surface
isatassumed
toandberotatingsmoothin!+)θ(6h
"mj
(3h
−
6λ
h
h
+
(6h
−
2γ
h
h
−
4h
hljmi
h
mmi
0
ml
mli
i
mmij
1
mi
ml
−
6λ
h
h
)θ
−
2γ
h
h
−
4h
hvmiisTheh)θNullhypothesis:
v.s.
Alternative
hypothesis:
hypothesis
µ
∈
̸
H.
bou
mmi
0 ml
mli iµmmij
1of mimj q spacemli θl
the null hypothesis∈
Hagainstthealternativeareorthgonal
coordinates
thetangent
andthecoordi∂H = (u, v) | v = −h(u), u ∈ R∂H = B(h)normal to the tangent space. The smooth function h(u) is then28
paper. in thisurface∂Hof
H
isassumedtobesmoothm 4 (Asymptotic expansion of double bootstrap probabilit

Chi-square test (very conservative)
distance2 ~ chi-square distribution with df=2過大
p=0.1851.836

過大
p=0.206
1.778

29

Multiple Comparison (conservative)
X ~ N 3 (η, I 3 / n)⇔

Y ~ N 2 ( µ, I 2 )η1 ≥ η2 , η1 ≥ η3⇔v ≤ − |u |/ 3

u = n / 2(η3 − η2 )v = n / 6(η2 + η3 − 2η1 )妥当
p=0.073過大
p=0.06830

Singed LR test (rejecting too much)
distance ~ N(0,1)

過小
p=0.033
1.836

妥当
p=0.038
1.778

31

for any σ 2 > 0. The bootstrap probability of H for a given y is defined as

ocedures are utilized for computing approximate p-values. W
proceduresare (y,utilized
forσ2computing
approximatep-v
BP(H|y).
We
ignorethesampling
Y for∼ N estimating
σ I )
(2)
BP (H|y) = P (Y ∈ H|y),
∗

2

σ2

∗

replicatesBPσ2 (H|y)generate= P (Y ∗B∈ bootstrap
H|y),

(3)

q+1

q+1

∗ (=Bayesian PP)
Bootstrap
probability
n0. Thegeneratebootstrap
replicateYofbybootstrap probability
for indicatesa given y theis definedas Y∗ respect to (2). In practice, we mwhereof PH(·|y)
probability
with(3) bootstrap
in the theoreticalargument
For a
can generate
replicate
Y of throughout.
Y by
∗

2

2
Y (1),∼ Ntheσ Ishould
Y in
scale
σ
= 1 in (2). The
q+1 (y,q+1Y) , . . be∗
2 .,Y ,∗1∗B

Y ∼ Nq+1 (y, σ Iq+1 )) indicates the probability with respect to (2). In practice, we maybootstrap replicatesand compute the frequency
BP
1 (H|y) is often used as an approximate p-valu
The bootstrap probabilityof H for #{Ya giveny is defined as∈ H, b = 1, . . . , B}
! σ2 (H|y) =BP∗b

Y , . . . , Yprobability,
0. The bootstrapof H for aB given y is define∗1

e the frequency∗Bsense explained in Section 2.5, and there are severa
estimating BP (H|y). ∗We ignore the sampling error of O(B) and use oBPσfor2 (H|y)= P (Y ∈ H|y),
σ2

−1/2

∗
(3)∗bin∈theargument throughout.
For approximating the distribution#{YH,btheoretical= 1, . . . , B}!
2
BPσ2 (H|y) =σ B
Y in (1), the scale should be σ 2 = 1 in (2). The ordinary bootstrap probabil−1/2
g BPσ2 (H|y). We ignore
thesampling
error
of
O(B
) and use probability
only
BP isBPinterpretedastheBayesianposterior
H
(H|y)
is
often
used
as
an
approximate
p-value. of
However,
it is biased in1

BP (H|y) = P (Y ∈ H|y),dicates the probability with respect to (2). In practice, w

if the probability
prior distribution of mu
is uniform.
indicates thewithrespect to (2). In prac

eoretical argument throughout. For approximating the distribution ofsense explained in Section 2.5, and there are several attempts to adjust the bias

strap replicates2

Efron and Tibshirani (1998)e scale should be σ = 1 in (2). The ordinary bootstrap probability

ootstrap replicates532

BP is even worse
過小
BP=0.019
p=0.033 (signed LR test) 妥当
BP=0.037
p=0.038 (signed LR test)33

2

Double bootstrap

alibrating BPσ2 (H|y), we firstconsider the projection of y onto the bou
BP 2 (H|y + ) = P (Y ∗∗ ∈ H|y + ).
σ

Double bootstrap probability

σ 2 (H|y), we first consider the projection of y onto the boundary
ecalibrating∂H. The BP
projectionis the point on ∂H that is closest to+ y;
We calibrate BPσ2 (H|y) by the distribution of BPσ2 (H|Y ). The doubl

face ∂H. The projection is the point on ∂H that is closest to y;
of y onto
thedefinedboundary=arg
min ∥yassurface:− µ∥,probability of HProjection
for µ̂(H|y)
a giveny isµ∈∂Hµ̂(H|y) = arg"min ∥y − µ∥,

#
µ∈∂H
!
q+1
2 = P BP 2 (H|Y + ) ≤ BP 2 (H|y) | µ̂(H|y)
+
2−
2µ
(H|y)
∥y − µ∥2 = !DBP
(y
)
.
We
generatebootstrap
replicate
Y
by ,
σ
,τ
σ
σ
i
i
i=1

ere ∥y − µ∥2 =q+1Adjusting BP
2 using resampling from the projection
+
(y
−
µ
)
.
We
generatebootstrap
replicate
Y
by
i
i
i=1

+
2 When σ 2 = τ 2 = 1, this is t
where the probability
iscalculated
by
(4).
∼
N
(µ̂(H|y),τ
Iq+1bootstrap)
H.YSHIMODAIRA/multiscale-double
q+1
+
2

Y

∼ Nq+1 (µ̂(H|y), τ Iq+1 )(4)

probability
ofbootstrap
H for a given
y is defined
as 1,1 (H|y).
doubleprobability
DBP
2
+
me
τ
>
value
boot
2 0. For each observed !
+y , we generate a second-level
some τ > 0. For each observed value y , we+generate a second-level "bootstrap
DBPτ 2 ,σ2 (H|y) = Pτ 2 BPσ2 (H|Y ) ≤ BPσ2 (H|y) | µ̂(H|y) .

ate
licateThevariances are usually σ 2 = τ 2 = 1 and we simply denote DBP or DBP(H|y2.3 Multiscale bootstrap
DBP1,1 (H|y). Efron and Tibshirani(1998) +called
DBP as a calibrated confidence leve∗∗2
2 )YY ∗∗∼∼accurate.
NN(y(y, +
σ , Iσq+1
q+1
Iq+1
)
q+1
mentioned that DBP is third-orderHall (1992), Efron and Tibshirani(1998)
2
We
will show
later in Section
5.3 that theBP
doublebootstrap
probabilityforσ> 0, τ
2Another
approachfor calibrating
(H|y)
is
the
multiscalebootstrap
σ34
scomputing
expressed asymptotically as
mputingσ =τ =1+
∼ N2 (µ̂(H|y),Y + ∼ N2Y(µ̂(H|y),I2 ) I2 )contour
surfaceof
BP=0.019$%$%+
DBP(H|y) = P BP(H|Y) ≤ BP(H|y) | µ̂(H|y)+|y) = P BP(H|Y ) ≤ BP(H|y) | µ̂(H|y)
Computing DBP:

Y + ∼ N2 (µ̂, I2 )Y + ∼ N2 (µ̂, I2 )$%DBP
$ = P BP(H|Y ) ≤ 0.019% | µ̂
+DBP = P BP(H|Y + ) ≤ 0.019 | µ̂

DBP=0.049
approach for calibrating BPσ2 (H|y) is the multiscalebootstrap. The

2 (H|y) is the multiscale bootstrap. Ther calibrating
BPσofootstrap
probability
Shimodaira (2008) is defined byobability of Shimodaira& (2008)
is defined' by−1
NBPσ2 (H|y) = Φ σΦ (BPσ2 (H|y))&−1

BP(H|y)=0.019

(6)

'

NBP
(H|y)
= Φ σΦ
(BPσ2 (H|y))Φ(·) isσ2thecumulative
distributionfunction of the standard normal(6)dis-0, 1), and we denote the density function as φ(·). For several σ 2 values,umulative distribution function of the standard normal dis-

σ 2 , we compute BP 2 (H|y), i = 1, . . . , S, and extrapolate NBP 2 (H|y)

35

− Φ β0 − β1 − β2 + β3

(28)

DBP adjusts the bias of BP!
"
PV(H|y) ≃ 1 − Φ β0 − β1 − β2 + β3P(BP<0.05)

"

妥当

2
2 (H|y) ≃ (β0 − β1 − β2 ) − β3 σ
DBP1,σP(DBP<0.05)µ ∈ ∂Hθ = 0.0, 0.5, 1.0, . . . , 3.0µ = (θ, −h(θ)) ∈ ∂H
DBP = PV + O(n−3/2 )DBP=0.05

DBP=0.049
BP=0.019

BP=0.05

妥当
DBP=0.039
BP=0.037

DAU = PV + O(n−2 )20

er terms of double bootstrap probabilitiesµ̂(H|y) = (0, −h(0)) ∈ ∂H

36

1

simpler.

Introduction
DAU = Φ̄(g(−1))σ viaApproximately unbiased
p-valuesnbootstrapσ =Another mapproach f
1 Introduction
AU = Φ(−BZ(−1))
Multiscalebootstrap
X = {x , . . . , x }AU = Φ(−BZ(−1))

σ2

2

2

1

n

Another approach for calibrating BP
is the mu
malizedbootstrapp
σ 2 (H|y)m outof n bootstrap : Politis and Romano(1994), Bickel et al. (1997)
bootstrapY

DAU = Φ(−DBZ(−1))

DAU= Φ(−DBZ(−1))
malized bootstrap probabilityof Shimodaira(2008) isX ∗ = {x∗1 , . . . , x∗m }X = {x1 , . . . , xn }AU = Φ̄(f (−1))

Y

Y

∗

$

−1
NBPσ2 (H|y)= Φ̄(f
Φ2 σΦ
(BPσ2 (HAU =(−1))
for σ > 0. Φ(·) is t

The idea∗of multiscalebootstrap
Shimodaira(2002,
2004,
2008)
DAU = :Φ̄(g(−1))+++
X = {x∗1 , . . . , x∗m2}
X = {x1 , . . . , xn }for σ > 0. Φ(·) is the cumulative
distributionfunctiontribution
N (0,1), anDAU= Φ̄(g(−1))

nY += 1), and we denote the density
tribution σN (0,say,2 σ12 , n. . .function, σS2 , we ca
m
σ =
+
∗∗
X + = {x+X ∗∗ = {x∗∗
,
.
.
.
,
x
}
1 , . . . , xn }
m
1
m
2
2
2
2 (H|y),say, σBP. 2. (H|y), σS , weBPi=1, .More. . , S,nandextrapolate
BPσto
We computeBP forAnother approachfor calibrating
iscomputethe multiscaleThe
tobootstrap.
σ =
−1.
s
1, . σi+YY ∗∗(equivalently m = −n )2to σ = −1. More(2008)
specifically,
we fitbya model
malized bootstrap probability
isdefined∗∗∗∗ of Shimodaira
∗∗
2
Y∗

2

X = {x , . . . , xm }n/m
The BP with m 1= -n is denoted
as AUσ ( ==Approximately
Unbiased)Y∗∗

Another approachfor calibrating %BPσ2 (H|y) is the mul2$ !
−1 σ 2 NBPσ 2 (H|y) = 1 − Φ(β0 + β1 σ +BP
NBPσ2 (H|y) = Φ σΦ (BPσ2 (H|y))37

e- x 2 /2

0.0

0.1

f ( x) =

1
2p

0.2

0.3

0.4

正規分布の上側確率

-2

0

2

4

-4

-2

0

2

4

F ( z ) = ò f ( x)dx
F (1.64) = 0.05

0.0

0.2

z

0.4

¥

0.6

0.8

1.0

-4

1.64

38

2.3
Multiscalebootstrap
Multiscalebootstrap
2.3
Multiscalebootstrap"#
m=-n への外挿2! ∞2σ 2 ) = Φ̄−1σDBP
(H|y)
σ ,1
x2
1
2

! ∞
! ∞ 1 − x22− 2Φ(−z) =√ 1 e − 2x2 dΦ̄(z)e=√Φ̄(z) = Φ(−z) =dxΦ̄(z)= Φ(−z) = z √2π e
d
2π
z
2π
z"
#
!
∞
2
2
"
2
−1
1 −)x = −σΦ BP 2 (H|y) #σ Φ̄(z) = 1 − Φ(z) =AU = Φ(−BZ(−1))
BZ(σ
√BZ(σ
e 2 )2 =dx−σΦ−1 BPσ 2 (H|y)
σ
2π
z
"
"
#
2
−1 "
multiscalebootstrap
2
2
2
−1
DBZ(σ
)
=
−Φ
DBPσ2 ,1 (H|y)
σ
=
−1
σ
=
1
2 (H|y)
2
−1
BZ(σ
)
=
−σΦ
BP
σ
AU==Φ(−DBZ(−1))Φ(−BZ(−1))
DAU
DBZ(σ ) = −Φ DBPσ2 ,1 (H|y)
ブートストラップ確率
2.3 Multiscalebootstrap #"#BP=0.019"2−1 "
2 (H|y) #
f
(σ
)
=
σ
Φ̄
BP
2
−1
σ
DBZ(σ ) = −Φ DBPσ2 ,1f(H|y)
2
−1!
(σ
)
=
σ
Φ̄
BPσ2 (H|y)"2−1
DAU
AU==Φ(−DBZ(−1))Φ̄(f (−1))BZ(σ)bootstrap= −σΦ" BPσ2 (H|y) #double"2#−1 "マルチスケールブートストラップ2
g(σ
)
=
Φ̄
DBP
−1
σ 2 ,1 (H|y) #!
"f (σ ) = σ Φ̄ BPσ2 (H|y)
2
−1
2
g(σ
)
=
Φ̄
DBP
(H|y)
2
−1
AU=0.043
DBZ(σ ) = −Φ DBPσσ2,1,1 (H|y)
DBP=0.049
AU = Φ̄(f(−1))
DAU
Φ̄(g(−1))
"
#
2
σ
−1
2
4
DAU=0.055σ Φ̄
BPσ2 (H|y) ≃ βダブルブートストラップ
0 + β1 σ + β2 σ2σ 2σ
マルチスケール・ダブルブートストラップ"# = Φ(−BZ(−1))
AU2ing DAUBPσ (H|y)is thebootstrap.
The nor2 multiscale−1= Φ̄(g(−1))
2 ,1 (H|y)
g(σ
)
=
Φ̄
DBP
σ
2 (H|y) is the multisc
Another approach for calibrating AUBPσ=Φ(−BZ(−1))
We apply the multiscale bootstrap to DBP for getting DAU

Extrapolation to sigma =-1

Equivalently, we could say applying double bootstrap to AU for getting DAU

39

DAU = Φ(−DBZ(−1))

contour curves of p=0.05
BP, AU, DBP, DAU

40

Theorem 3 (Fourth-order accurate p-value). For the H and y = (0, λ0 −ＤＡＵの棄却確率は誤差が小さい
Theorem 1, an approximately unbiased p-value of fourth-order accuracy

Rejection probabilities P(p<0.05)
asymptotically as

Error: DAU < {DBP, AU}! < BP

PV(H|y) ≃ 1 − Φ λ0 − γ1 + λ0 γ2 + 3γ4 − γ1 γ2 − 43 (2 + λ20 )γ3!
"= 1 − Φ β0 − β1 − β2 + β3
BP

AU

!
"
PV(H|y) ≃ 1 − Φ β0 − β1 − β2 + β3

−1

!

"

Φ̄ DBP1,σ2 (H|y) ≃ (β0 − β1 − β2 ) − β3 σ 2
DAU

DBP

θ = 0.0, 0.5, 1.0, . . . , 3.0µ ∈ ∂Hµ = (θ, −h(θ)) ∈ ∂H

41

棄却確率 (0.05に近いほど良い）
H. SHIMODAIRA/multiscale-double bootstrap
Table 2
過大
Rejection probabilities(in percent) at significance level α = 5%.
u
BP
AU2
AU3
DBP
DAU
MCB

0.0
13.39
7.655
6.609
6.619
6.476
5.000

0.5
8.894
5.171
4.718
4.590
4.660
3.340

妥当

1.0
6.678
4.459
4.493
4.202
4.481
2.880

1.5
5.676
4.447
4.746
4.364
4.746
2.783

2.0
5.253
4.628
4.982
4.610
4.981
2.768

2.5
5.086
4.801
5.080
4.795
5.084
2.766

妥当

3.0
5.027
4.912
5.081
4.905
5.092
2.766

過小

ometry of smooth surfaces

section, we discuss only geometry of smooth surfaces via simple but tedious c42
without any probability argument. The results will be used in Section 5 forderNew Discovery?

58%

2%

43

ベイズの事後確率
1=human 3=rabbit 4=mouse

1=human 3=rabbit 4=mousep ³ 0.05

ブートストラップ法（m=n)で計算

44

Shimodaira-Hasegawa testは安全すぎる

SH ³ 0.05

対数尤度に多重比較法を適用してモデル選択の多重性を調整
45

頻度論のp-値 (近似的に不偏 AU)
1=human 3=rabbit 4=mouse

マルチスケール・ブートストラップ法（m=-n)で計算
46

曲率を反転するとベイズは頻度論
頻度論のp-値

ベイズの事後確率- b1

Jerzy Neyman(1894-1981)

qb

1

b0

H

qˆ

b0

b1

qˆ
Thomas Bayes(1702-1761)

頻度論

ベイズ統計学AU = Φ ( β0 − β1 ) + O(n−3/2 )BP = Φ ( β0 + β1 ) + O(n−3/2 )=47

Rescaling the whole picture
Shimodaira (2002)

s 2 =1s 2 =1´

1sBP = Φ ( β0 + β1 ) + O(n−3/2 )é b0ùBP(s ) = F ê + b1s ú + O(n -3/2 )
ës
û
2

{

(

NBP(σ 2 ) = Φ σ Φ −1 BP(σ 2 ))} = Φ ( β + β σ ) + O(n )
2

0

1

−3/2

48

ベイズと頻度論のギャップを埋める{(

NBP(σ 2 ) ! Φ σ Φ −1 BP(σ 2 ))}

正規化したブートストラップ確率(NBP)
スケーリング則NBP(σ 2 )AU = NBP(−1)
Jerzy Neyman(1894-1981)

Thomas Bayes(1702-1761)

BP = NBP(+1)

ベイズ統計学頻度論
s 2 = -10s 2 =1m = −n

m=∞

m= n

σ2 =

n
m

AUの計算手順： いくつかのm>0でBPを計算して，NBPに変換してからm=-nへ外挿49

We can’t compute nBP(-1) for cones
smooth surface

é b0ù
BP(s ) » F ê + b1s úësû
2nBP(s 2 ) » F ( b 0 + b1s 2 )

cone

éb
ù
BP(s 2 ) » F ê 0 + b1 úës
ûnBP(s 2 ) » F ( b 0 + b1s )50

Taylor expansion of nBP using k terms
Shimodaira (2008)

é k -1 (s 2 - s 2 ) j ¶ jy (s 2 ) ù
0
nBPk (s ) = F ê å
ú
2 j
j
!
¶
(
s
)
êë j =0s 02 ú
û2y (s 2 ) = b 0 + b1s

AU k = nBPk (-1)

Taylor expansion（k=1,2,3,4）51

Scaling-law and geometry
Shimodaira (2008)()

NBP(σ 2 ) = Φ β0 + β1σ 2 + β 2σ 4 + β3σ 6 + + O(h2 )distancemean curvature

p-value = Φ ( β0 − β1 + β 2 − β3 +) + O(h2 )j!
∂ 2 j h(u)β j (u) = j∑
2 j! j1++ j p−1= j j1 ! j p−1 ! ∂u12 j1 ∂u 2p−1j p−11v
u=(u1,..., up-1)

v = - h(u )Using a new asymptotic theory with “nearly flat surface h” instead of large n.
52

Extrapolating NBP to m=-n
NBP(σ ) ! Φ {σ Φ ( BP(σ ) )}2

Plotting :

−1

2

(

σ Φ −1 BP(σ 2 )

)

2
2
4
2( k−1)Fitting a model : ψ (σ ) = β0 + β1σ + β 2σ ++ β k−1σ

BP=0.0078
k=2: AU=0.057
k=3: AU=0.052
p=0.05

-1

1

s2
53

What NBP is calculating ?
Shimodaira (2008)

s 2 = -1s 2 =1()

NBP(σ 2 ) = Φ β0 + β1σ 2 + O(n−3/2 )curvature

NBP(−1) = p-value + O(n−3/2 )54

What is an unbiased test?

{(“similar” on the boundary)}P p (qˆqˆ) < a | q q = a ,q Î ¶q H , 0 < a < 1

"

qˆ
q

qˆ

qˆ*

55

Matching priorsθ ↔ (ϕ , λ )ϕ ∈, λ ∈p−1

ϕ

H = {θ : ϕ ≤ ϕ 0 }ϕ = ϕ0Tibshirani (1989) assumes orthogonality Iϕλ (ϕ , λ ) = 0

π (θ ) = g( λ )I (ϕ , λ ) | J (ϕ , λ ) |1/2
ϕϕd log π (θ )= −2β1dϕ−1(curvature of the boundary surface of H)β1 ∝1/ radius

Reweighting the bootstrap by this matching prior makes BP=p-value

56

Correcting a deviation from normality
Efron (1987), Efron and Tibshirani (1998), Shimodaira (2004)

1 E ((qˆ - q )3 )a=
6 E ((qˆ - q ) 2 )3/2

acceleration constantnBP(s 2 ) = F éë( b 0 - b1s 2 ) - a(2 b 02 + s 2 ) ùû + O( n -1 )p -value = F éë b 0 - b1 + a (1 - b 02 ) ùû + O(n -1 )F -1 ( p-value) - F -1 (AU(-1)) = a b 02 + O(n -1 )57

Two-step multiscale bootstrap
Shimodaira (2004)

Dn

Dm*11

Dm*12

Yes

Dm*21

Dm*22

No

Dm*10000
1

Dm*10000
2

Yes#{Yes}
10000= BP(s 2 ) + O p (n -1/2 )BP(s 12 , s 22 ) =with s2= s 12 + s 22

s F -1 (BP(s 12 , s 22 )) = s F -1 (BP(s 2 )) + n -1/2 as 12s 22s -4 ( b 0 - s 2 ) + O p (n -1 )c.f. Three-step multiscale bootstrap estimates six geometrical parameters andgives third-order accuracy
58

Many hypotheses: FDR, visualization
qˆ
H1

H2

H3
H4

p1 , p2 , p3 , p4 ,...

59

Taylor seriesas
Bymay
taking
at aandpoint∂H andaxes properly,, uq tangent
notthebeoriginon ∂H,
theonaxes
uq1 , . rotating. . , uq arethe
slightlyrotated ufrom
1 , . . . the
q
q
q
q
#
#
#
#
#
normaltothetangent
space.
The
smoothfunction
h(u)
is
then
represented
innormal to the tangent space. The smooth function h(u) is then represented in the
the
q
q
q
q
q hij ui uj +
h(u)
=
h
u
u
u
+
·
·· ,
i j k
# # of the tangent### and v is the coordinate ijk
are orthgonal
coordinates
spaceof thedirection
space.
i=1 j=1 hijk ui uj uki=1
h(u) series= as+ ·j=1· · , k=1(9)TaylorTaylorseriesas hij ui uj +Asymptotic theory of
where the th
are denoted as
Definition 2 (ClassS).
Let h(u) be accuracy
a smooth function of u ∈ R . The function h
4 derivatives
order
$
$
i=1 j=1i=1 j=1 k=1(9

qq Theqq smooth functionqqqq
qqh(u) is then represented in thenormal to the tangent space.
#
#
#
#
#
#
#
#
#
#
q
h(u)
hhijijuuiiuujj +
hhijk
(9)
h(u) == as+
uiiuujjuukk ++······ ,,(9)ijku
where the derivatives are denoted
i=1
j=1
i=1
j=1
k=1
i=1 j=1i=1 j=1 k=1
Taylor series as1 ∂ 2 h(u) $1 ∂ 3 h(u) $hijas= 1 asymptotically=iswheresaid to be
class
S are
if $$itdenoted
is expressed$ , $$ hijk as$.
1 ∂ 2 h(u)
∂ 3∂u
h(u)
derivatives
where thethederivatives
aredenoted
as
2
∂u
6
∂u
∂u
∂u
0
0
i
j
i
j
k
q
q $ ,
q
q
q
hij = #
h=.
$ijk # # ##2 ∂u ∂u 0
i ∂uj ∂uk 0$$$ 6 ∂uto
= i wej allow
h1ij1∂u∂2i2h(u)+hijk∂Su33h(u)uj uk +
· · , ∂ 3so(9)origin$$ ·below
More h(u)generally,hbelongs
class
that$$ thejidefined
2u
1
h(u)
h(u)
1
$
$
$
1
∂
h(u)
1
∂
h(u)
−1/2
−1
$
$
$ h(u)$u huijk+
3 S
hhj=1
=
,,=
h
=
..hijkl ui u
hij≃
=
O(n
)
=
O(n
)
$$to
$$=
$
$
i=1
i=1
j=1
More generally,wehbelongs
class
defined
below
so
that
the
origin
k=1
=
h
=
ij
ijk
h
+
h
u
+
h
u
u
+
h
u
u
u
,
(10)
1 ∂ 2 h(u)
1
∂
h(u)
ij
ijk
0
i
i
ij
i
j
ijk
i
j
k
j
k
l
$ allow
$
−1
∂u
6
∂u
00 6
2)22∂u∂u∂u∂u
0∂ui i∂u
∂u
6
∂u
∂u
∂u
jj 00=
ii∂u
jj∂u
kO(n
j
i
j
k
i∂u
k
hij may= not be$ on=∂H,
O(n−1/2h=)$ijk u , . . . , u are slightly rotated from the tangent
and the axes
1 6 ∂ui ∂uq j ∂uk 0
2 ∂ui ∂uj 0$$$$
where the derivatives11aredenotedas∂∂22h(u)11 ∂∂33h(u)$$$belowfrom
h(u)
h(u)
−1/2
−1
$
−1/2
−1
More
we
belongs
definedsothemaynotgenerally,
be on ∂H,
thehaxes
u1the, . . to.)), uclassslightlyrotatedthe astangentAsn=becomescoeﬃcients
approachzero
hij origin=O(n−1/2
hhijijand=O(n)) that$larger,$$ =q hhare=allow= O(nO(n=S=asymptoticlly
O(n
ijk
$
ijk =
−1
−1/2
−1
−1/2
22 ∂u
66∂u
00 approachthethecoeﬃcientsh0 = O(1),), jjh∂u=0ij0 O(nspace.
∂uii∂u∂ujjare∂uii∂u∂u∂u
As n becomeswherelarger,coeﬃcients
zerohasymptoticllyas= O(n ), h)ijk = O(n ),i = O(nijkkh$
3
1 ∂ 2−1h(u)1 ∂derivatives
h(u) $$ are O(n−(k−1)/2−1/2
$
and
h
=
O(n
).
The
k-th
order
)−1/2
for)) ktangent≥−2
1, becausAsnn∂H,
becomeslarger,thecoeﬃcientszero
asymptoticlly
asijk
space.
As−1
becomeslarger,thecoeﬃcientszero
asymptoticlly
as hh
= O(n
O(n
ij
h).ijthe==. rotated
−3/2
−3/2
may
not
be
on
and
axes
u$1 , . . happroach.approach, uequality
areslightlyfromthe−(k−1)/2
$ up
ij =
ijk
qO(n
h
=
O(n
Here
≃
denotesthecorrecttoO(n)
erring
O(n
).
and
h
=
O(n
).
The
k-th
order
derivatives
are
)
for
k
≥
1,
becauseijkl
ijk∂uiLet∂uj h(u)6 ∂ui ∂u0曲面の微分係数j ∂uk 0 of u ∈ Rq . The function h
Definition 2 (Class 2S).
be a smoothfunction√−1−(k−1)/2−1). The k-th order derivatives are!
−(k−1)/2and=O(n)) for≥becauseのスケーリング則!√and hhijk= O(nO(n). The k-th
O(n
for
≥ 1,
1,factor
because n.
coordinates
u1 , .order. . , uqderivatives
as well asareh(u)scaledbykk theijk theq areq
the coordinates
well as h(u)suchare asscaledn. ij ui uj are used,
Theusummation
hij uby= factor
where the1 , . . . , uq as convention
space.
i uj thei=1as j=1 has
q = O(n−1/2 )
is
said
to
be
class
S
if
it
is
expressed
asymptoticallyAsnbecomeslarger,thecoeﬃcients
approachzero
asymptoticllyhij√
Definition 2 (Class S). Let h(u) be a smooth function of u ∈ √R . The function hthethe coordinates
coordinates uu11,,......,,uuqq as
as wellwell asas h(u)h(u) areare scaled
scaled byby thethe factor
factor n.
n.

free indeces i, j, k, l run through 1, . . . , q if appeared twice in a formula.
and hijk = O(n−1 ). The k-th order derivatives are O(n−(k−1)/2 ) for k ≥ 1, becauseqis said to be2 class
S S).
if it
is ≃expressed
asymptotically
h(u)
h
ui +
ui uashijkluui∈
uj uR
(10)
0 + hibe
ij ui uj + hijk
j uk + of
k ul., The function
Definition(ClassLeth(u)ahsmoothfunctionh√the coordinates
u1 ,care. . . , uofq as
wellup
as to
h(u)
are
scaledby theO(nfactor−3/2−2
We taketerms
O(n) ignoring) n.
is said to be whereclass theS coeﬃcients
if it is expressedas
10 asymptotically 10
areh(u) ≃ h0 + hi ui + hij ui uj + 10
h10ijk ui uj uk + hijkl ui uj uk ul ,(10)h0 = O(1), hi = O(n−1 ), hij = O(n−1/2 ), hijk = O(n−1 ), hijkl = O(n−3/2 )3.2 Asymptotic expansion of the bootstrap probability
h(u) ≃ harehi u=
hij uihuj =+ O(nhijk−1
ui),
uj u
+
ukhu
, = O(n−1 ),
(10)
−1/2
0 +h
i +
ijkl ui uj ),−2l
where the coeﬃcientsO(1),hk −3/2=hO(n. Here ≃ denotes
0 the equality icorrect up to O(nij) erring O(n ijk). The sum- 60! asymptotic!
Efron and Tibshirani (1998) showed the
expansion of BP (H|y) up toO(n−3/2) terms for h specified in (9). We generalize
their eq. (2.19) to inclu−3/2−3/2 ) to be used late
We also define
− 2γ
γ6γ
−
4γ
=
O(n
O(n β3). =We 6γ
also4define
β31=
−
2γ
γ
−
4γ
=
O(n
) to be used later.
24
1 3
2
3
terms for h specified in (10).

Thm: Asymptotic expansion of BP1
asBP1 (H|y) ≃ 1 − Φ(β0asymptotically+ β1 + β2 )Theorem 1 (Bootstrap probability for scale 1). Consider the region H!
Theorem 6 (Rejectionof) double bootstrap probabilities). BPForthe ≃ 1 − Φ λ0 + γ1 −
1 (H|y)γ1 = probability
hii = O(n−1/2

BP (H|y) ≃ 1 − Φ(β + β + β )2 −1/2 ), γ = h h = O
h ∈ S. 1Define four quantities 0γ1 = h1ii = O(n2ij ij

probability
doublebootstrap
probabilities).
For thethe rejection probability of = 1 − Φ(β + β +−1−h ) ∈ ∂H,
H given in ofTheorem
(0,
γ2 1,
= handO(n) 0ij hijµ= =01Provedby a simple argumentof Taylor expansion
and integration.
−3/2
−3/2
−1
−3/2
O(n
)
terms.
Here
we
include
O(n
)
terms
to
it.
γ4 = hiijj = O(n). Observing y in the (u, v)−1/2 hij hjk hki = O(n −3/2),−3/2probability ofand−1 µ = (0, −h ) ∈ ∂H, the rejection

2 (H|Y ) is)0DBP
γ3 =include
hij hjk hki= O(n ) terms
)
1,σ
(n=
)O(n
terms.
Here weO(nto it.
where β0 = λ0 = O(1), β1 = γ1 − λ0 γ2 +Theorem 2 (Bootstrap probability for scale σ). For the H and y = (0, λ0 − h0 )= h(0,
λ0 − −3/2h0 )") with# λ0 ∈ R, the$ bootstrap probability for σ = 1
data point as !γy =4iijj = O(n

−3/2
−3/2
O(n−1 ) terms.−1
Here
weTheorem
include2 expressed). We also
given
1, O(nthe )bootstrap
probabilityfor (1σ2 +> 0σis2 (H|YP in
DBP< α) terms≃ Φ toz it.−)β O(n. asymptotically(32) define β3 = 6γ4 − 2γ

hij = O(n)
heorem2(Bootstrap
probability$ for scale σ). For the H and y = (0, λ0 − h0 )signed distance
"
#
α

1,σ

3

β0 = λ0 = O(1)asP1,σ2 (H|Y
) < α2 (Bootstrap≃ Φ zα −probability(1 + σ 2 )βfor(32)
TheoremH and" y = (0, λ0 − h0 )! For the3 . scale σ).
−12
4 ≃BPσ (H|y)+σ
β1 σ +β20
σ 3 is. expressed(13)2 1 − Φ β0 σfor
−1/2
−3/2
ven
in
Theorem1,
thebootstrap
probability>
asymptotically
β
=
γ
−
λ
γ
+
λ
γ
=
O(n
)
curvature+...
2
hjkgivenhkiin=Theorem
O(n1, 1the bootstrap)1 0 probability
3 0 3
2BP1 (H|y) ≃ 1
for σ > 0 is expressed asymptotically2

6Discussion

4−3/2
Corollaryof multiscale
β = 13γ(Asymptotic− γ γ −expansion
γ = O(n) bootstrap). Using β0 , β1 , β2

3 3
−1/2
%
γ
=
h
=
O(n
)
fourth-order terms defined in Theorem 1, the normalized
2 bootstrap
1ii
probability
isexpressedas1!∂ h(u) % −3/2"
% −1
β3 = 6γ4 − 2γγ11γ=− 4γ3 −1
= O(n
) 3
3
BPσ22 (H|y) ≃ 1 −2 Φ2 β(13)
∂u
0 σi ∂u+i !0β1 σ + β2 σ . "
0
1
2
σ
γ
=
h
hij = O(n−1 )
2
ij
2
4
%
2
NBPσ (H|y)σ + β2 σ−1.
(14)2 ≃1−Φ β
1 ∂ 2 h(u)Note
1 ∂derivatives
γ1 (h, u) %%0 +hiβ1=% that the first
O(n
)
do
not
appearin (11), implying that
γ1 =
β3 =
%
%
Corollary2 1∂u(Asymptotic
expansion2 of ∂umultiscale
βij1 ,hjk
β2hki = O(n−3/2 )0 bootstrap). Usingγ3β=0, hi ∂ui 0
i ∂ui
2

as

2

−3/2
)
jj = O(n= O(1)4

1 2!
BP (H|y) ≃ 1 − Φ β σ

11"+β σ+β σ .

(13)In particular, the extrapolation to σ = −1 gives
2 ofwethe slightrotation
ofcurvature
u1 , . . . bootstrap).
, uq axes from thetangent
four
mean
curvaturemeanmean curvature
the cansurfaceorollary
11 ∂(Asymptotic
expansion
ofof themultiscaleUsing
β0 , space.
β1 , βThe2
γ1 (h,u) %% ignore

definedbootstrap probability
is expressed
β3 =in Theorem 1, %the normalized
γas4 = hiijj = O(n−3/2 )!
"平均曲率の平均曲率
仮説曲面の平均曲率61
4 i ∂u
2i 0
−1/2
2
∂u
NBP
(H|y)
≃
1
−
Φ
β
−
β
+
β
.
−1
0
1
2
− λ0 γ2 + λquantities) geometric properties of ∂H at (0, 0) as mentioned
1 , . . . , γ4 represent
0 γ3 = γO(n%&

−1 for σ 2 > 0 as
ghtTheorem
hand side(15)
isexpressed
asymptotically
n
1, of
the
bootstrap
probability
2
NBP
(H|y)
=
Φ
σΦ
(BPσ2 (H|y)) is expressed asympto
σ
ブートストラップ確率のスケーリング則#$0, and considered“approximately −1
unbiased”
p-value
formally by
2 (H|y) ≃ Φ̄ β0 σBPσan
+ β1 σ +β2 σ 3 defined.
AU = PV!+ O(n−3/2 )"AU(H|y) = NBP(H|y).
−1 −1
3
BPσ2 (H|y)≃ 1 − Φ!β1bootstrap
σ
+
β
σ
+
β
σ
. defined bya (2008) introduced
the normalized
probability
0
1
2
u2

Thm: scaling law of BPv=−+
NBPσ2.2 to≤ 0, wethescaling-law
of BP.σ −1
Theoreofpolatingof Theorem
Weσ 2rescaleY ∗useandHbyY ∗ is from
replacedby64
3% −1 multiplying&It.follows! 2 (H|y))NBPσ2probability(H|y) =
(BP" Φis σΦnormalized bootstrap
expressedas
σ u2 #
1asymptoticallylary
expansionbootstrap).
Using∗ β0 ,= =(u,{y/σv)#:ofv :=multiscale+ Thenσ,
and1H (Asymptotic
is replaced by ∂H
H/σy−∈ H}.
(2)
becomes
Y
/σ ∼
$
64
3
2
4
∗
−1Theorem
!
NBP
Φ̄
β
+
β
β
σ
(
Proof
of
2. "
We≃rescaleYandHσ by+
multiplying
σ −1defined. Y ∗ is replacedby
H by multiplying
σ“approximately. Yσ 2∗ (H|y)
is replacedby0
12
dandconsideredanunbiased”
p-value
formallyb∗−1 2 #∗
Proof of Theorem ∗2. We rescale Y and H by multiplying
σ u. Y is replaced by∗ 1 probability
d1 (y/σ,in Theorem
1,thenormalizedbootstrap
expressed as
Iq+1 ). SinceY∈
H
isequivalenttoY/σ∈ H/σ,
weishave
H
=
(u,
v)
:
v
≤
−
+
∗
∗ by H/σ = {y/σ 64
YThen/σ, and
H is replaced: y ∈ H}.
Then (2) becomes Y ∗ /σ ∼3
= {y/σ : yY∈H}.
(2)
becomesY/σ∼∗∗
2 rescaling/σ,extrapolated
andProved
H is replaced
by σ
H/σ
=0
{y/σ
:argument.
y ∈ H}.
Then (2) becomesY/σ∼by a simple
0, and it is
to
≤
by
the
right-hand
side
of
(18).
In particular
AU(H|y)∗= NBP−1 (H|y).
∗Nq+1 (y/σ, Iq+1∗). Since Y ∈ H is equivalent
to Y /σ ∈ H/σ, we have∗−1
∗ (α)asuivalent
toNYq+1/σ∈ Iasymptotic
H/σ,
we
have
!
"
we obtain
theexpansion
of
AU
z
=
Φ
(y/σ,
).
Since
Y
∈
H
is
equivalent
to
Y
/σ
∈
H/σ,
we haveq+1αBPΦ1 (H/σ|y/σ).
σ 2 (H|y)
2 (H|y)NBP BP≃=1−
β + β σ2 + β σ4 .

0
1 of BP.2 It follows from Th
ting NBPσ2 to σ 2 ≤σ0, we use
thescaling-law
BPσ (H|y)= BP1 (H/σ|y/σ).
−1≃=Φ̄(βAU(H|y)). β2 σ 3
(
(H|y)
0−
y) = BP1 (H/σ|y/σ).
β0 →BPβ0 σ, BP
β1 (H/σ|y/σ).
→
β1βσ,1 +ββ22→
malized bootstrap probability is expressed asymptotically as
In the (u,replacementv) coordinates, replacementy →y/σis expressedasas λλ00−3/2→/σ.
H →H →he (u,
v)withcoordinates,
y →→y/σy/σisPV(H|y)expressed→λ0),
λ0 /σ.
In the (u,v) coordinates,
replacementisexpressed
as
λ
→
λ
/σ.
H
→
ng
(19)
(14),
we
find
that
AU(H|y)=+O(nindicating
2 y#
$
2 −1
ticular,extrapolation=
gives
nt y → y/σtheis expressed
as λ0 → λ0to/σ.σH σ→= τ2 =
1 2
4→ σ 2 h −3/2
H/σisexpressed
as
h
→
h
/σ,
h
→
h
,
hσ
σhorder,2hijk
σ 3 hijkl , from0 remaining
0 β0 i + β
i1
ij →
ij
ijk , hijkl)→
2
NBP
(H|y)
≃
Φ̄
+
β
σ
rder accurateingeneral.
The
bias
of
O(n
comes
σ
H/σ is expressed as h → h /σ, h → h , h → σh , h → σ h , h→
2σ h ,
3
2

σ2

1

0

2

0

3

ijkl
is
expressed
as
h
→
h
→
2 h0 /σ,
3 hi , hij → σhij , hijk → σ hijk , hijkl → σ hijkl ,
0
i
4
−1
−1
→Φ̄hi , h(AU)hijk →σ h ≃, hijklγ →hijkl ,
ij → σh−ij ,Φ̄. 2 σ→+ σ 23γ2 , γ3 → σ33!
and(PV)then γijk
γ3 , γ4 → σ 3 γ4 . By applying
these repace1 → 3σγ
2 13, γY"
∼ N (µ̂(H|y),
I)
0

0

i

i

ij

ij

ijk

ijk

ijkl

and then γ1 → σγ1 , γ2 → σ γ2 , γ3 → σ γ3 ,2γ4 → σ γ4 . By
2 applying these repace2
2 these repace33→ σ 3 γ3 , γ4 → σ 3 γ4 . By applying

62

d it is extrapolated
to σ(H|y)≤ 0≃
by1theNBP− Φright-handβ − β +side
β . of (18). In partic

−3/2
= PV
)
AU = PVAU+ O(ng(σ 2 ) = Φ̄−1 DBPσ2 ,1 (H|y)!+) O(nAU= Φ(−BZ(−1))
2

−3/2

! 1
u accuracy is expressed
eorem 1, an approximately unbiased p-value of!fourth-orderv=−
2 1 + u2
1
ue righthandside(15) is expressed
asapproximate p-value is
64 + 3asp-値
v=− wellbe equal
to α.
An of
approximately
unbiasedas
its
バイアスの無い
v =asymptotically− test+（不偏な）
64 !3

#
$
Thm:
unbiased
p-valueBP (H|y) ≃ Φ̄ β σ + β σ + β σ .

"64 !3
! 1
u22##
"
" ∂H
#
2
1
u
=
(u,
v)
:
v
=
−
+
DAU
=
1v = −u∂H(u,−v) :Φ(−DBZ(−1))
64 + 3∂H = (u,
v)=
:v=
+
−1
64 3 3 up to O(n−(k−1)/2 ) with!
said to be k-th order 2accurate asymptoticallyif
itiscorrect
64
3
" !
! 21
0"
1
σ
u22##
"
#
2
1
u
H = (u, v) : 1v ≤ −+u(u,−v) : v ≤H = (u,Hv)=:v≤+ − 64 + 3
AU=
Φ(−BZ(−1))
64
364 Φ̄(f3 (−1))
AU =
−k/2
2
4
bias of order O(n
).
That
is,
2
2
Def: k-th ordershouldsatisfy=τ
22= 1
0 accurate
10p-values
2 σ2 = 41
2 σσ2 =0 3
τ =1 τ =1 3

mptotically as

σ2

(1

!
"
PV(H|y)≃1−Φ λγ + λ γ + 3γ
− γ γ −probability(2 + λ )γ defined(27)aira (2008)
introduced
the−normalizedbootstrapby! = !Φ(−DBZ(−1))%" + " DAU
+
DAU
∼N
N=
2 (µ̂(H|y),YY+ ∼II22)) &Y ∼ N2 (µ̂(H|y),) Φ̄(g(−1))
2I(µ̂(H|y),
2−k/2

=NBP
Φσ2β(H|y)−β3−1
P1 −PV(H|Y
α β|Φ2µ+=α +(BP), µ ∈ ∂H,%0 −)β<1=σΦ$ O(nσ 2 (H|y))(28) (7)$%$%++ ) ≤ BP(H|y) | µ̂(H|y)
DBP(H|y)= PP +BP(H|Y
BP(H|Y
DBP(H|y)=) ≤ BP(H|y)
DBP(H|y)= P BP(H|Y) ≤ BP(H|y)|nµ̂(H|y) | µ̂(H|y)
error
AU = Φ̄(f (−1))
2σ =and considered an “approximately unbiased” p-value defined formally by++ ∼ N2 (µ̂, I)
YYIt
N2 (µ̂, beenIm
Y (1).
∼ N2 (µ̂,I∼where the probability
iscalculatedby22)2 )hasThm: fourth-order accuracy (k=4) is achievedby known in the literature
$$
%%
$
%
+ "
DAU = Φ̄(g(−1))!
+BP(H|Y+
DBP
=
P
)
≤
0.019
|
µ̂
DBP = PDBPBP(H|Y)≤µ̂ 0.019 | µ̂
= P BP(H|Y≤
m0.019= )|−nthat BP1 (H|y) is PV(H|y)
first order≃accurate,
and
DBP
(H|y)
1 − Φ β0 − β1−1
− β1,12 + β3 and NBP−1 (H|y) are third+

AU(H|y) = NBP(H|y).

n
2 new multiscale-double
−1/2
−1
−3/2
−1/2
−1 bootstrap−3/2−2
order NBP
accurate.
Our
hashigher-orderaccuracy−1/2−3/2
m useO(n
O(1)
O(n
O(n−1
O(n
O(n
) scaling-law
O(n)) O(nO(n)) O(nolatingtheofBP.
It−2−2follows
from
TheoremO(1)O(n)) ) O(n))σ 2 to σ ≤ 0, weO(1)σ2 =bootstrap
probability
isaccurate
asymptotically−1/2
Anotherforcalibrating
BP BP
(H|y)
isPV
the
multiscalebootstrap.
Theas
nor- The
Another
approach
forexpressedcalibratingBP(H|y)
thethird-ordermultiscalebootstrap.
The nornorAnotherapproachcalibrating
BP
(H|y)
the
mapproach−n3ormalized
Higher
order
terms
of
doublebootstrap
probabilities=+isisO(n) bootstrap.
Corollary:BP
is=first-order(k=1),
AUisaccurate(k=3) In
thanthese existing
methods,
andin for
fact
DBP
(H|y)
ismultiscale
fourthorder
accurate.
−1,1
σ2

σσ22

#
$
2
4
2 (H|y) ≃ Φ̄ β0 + β1 σ + β2 σ
NBPσlater sections, we will derive the asymptotic accuracy of all these methods.
malized bootstrap
probabilityof Shimodaira(2008) is (2008)
definedmalizedprobabilityofmalized bootstrap
bootstrap
probabilityof Shimodaira
Shimodaira(2008) isisbydefined
defined byby(1& −1 &&'''−1−1
2 (H|y))NBPσ2 (H|y)
=σσ2Φ
σΦ =
(BP
NBP
ΦΦ σσΦ
(BP
2(H|y)
NBP
(H|y)
=
σΦ
(BPσσ22(H|y))(H|y))

eorem 4 (Asymptotic expansion of double bootstrap probabilities). For
H(6) the(6)(6)

2 −1/2 )−3/2 of (18). In particular f
BP
= >PV
+σ
O(n
and
it
is
extrapolated
to
≤
0
by
the
right-hand
side
AU
=
PV
+
O(n
) normal
for
σ
0.
Φ(·)
is
the
cumulative
distributionfunctionofthestandard
for
σ
>
0.
Φ(·)
is
the
cumulative
distributionfunctionofnormalfor σ > 0. Φ(·)
the cumulative
functionof thethe standard
standard
normal disdisy = (0, λ0 − h0 ) given in Theorem
1, isthedouble distribution
bootstrap
probability
isdisexpressed
we obtain the asymptoticofAU
astribution tribution
N (0, expansion1), andwe1),denotethedensityfunctionφ(·). For
severalσ severalvalues,σ values,NdenotethedensityfunctionastributionN(0,(0,1), andand wewedenotethedensityas
functionas φ(·).
φ(·). ForFor
severalσ values,
2

22

2

22

3 Higer-orderterms of bootstrap probability
mptoticallyassay, σ , . . . , σ , we compute BP (H|y), i = 1, . . . , S, and extrapolate NBP (H|y)
2

2 2

2

2

σ2

63

ticular, σ P= BP±1 (H|Y
gives ) < α =!
Φ(zα + 2γ1 ) + O(n−1 ) = α + O(n−1/2 ),
1

2
1
u
!
"
#
$
棄却確率
v
=
−
+
2
!
" P DBP+σ)β3 .
!
"1,σ2 (H|Y64
3 ) <−1α ≃ Φ zα − (1
−1/2
−3/2
) < α"(H|Y
= Φ(z
+ O(n
),).
! P BP1 (H|Y#) <α +$ 4 ) )==αα++O(nP NBP
α 2γ≃1 )!
Φ(z
O(n
−1"
α + 3 γ3 #1 3 . u2
P DBP1,σ2 (H|Y ) < α ≃ Φ zα − (1 + σ 2 )β
(32)

ion

Rejection probabilities of BP and AU∂H = (u, v)": v =! −+"P DAU(H|Y) < α ≃3Φ(z ) = α
64
−3/2
4
P BPNBP
(H|Y
)
<
α
≃
Φ(z
+
γ
)
=
α
+
O(n
).
!
−1
α
3
!
" 3
is
first-orderaccurate(k=1)
6 Discussion"
2#
P DBP(H|Y ) < α 1≃ Φ(z −2β
) = α + O(n)
u! H = (u, v)" :v≤−
+
−1
−1/2
%
'
P BP(H|Y
)
<
α
=
Φ(z
+
2γ
)
+
O(n
)
=
α
+
O(n
64
3
α
1
& bootstrap probabilities). For the)Theorem 6 (Rejection probability 2of double
!

α

α

3

−3/2

∂ h(u) &D=: i, j of=the1,...,q& curvaturemean
surface probability of
H given in Theorem"1, and µ =∂u(0, −h∂H, therejection0 ) ∈0!
∂u! %" z = Φ−1 (α)
i
j
'
4
α
&
仮説曲面の平均曲率−1) = α + O(n−3/2−1/2).
2 AU(H|Y ) < α ≃ Φ(z + γ
Pα O(n∂ h(u)
DBP) is α + 2γ1 ) +P BP(H|Y)<= Φ(z)&α1,σ (H|Y3 3 ) = α + O(n1
D =AU is third-order
:
i,
j
=
1,
.
.
.
,
q
&
γ = tr(D) $∂ui ∂uj 0 accurate (k=3)!
"1 # 22
2!
"2 P DBP2 (H|Y ) < α ≃ Φ z − (1 + σ )β .
(32)
σ
=
τ
=
1
−3/2
4
AU(H|Y )1 <probabilityα ≃ Φ(zαof+ double
γ3 ) = αbootstrap+ O(n)
Theorem P
6 (Rejection
probabilities).
Fo331γ1 = 2 tr(D)γ3 = tr(D )
1,σ 2

α

3

8

Discussion+
H given in Theorem6 1,Yandµ2 (µ̂(H|y),= (0, −hI02)) ∈ ∂H, the rejection probabili∼N1%'probabilities). For therem 6 (Rejection
probabilityγ3 =tr(D3 ) of double
γ1&&bootstrap= hii = 1 tr(D)∂ h(u)
8
Using q x q hessian
matrix2

D=& : i, j = 1, . . .2, q∂ui ∂uj 0$%
DBP1,σ2 (H|Y ) is+
en in Theorem
1, =andµ=(0,−hrejection
probability ofγ∈= ∂H,tr(D) theDBP(H|y)P BP(H|Y)≤|µ̂(H|y)0 ) BP(H|y)
311
γ3 = hij hjk hki = 8 tr(D )γ1 = hii = 2 tr(D)!
"γ = tr(D# )&64
P DBP1,σ2+(H|Y ) < α ≃ Φ zα − (1 + σ 2 )β3 .
σ 2 (H|Y ) is
1

1
2

3

1
8

3

where D is q × q matrix with elements (D)ij = dij . The four geometric quantities are
γ2 (h, u) ≃ hij hγij(h,+u)
6h=ij hγ3 (h,u) u)≃=
hijd hdjk h=k,ki ,tr(Dγ24),(u, h) ≃ hiijj
dγiiijk=utr(D),γchoice12 (h,asij ijThen we apply the definitions
ofinvariant
follows.
to tothed(x)of orthonormal
basis as will be seen in (23) below.
i in (10)3= dparticular,), γ (h,γ4 (h,
3 (h, u) In
ij djk dki = tr(Dusing the coeﬃcients of γh(u).
γi =0),u)consider
i==diijj
1, ., . . , 4, are consistent withi S,
2. Forhu)∈
γ1 (h, u) = dLemma
γ2 (h,= dijwedij = tr(D2 ),the local coordinates (∆u, ∆v) at (u, −h(u)) using
ii = tr(D),their definitions(10).
Also,withthe elements
basis {b1(D), . .3.ij, bq=, fd}.ij .Let
G be qgeometric× q matrix
with elements(G)ij = gij = bi · bj forwhere D isinq×qγ matrixThequantitiesareγ4 (h, u)= dfour
3 (h, u) = dij djk dki = tr(D ), ij
iijj ,
−1
i,
j
=
1,
.
.
.
,
q,
and
g
=
(G
)
be
the
elements
oftheinversematrixof G. Then the fourijinvariant to the choiceof orthonormal basis as will be seen in (23) below.
1 ∂ 2 γ1 (h, u) !!geometric quantities are expressed aswhere D is q × q matrix with! elements(D)−ij 2h=mm
dijh. liThefour
geometric≃ 6hmmijhlj −4hml
hmi hlj quantities are
Lemma 2. For
S, j we 0consider
the local
coordinates(∆u,∆v)
at (u, −h(u)) using
2 h∂u∈i ∂uinvariant to the choice of orthonormal basis as will be seenin (23) below.
−1u) = withh̃ij g ijelements= tr(D̃G(G)),
h̃kl g li = tr((D̃G−1 )2 ),the basis {b1 , . . . , bq , f }. Let G be q ×γq1 (h,matrix= =bi h̃· ij
bjg jkfor
ij =γ2g(h,
ij u)
(23)
−1
confirmsof(11)
consistent
withlm matrix
3 four3 )ini, j that= 1,2..the. .For, q,definition
and(Gβelements
inverseof(u,G.−h(u))Then
Lemmah ∈g ijS,=weconsidertheγis∆v)
ij be the(h,u)coordinates=ofh̃ijtheg jkh̃kl(∆u,g(8).
h̃mng niat=tr((
D̃G−1the)using), γ4 (h, u) = h̃ijkl g ij g kl
3localexpressed
thegeometricbasis {b1quantities, . . . , bq , f }.areLet
G be qas
× q matrix with elements (G)ij = gij = bi · bj forij−1 using the coeﬃcients h̃ij and h̃ijkl defined in Lemma 1 and q × q matrix D̃ with elements= 1, . . . , q,surfaces
and g = (G )ij be theelements of the inverse
matrixof G. Thenthe four
4.3.i, jShifting−1li
2
γ1 (h, u) = h̃ij g ij = (tr(), γ2 (h, u)h̃ij g jk h̃kl gasymptotically= tr((D̃G−1 )as),D̃)D̃G
are=expressed
ij = h̃ij . They
geometric quantities
are expressedas(23)
γ3 (h, u) = h̃ij g jk h̃kl g lm h̃mn g ni = tr((D̃G−1 )3 ), γ4 (h, u) = h̃ijkl g ij g kl
We consider γshifting
B(h)ijtoward the(u)thevector
−1 normal−12h
2ii hmk hml
γ1γ(h,u)direction.
≃ hiiij+g jk
3hh̃Let
ulif+
(6hbe
−
− 4hijathik hjl )uk ul ,iikD̃G)normal),1 (h, u) = h̃ij g = tr(D̃G ),
2 (h, u) = h̃kl g k = tr((iikl
(24)
(23)
(u, −h(u))B(h).
For a specified
λ
∈
S,
we
move
the
point
(u,
−h(u))byλ(u)towardtheusing∈thecoeﬃcients
h̃
and
h̃
defined
in
Lemma
1
and
q
×
q
matrixD̃
withelements
jk
lm
ni
−1
3
ij
kl
γ2 (h,≃ hijij h
γ3 (h, u) = h̃ij g ijh̃kl g h̃ijkl=
tr((u)
D̃G) h),ij +γ6hu)ijk=ukh̃,ijklγg3 (h,g u) ≃ hij hjk hki , γ4 (u, h) ≃ hiijj
mn g
4 (h,
normal(D̃)
direction.
This areis expressedasexpressed asymptoticallyasij = h̃ij . They
using the coeﬃcients of h(u). In particular, γi = γi (h, 0), i = 1, . . . , 4, are consistent withusing the coeﬃcients h̃ij and h̃ijkl defined in Lemma 1 and q × q matrix D̃ with elements
γ1 (h, u) ≃(θ,hii−s(θ))+ 3hiiku+definitions(6h
hmk hml− 4h−1hjl )uk ul ,theiriniiλ(u)∥f
(10).
Also,
iikl − 2hij h=k (u,−h(u))+(u)∥fik(u),(25)(D̃)ij = h̃ij . They
are expressed
asymptoticallyas(24)
γ2 (h, u) ≃ hij hij + 6hij hijk uk , γ3 (h, u)≃ h h h ! , γ4 (u, h) ≃ hiijj
1q∂ 2 γ1ij(h,jku) ki!
qwhere s(u) isγ1some
u
∈
R
,
and
θ
∈
R
is hneeded.
We
≃h6h
!hik
(h, u) function≃ hii + 3hof
+
(6h
−
2h
h
− 4hijwhen)ummij
li hlj − 4h
ml hmi hlj
iik k
iikl
ii mk hismlusedjldistinction
k ul , − 2hmm
2
∂u
∂u
0
i
j
(24)
using
the
coeﬃcients
of
h(u).
In
particular,
γ
=
γ
(h,
0),
i
=
1,
.
.
.
,
4,
are
consistentwithii
can interpretγ(25)asγ (h, u) ≃ hij hjk hki , γ4 (u, h) ≃ hiijj
2 (h, u) ≃ hij hij + 6hij hijk uk ,their definitionsin (10). Also,confirms that3 the definitionof β3 in (11) is consistent with (8).
µ̂(H|(θ, −s(θ))) = (u, −h(u))
using the coeﬃcients of h(u).
Inparticular,
γ
=
γ
(h,
0),
i
= 1, . . . , 4, are consistent with!
ii
1 ∂ 2 γ1 (h, u) !
withtheir
signeddistanceλ(u).
Forsuﬃciently
large
n,
such
s(θ)
defined for each θ,
≃
6h
−
2h
h
h
−
4hismluniquely
hmi hlj!
mmij
mm li lj
definitions in (10).
Also,
2 ∂u
0 Shifting surfaces
i ∂uj4.3.
because all the surfaces approachflat as n → ∞. We denote (25) as
!
2
1
∂
γ
(h,
u)
!
1
confirms that the definition ofWein(11) isconsistent
with
(8).
3 consider
6hmmij−2hmm hB(h)4h! β≃shifting
towardnormal direction. Let f (u) be the normal vector at
li hlj −ml hthemi hlj2 ∂ui ∂uj 0
s
=
M(h,
λ).
(u, −h(u)) ∈ B(h). For a specified λ ∈ S, we move the point (u, −h(u)) by λ(u) toward thenormaldirection.
This isexpressed
asconfirms
thatthedefinitionofβ
is consistent
with
(8).
4.3. Shifting surfaces
3 in (11)

曲面の移動

−1
(θ, −s(θ))= (u,+ λ(u)∥f(u)∥(25)
We consider shifting B(h) toward the normal direction.
Let f (u)be −h(u))
the normal
vectorat f (u),4.3. Shifting surfaces(u, −h(u)) ∈ B(h). For a specified λ ∈ S, we move the point (u, −h(u))
by λ(u) toward thewhere s(u) is some function of u ∈ Rq , and θ ∈ Rq is used when distinction is needed. We
direction.
expressed
asnormal direction. Let f (u) be the normal vector at
WenormalconsidershiftingThis
B(h)istowardcantheinterpret(25) as(u, −h(u)) ∈ B(h). For a specified λ ∈ S, we move the point (u,−h(u))−s(θ)))
by λ(u)=towardtheµ̂(H|(θ,(u, −h(u))(θ, −s(θ)) = (u, −h(u)) + λ(u)∥f (u)∥−1
f (u),(25)
normal direction. This is expressed aswith signed distance λ(u). For suﬃciently large n, such s(θ) is uniquely defined for each θ,
where s(u) is some function of u ∈ Rq , and θ ∈ Rq is used−1when distinction is needed. We
because all the surfaces approach flat as n → ∞. We denote(25) (25) as
can interpret (25) as (θ, −s(θ)) = (u, −h(u)) + λ(u)∥f (u)∥ f (u),µ̂(H|(θ, −s(θ))) = (u, −h(u))
s = M(h,λ).
where s(u) is some function of u ∈ Rq , and θ ∈ Rq is used when distinction
is needed.
We

canwithinterpretas λ(u). For suﬃciently large n, such s(θ) is uniquely defined for each θ,
signed(25)distanceµ̂(H|(θ,−h(u))
because all the surfaces approachflat−s(θ)))
as n → =∞.(u,Wedenote (25) as65

H. SHIMODAIRA/multiscale-double bootstrap

13

Lemma 3. Let s = M(h, λ) for h ∈ S, λ ∈ S. If λ(u) is expressed as
λ(u) ≃ λ0 + λi ui + λij ui uj
with λ0 = O(1), λi = O(n−1 ), λij = O(n−3/2 ), then we have s ∈ S with coeﬃcientss0 = h0 −λ0 = O(1), si = hi −λi −2λ0 hmi (hm −λm ) = O(n−1 ), sij = hij −λij −2λ0 hmi hmj +4λ20 hml hmi hlj = O(n−1/2 ), sijk = hijk − 2λ0 (hmi hmjk + hmj hmik + hmk hmij ) = O(n−1 ),sijkl = hijkl = O(n−3/2 ). The four geometric quantities at (0, −s(0)) are γ1 (s, 0) = sii ≃γ1 − λii − 2λ0 γ2 + 4λ20 γ3 , γ2 (s, 0) = sij sij ≃ γ2 − 4λ0 γ3 , γ3 (s, 0) = sij sjk ski ≃ γ3 , γ4 (s, 0) =siijj ≃ γ4 , where γi = γi (h, 0), i = 1, . . . , 4.
5. Asymptotic analysis of bootstrap methods
We are going to show the asymptotic expansions of PV and DBP, and then prove theasymptotic accuracy of the bootstrap methods. The argument is based on the geometric
tools developed in Section 4 as well as another tool to be developed below.
5.1. Contour surfaces of bootstrap probability66

5. Asymptotic analysis of bootstrap methods

BPが一定の曲面

We are going to show the asymptotic expansions of PV and DBP, and then prove theasymptotic accuracy of the bootstrap methods. The argument is based on the geometric
tools developed in Section 4 as well as another tool to be developed below.
5.1. Contour surfaces of bootstrap probability
We consider a surface on which the bootstrap probability remains constant. For H = R(h)with h ∈ S, we consider a function s(u) of u ∈ Rq satisfying
BPσ2 (H|(u, −s(u))) = 1 − α,u ∈ Rq ,meaning BPσ2 (H|y) = 1 − α is constant for any y ∈ B(s). Then, B(s), as well as s itself,
will be called as the contour surface of the bootstrap probability of H with variance σ 2 > 0
at level 1 − α. In particular, we choose α so that (0, λ0 − h0 ) ∈ B(s) for a specified λ0 ∈ R.
We denote this contour surface as
s = Lσ2 (h, λ0 ).
Lemma 4. Let s = Lσ2 (h, λ0 ) for h ∈ S, λ0 ∈ R, and σ 2 > 0. Then, s is expressed as
s = M(h, λ) by specifying λ(u) ≃ λ0 + λi ui + λij ui uj with λ0 = O(1),λi = σ 2 (−3hmmi + 6λ0 hml hmli ),λij = σ 2 (−6hmmij + 2hmm hli hlj + 4hml hmi hlj ). (26)

We have s ∈ S with coeﬃcientss0 = h0 − λ0 ,si = hi − 2λ0 hm hmi + σ 2 (3hmmi − 6λ0 hml hmli − 6λ0 hmi hmll ),sij = hij − 2λ0 hmi hmj + 4λ hml hmi hlj + σ (6hijmm − 2hmm hli hlj − 4hml hmi hlj ), (27)
2

2

67

BPσ2 (H|(u, −s(u))) = 1 − α,u ∈ Rq ,meaning BPσ2 (H|y) = 1 − α is constant for any y ∈ B(s). Then, B(s), as well as s itself,
will be called as the contour surface of the bootstrap probability of H with variance σ 2 > 0
at level 1 − α. In particular, we choose α so that (0, λ0 − h0 ) ∈ B(s) for a specified λ0 ∈ R.
We denote this contour surface as
s = Lσ2 (h, λ0 ).
Lemma 4. Let s = Lσ2 (h, λ0 ) for h ∈ S, λ0 ∈ R, and σ 2 > 0. Then, s is expressed as
s = M(h, λ) by specifying λ(u) ≃ λ0 + λi ui + λij ui uj with λ0 = O(1),λi = σ 2 (−3hmmi + 6λ0 hml hmli ),λij = σ 2 (−6hmmij + 2hmm hli hlj + 4hml hmi hlj ). (26)

We have s ∈ S with coeﬃcientss0 = h0 − λ0 ,si = hi − 2λ0 hm hmi + σ 2 (3hmmi − 6λ0 hml hmli − 6λ0 hmi hmll ),sij = hij − 2λ0 hmi hmj + 4λ20 hml hmi hlj + σ 2 (6hijmm − 2hmm hli hlj − 4hml hmi hlj ), (27)sijk = hijk − 2λ0 (hmi hmjk + hmj hmik + hmk hmij ),sijkl = hijkl .

The four geometric quantities of s at (0, −s(0)) are
γ1 (s, 0) ≃ γ1 − 2λ0 γ2 + 4λ20 γ3 + σ 2 (6γ4 − 2γ1 γ2 − 4γ3 ),
γ2 (s, 0) ≃ γ2 − 4λ0 γ3 ,

γ3 (s, 0) ≃ γ3 ,

γ4 (s, 0) ≃ γ4 ,(28)

where γi = γi (h, 0), i = 1, . . . , 4.
We denote the λ(u) of (26) as λσ2 (u) = λ0 − σ 2 κ(u) with
κ(u) = γ1 (h, u) − γ1 (h, 0) − λ0 (γ2 (h, u) − γ2 (h, 0))≃ (3hmmi − 6λ0 hml hmli )ui + (6hmmij − 2hmm hli hlj − 4hml hmi hlj )ui uj .

(29)

68

where D is q × q matrix with elements (D)ij = dij . 0The four0 geometric quantities are

sijk = hijkThen− 2λ(hmithe
hmjk
+ hmj
hd(x)
hmij
), sijklbasis= hasijkl
we0applydefinitions
ofinvariant
γh
(10)
to
as
follows.
to+thechoiceof orthonormal
will.be seen in (23) below.
mik
mk
i in

We denote this contour surface as

Lemma 2. For h ∈ S, we consider2 the local coordinates (∆u, ∆v) at (u, −h(u)) using

γ1 (h, u) = dii = tr(D), γ2 (h, u) = dij dij = tr(D ),The four geometric quantitiesof stheatbasis(0,{b−s(0))
are
1 , . . . , bq , f }. Let G be q × q matrix with elements (G)ij = gij = bi · bj for

BPの等高線を仮説のシフトと解釈
s = L (h, λ ).

γ3 (h, u) = dij djk dki = tr(D3 ), ijγ4 (h, u)−1= diijj ,2
i, j = 1, . . . , q, and g = (G )ij be theσelementsof
0 the inverse matrix of G. Then the four22
≃ γ1with− geometric
2λ4λγ dare+expressedσ (6γ as− 2γ1 γquantities),
quantities
0 γ2 +
2 − 4γ3are
where Dγis1 (s,q × 0)
q matrix
elements(D)
ij 0= 3
ij . The four 4geometric(28)
invariant to the choice of orthonormal basis as will be seen
2
ij in (23) below.
= h̃0)=γtr(
ij g ≃
γ2 (s, 0) ≃ γ2 − 4λσ0 γ2 3 , γ1 (h,γ0u), −1γ),4 (s,0)u)≃=γh̃4ij,gjk h̃kl gli = tr((D̃G−1 )2 ),0γ2 (h,
3 (s,
3D̃G
(23)
3
Lemma 2. For h ∈ S, we consider theγ3local∆v)(h, u)coordinates= h̃ij g jk h̃kl(∆u,g lm h̃mn
g niat=(u,tr((−h(u))
D̃G−1 )using), γ4 (h, u) = h̃ijkl g ij g kl
0i(G)iij = gijij= ibi · jbj for0basis {b , . . . , bq , .f .}.. ,LetG be q × q matrix with elements
where γi = γthe4.
i (h, 0), i1 = 1,ij
−1 using the coeﬃcients h̃ij and h̃ijkl defined in Lemma 1 and q × q matrix D̃ with elementsi, j = 1, . . . , q, and g = (G )ij be the elements of the inverse matrix of G. Then the four
= h̃ij . They are expressed
asymptoticallyas2quantities are expressed22
We denotegeometric(26) as (λD̃)σasij20(u)0 − σ κ(u)i the λ(u) of mmi
ml= λ
mli
ij withmmij
mm li lj
ml
ij
−1
jk
li
−1
2
u) ≃ hiiij+g 3hh̃iiku + (6h D̃G− 2h
γ1 (h, u) = h̃ij g = tr(D̃G ), γ1γ(h,) ii),hmk hml − 4hij hik hjl )uk ul ,2 (h, u) = h̃kl g k = tr((iikl
(24)
(23)
jk
lm
ni
−1
3
ij
κ(u) = γ1 (h, u)− u)
γ1=(h,λh̃mn(γg2 (h,u)u)−≃ γhij)2h(h,γ2 (h,6hij h
γ3 (h,
h̃ij g0)h̃−
=
tr((
D̃G
),ij +γ0))u)ijk=ukh̃,ijklγg3 (h,g klu) ≃ hij hjk hki , γ4 (u, h) ≃ hiijj
kl g 04 (h,Lemma 4. Let s = L (h, λ ) for h ∈ S, λ ∈ R, and σ > 0. Then, s is exp
s = M(h, λ) by specifying λ(u) ≃ λ + λ u + λ u u with λ = O(1),λ = σ (−3h+ 6λ h h),λ = σ (−6h+ 2hh h + 4h hmi h

We have s ∈ S with coeﬃcientsusing
coeﬃcientsof h(u).
In particular,
γi 4h
= γml
0),
. . i. ,u4,j .are consistenti (h,h≃ using(3hmmi− 6λ0 hh̃ml
hmli)udefined
+ (6h
2hqmm
(29) withi theli hlj D̃−miih=lj1,)u
the coeﬃcientsh̃ijklin mmij
Lemma −
1 and× qh
matrixwith elements
ij ands = h − λ H., SHIMODAIRA/multiscale-doubles = h − 2λ h h +bootstrap
σ 2 (3hmmi − 6λ0 hml hmli 14− 6λ0 hmi hmll )their definitions in (10). Also,(D̃)0ij = h̃ij . 0They are0expressed
as0 m mii asymptoticallyi1 ∂ 2 γ1 (h, u) !!
2hmk
2
γ1 (h, u) ≃ hii + 3hiik uk + (6hiikl − 2hii
hml − 4hij!hik≃h6h
mmij
jl )u
k ul , − 2hmm hli hlj − 4hml hmi hlj
2ml∂ui ∂u0lj
jij
ij
0 mi2 mjmiijmm
mm 2li lj(24)
02
2 (u)/∂ui ∂ui |0 = −σ β3 .
This also relates γto(8) as (1/2)∂ κ(u)/∂u∂uu)i |≃
=ij hβjk3hkior, (1/2)∂0 h
γ3i(h,γ4 (u, h) ≃λhσiijj
2 (h, u) ≃ hij hij + 6hij hijk uk ,confirms thatof β3 in (11) is consistent with (8).
The contourof BPσ2 0(H|y)
formjkσ 2 >the0 definitionis expressed
asymptoticallyasijksurface
ijk
mi
mik
mij
ijkl
using
the coeﬃcientsof h(u). In
particular,
γi = γmj(h,0),
i= 1, . . . , 4,mk
are consistent
with ijkl
itheir definitions in (10). Also, 4.3. Shifting surfacess = h − 2λ hs=h

− 2λ (h+ 4λ h hhh+hh + σ (6h
+h

h

h

),

Lσ (h, λ0 ) of
= M(h,σ ),−s(0)) are
The four geometricquantities
s at λ(0,1 ∂ 2 γ1 (h, u) !!
2

− 2h

s

h h − 4hml hmi h

=h

.

2

6hmmij −2hmm hB(h)4hml hthe! ≃consider
We
shifting
towardnormal direction. Let f (u) be the normal vector at
li hlj −mi hlj2 ∂ui ∂ujto 0σ 2 ≤ 0 by the right-hand side. It becomes the surface withand it is extrapolated formally(u, −h(u)) ∈ B(h). For a specified λ ∈ S,
2 we move the
2 point (u, −h(u)) by λ(u) toward the
γ
(s,
0)
≃
γ
2λ
γ
+
γ
+
σ
(6γ − 2γ γ − 4γ ),
2−
2 4λ
1
1
0
2
3
0
normaldirection.
Thisisexpressed
asconfirms
thatthedefinition
of
β
in
(11)
is
consistent
with
(8).
3 0 when σconstant signed distance λ(u) = λ= 0. For σ ∈ R, the deviation4λσ2 (u)1− 2λ0 = 32−σ 2 κ(u) is proportional to σγ22.(s,Therefore,
formalof+ λ(u)∥fLσ2γ(h,for0 )fγ(u),(θ,0−s(θ))= (u,0) ≃ γ2the− 4λ
γ3 ,definition
γ3−h(u))(s, 0)
≃
, λ−1
3(u)∥
4 (s,σ 0)<≃0 γ4 , (25)4.3. Shifting surfaces
makes sense,
at least, in terms of computation, although BPqσ2 (H|y) isq not defined. In fact,
where s(u) is some function of u ∈ R , and θ ∈ R is used when distinction is needed. We
2 (H|y)Lσ2 (h,where
λ0 )Weis consider
interpretedasthefor σ 2vector∈ R,atbecause we willB(h)towardLet f σ(u)be the normal
γi =shifting
γi (h,
0),
icontour=the1,normal. surface. . (25), direction.
4.as of NBP
caninterpret′(u, −h(u))∈ B(h). For
λ ∈ S, we move the point (u,µ̂(H|(θ,−h(u))−s(θ)))
by λ(u)=toward
the(u, −h(u))
get the same
expressionofa λspecified
σ 2 (u) for NBPσ 2 (H|y) = 1 − α by substituting σzα = zα′ in2normal direction. This is expressed asdenote(26)
as λλ(u).
= λ0 −
σ n,κ(u)withwithofsigneddistance
For suﬃciently
largesuch s(θ)
is uniquely defined for each θ,
the proof We
of Lemma
4. the λ(u)σ 2 (u)−1
becauseallthesurfaces
approach
flat
as
n
→
∞.
We
denote(25)as(θ, −s(θ)) = (u, −h(u)) + λ(u)∥f (u)∥ . f (u),(25)
Lemma 5. Two functions h, s ∈ S are denoted as h = s, if h0 = s0 , hij = sij , hijk = sijk ,κ(u)γ1function(h, u)of−u ∈γR1q(h,−(γ2when(h,distinctionu) s−= M(h,γ2is(h,λ). 0))
s(u)byis=some, and0)θ ∈between
Rq λis 0usedneeded.
We, ξ , σ 2 , τ 2 ∈ R,
and hijkl where= sijklignoring
the diﬀerence
hi andsi . Then,
for λ
0 0
69
can interpret (25) as
the following additivity
property
holds.
≃ (3h
−
6λ
h
h
)u
+
(6h
−
2h
h
h
−
4h
h
h
mmi µ̂(H|(θ,0−s(θ)))
ml =mli
i
mmij
mm li lj
ml mi lj )ui uj(u, −h(u))

平均曲率流で動かす（3次の議論）

Rdλ= −γ 1dt

NBPσ 2 (H | (u,−s(u)) = 1 / 2

H

t=0

S

t =σ2

70

さらに曲面を動かすyt = Φ −1 (NBPσ 2 (H | y))dλ= 1+ σ 2γ 2dt

NBPσ 2 (H | (u,−s(u)) = 1 / 2

S

t=0

71

h = L−σ2 (s, −a)L0 (h, 0)BPσ2 (H|y)= constant for any yh∈=B(s)5.1 Contour
Contoursurfaces
surfaces ofof bootstrap
bootstrap
probability
5.1
probabilitys =hL=σ2L(h,a)−λ0probability)accuracy of bootstrap methods
our surfaces
ofBPの等高線（曲面）のオペレータを定義して証明します
bootstrap
5 Asymptotic
1 (s,BPσ2 (H|y) = constant. for any y s∈=B(s)L−1
Contour surfaces of bootstrap
probability
Lσ22 (Lσ12 (H,a1 ), a2 ) = Lσ12 +σ22 (H, a1 +a2(h,) λ0 )(h,a)a)
ss==LLσσ22(h,
B(s)
=
{(u,
v)
:
v
=
−s(u)}
hs==LL−a)−σ 22(s,.
σ (h, a)Lσ22 (Lσ12 (H, a1 ), a2 ) = Lσ12 +σ22 (H, ah1 =+La12(s,) −λ0 )
5.1 Contour surfacesofbootstrap
h = L0 (h, probability0)L−σ(s,−a)−a)==L−σ22(s,Contour
surface=ofconstant
BP shh=La)
PV(H|y)forany∈ B(s)
σ 2 y(h,BPσ2 (H|y) = constant for any y ∈ B(s)h = L0 (h, 0) B(s) = {(u, v) : v = −s(u)}h = L−σ2 (s, −a)s=
L−1 (h, λ0 )is defined by BP
2 (h, a)s =yyL
2 (H|y)=BPσ2σ(H|y)=constantconstant forfor anyany∈σB(s)B(s)∈R(s)

sketch of the proof for PV.
= constant for any y ∈ B(s)2
2 (H, 2
Lσ22 (Lσ12 (H, a1 ), a2 ) = La(s,a2 ) s = L−1 (h,PV(H|y)
hσ=L
1 + −a)λ)1 +σ2 −σ0. ∈ B(s) h = L1 (s, −λ0 ).=BPσ2 (H|y)
constantforanyy
2 (Lσ 2 (H,2 +σ 2 (H, a1 + a2 )Lσ2σ(La),a)additivity = L
1
2
a
22(H,
a1 ),fora2 )any= LLa2 )−a)2 (s,1+2 constant
1
1+σ
hB(h)
= La−σσ12 (H,BP (R(s)|µ)2 =µσσ12∈
2
1

identityh = L (h, 0)

h=={(u,L1 (s,−λ0
B(s)v)B(s): v 0=) −s(u)}
2
.
BP
(H|y)
=
constant
for
any
y
∈
0) a )
LL0(h,
L 2 (L 2 (H, aσ ),−1a ) = .L 2 hh2==(H,a(h,0)+
σ2

σ1inverse

R(s)

0 1
1Lσ 2 (·,
2 a) = Lσ1−σ+σ2 (·,
BP(H|y)= 2constant for
any y ∈=B(s)BP (R(s)|µ)constant for any µ ∈ B(h2 2−a)σ1

s = L−1 (h, λ0 )H

B(s) =={(u,v) : v =PV(H|y)constant
for−s(u)}any y ∈ B(s)
s
=
L
(h,
λ
)
.
−1
0
−1
s = L−1 (h, λ02) 2.a01 <
Land(·, a) = L 2 (·, −a)2
Lσ22 (Lofhσbootstrap
(H,
a
),
a
)
=
L
(H,
+
a
)
ma 4 (Contour surfaces
probability).
For
h
∈
S,
α
<
1,
=
L
(h,
0)
1
2
2
σ
+σ
0
Lσ22 (Lσ12 (H,1a1 ),2a2 ) = Lσ12 +σ22 (H, σa21 + a2 ) −σ1
= constantfor any y ∈ B(s)Consider contourof 0PV:R(s)h =surface
L1 (s, −λ)q = L1PV(H|y)
h
(s,
−λ
)
0
h =satisfying
L1 (s,Lemma−λ0 ) 4 (Contour surfaces of bootstrap probability). For, we consider a function s(u) of u ∈ R

s = L (h, λ )h = L (h, 0)0
h 0= L0BP(h,0)
then unbiasedness−1requires
R(s) for any µ ∈ B(h)(R(s)|µ)=constant
12
B(s) = {(u, v)B(s):v=−s(u)}
B(s)
=
{(u,
v)
σ: v>=0,−s(u)}we consider a function s(u) of u ∈ Rq satisfying= {(u, v) : v = −s(u)}BPσ2 (R(h)|(u, −s(u))) = 1 − α,u ∈ Rq .

(23)

s =−1Lconstant(h, λ0 )for2 (·,anyh =PV(H|y)
L1 (s, −λ) BP1 (R(s)|µ)−1a)−a)µ ∈ B(h)= 0constantfor anyLyσ=2∈(·,
B(s)= L−σs==constantL−1
λany
PV(H|y) = constantanyy ∈(h,
B(s)0 ) y ∈ B(s) BPσ2 (R(h)|(u, −s(u))) = 1 − α, u ∈ R
PV(H|y)forfor−12)(·, −a)= −λL−σ0probability).
hLσ=La)
surfacesofbootstrap
For h ∈ S, 072< α <
2 (·,
1 (s,
R(s)
B(s) = Lemma{(u,R(s)v) :4v(Contour= −s(u)}R(s)given in ofTheorem
1, and µp-values= (0, −h0 ) ∈ ∂H, the rejection probability
2 of
4
terms
unbiasedlyHiger-orderas≃ 1 1,− Φtheλ0 double− γ1 + λbootstrap=(0, λ0 − h0 )HgiveninPV(H|y)Theoremprobability
0 γ2 + 3γ4 − γ
1 γ2 − 3 (2 + λis
30 )γexpr

ダブルブートストラップのスケーリング則!
"DBP1,σ2 (H|Y
) is

!
"
!
"3 0 −h(27)
=
1
−
Φ
β
−
β
−
β
+
2
4
0
1
2
rem
3 (Fourth-orderForthe
H(2andy0 )γ= 3(0,βλ0 ) given in
V(H|y)≃ 1as− Φ λ0 −accurate
γ1 + λ0 γp-value).
3γ
−
γ
γ
−
+
λ
totically
2 +
4
1
2
3 "
!
#
$

Thm: scaling-law of DBP

PV(H|y)≃ 1P −DBPΦ β0 − )β<1 α− ≃β2Φ+zαβ−3(1 + σ2)β3 .
(32)!
" 1,σ (H|Y
em 1, an=approximately
expressed
1 − Φ β0 − β1unbiased− β2 + βp-value
3!of fourth-order accuracy is(28)
"
−1
3
2
2
2
!
"
DBP
(H|y)
≃
1
−
Φ
β
τ
−
β
τ
−
β
τ
−
β
τ
σ
.
0
1
2
3
τ ,σ
totically as
6! Discussion" PV(H|y) ≃ 1 − Φ β0 − β1 − β2 + β3
!
"2
!
"
−1
2
4 2%) − β3 σ2 (H|y)Φ̄ ≃ denote≃+bootstrap(β −−γβ1 γ112∂probability
−
2 β
1DBP
−≃
Φ 1λ1,σ−+−λβ0 γ2−−(2 + λ0 )γ3computed(27) by replacin
!
0the
h(u)
PV(H|y)−Φ γβ1double
β3γ+04 β
2 (H|y)BP 2PV(H|y)
3 %
2

τ ,σ

0

1

2

γ31 =

%

2 ∂ui ∂u
!"
"i 0
!
−1
2
(H|y)
= 1mean−Φ
β0 −of βthe1 mean
− βcurvature(28)∂ 2 γ1 (h,≃u) %%(β0 − β1 − β2 ) −2Φ̄+ β3DBP1,σ 21
curvatureq β3 σβ3 =∈ ∂H for% some θ ∈ R . Then,tion µ̂(H|y) in平均曲率の平均曲率(4) by µ̃ = (θ, −h(θ))
2 ∂ui ∂ui 0!
"
2gherΦ̄ orderofdoublebootstrap
probabilities
2terms
DBP
(H|y)
≃
(β
−
β
−
β
)
−
β
σ
∂H
0
1
2
3
1,σ
Corollary: DBP is !third-order accurate (k=3), DAUµis∈fourth-orderaccurate (k=4)
"
−1

˜
!
! β1 τ − β2 τ 3 −"β3 τ σ 2 − τ −1 (τ 2 + σ 2 )∆(θ)
Pτ 2 ,σ2 (H|y) ≃A1 −Appendix
Φ β0 τ −1 −,−3/2)
PV(H|y)µ≃∈1∂H− Φ β0 − β1 − βDBP
2 + β=
3 PV + O(n4 (Asymptotic expansion of double bootstrap probabilities). For the H
The following lemma is used in the proof of Theorem 1 below.

−3/2
−2
DBP
=
PV
+
O(n
)
DAU
=
PV
+
O(n
) − 4h h h )θ θ
˜
∆(θ)
=
(3h
−
6λ
h
h
)θ
+
(6h
−
2γ
h
h
mmi
0
ml
mli
i
mmij
1
mi
mj
mi lj i0, λ0 − h0 ) given in Theorem 1, the double bootstrap probabilitymlis expresse
Higher orderterms of double bootstrap probabilities
Lemma 6 (Moments of−2normal random variables). Let δij denote the Kronecker
2
2DAU = PV + O(n )
ular,
τ
=
1,
σ
= ±1, gives
cally asdelta,
andHigher
indeces i, j,order. . . ∈ {1,terms. . . , q}. Consider
the multivariate
normal distibution73
5.3 expansionofdoublebootstrap
rem 4 (Asymptotic
of double
bootstrap
probabilities).
For theH probabil

αsubstituting
for λλ00−−hh
into(11).
(11).
More
s
1 −by
α by
substitutings sfor
forhhandand −h−h00 forMore
speci
0 0intosketch of the proof for DBP

stitutefor γγii.. Then
Then we
as as
α. α.
0 for
i (s,
ute
−λ−λ
λ0λ, 0γ, iγ(s,0)0)forweobtainobtain(27)
(27)
0 for

of Theorem 4. Let s = Lσ2 (h, λ0 ). Then,contour surface of BP
heorem 4. Let s = Lσ2 (h, λ0 ). Then,! τ 2 ,σ2 (H|y) = 1 − BPτ 2 (R(s)|µ̃).
DBP! τ 2 ,σ2 (H|y) = 1 − BPτ 2 (R(s)|µ̃).
DBP

The proof completes by applying the asymptotic expansion of BP to R(s)! τ 2 ,σ2 (H|y) with µ̃ = (0, −h0 ) for θ = 0. Applying The
t compute DBP! τ 2 ,σ2 (H|y) with µ̃ = (0, −h0 ) for θ = 0. Applying T
ompute DBP

2(R(s)|µ̃), we get(s)|µ̃), we get!
"! τ 2 ,σ2 (H|y) ≃ Φ β0′ τ −1 + β1′ τ + β2′ τ 3 ,DBP!
"! τ 2 ,σ2 (H|y) ≃ Φ β0′ τ −1 + β1′ τ + β2′ τ 3 ,
DBP

74
′
′
′
β0 , β1 , β2 are defined by replacing γi by γi (s, 0) of (25) and λ0 by −λ0 ,Theorem 6 (Rejection probability of double bootstrap probabi
en in Theorem 1, and µ = (0, −h0 ) ∈ ∂H, the rejection probabil
棄却確率
ectionHprobabilityof double
probabilities).
given in Theorem
1, bootstrap
and µ = (0,−h0 ) ∈ ∂H,For
the therejection
probabilities of DBP and DAUisσ 2 (H|Y )Rejection
em 1,DBP
and1,σµ= )(0,is−h0 ) ∈ ∂H, the rejection probability of
2 (H|Y!
"#$2
P DBP1,σ2 (H|Y
)
<
α
≃
Φ
z
−
(1
+
σ
α "
!
# )β3 .
$!

2
! P DBP1,σ2 (H|Y" ) < α ≃ mean
Φ curvaturezα − of(1the+σ)β3 .
mean curvature
P DAU(H|Y ) < α ≃ Φ(zα ) = α 平均曲率の平均曲率

DBP is third-order"accurate#(k=3)$! ) < α ≃ Φ"zα − (1 + σ 2 )β3 .
P DBP1,σ2 (H|Y

−3/2

P DBP(H|Y ) < α ≃ Φ(z − 2β ) = α + O(nDiscussion
6 Discussion
α

3

)

(32)

is fourth-orderaccurate (k=4)%
Theorem 6DAU
(Rejection
probabilityof 1double
probabilities). For the∂ 2 h(u)bootstrap
%
!
γ"1 =
%
%
2
2
∂u
∂u
0
1
∂
h(u)
P
DAU(H|Y
)
<
α
≃
Φ(z
)
=
α
i
i
%
H given in Theorem 1, and µ = (0, −hα0 ) ∈ ∂H, the rejection
probability ofon

γ1 = %%
1 ∂ γ1 (h, u) %2 ∂ui ∂ui 0
!
"
β3α =
DBP1,σP2 (H|Y
) is
%−3/2
% 2β
%
DBP(H|Y) < α1 ∂≃2 h(u)
Φ(z
−
)
=
α
+
O(n
3
% 2 ∂ui ∂ui 1 ∂0 2 γ)1 (h, u) %
β3 =
γ1 =
%
%
!
"
#
$
2 ∂ui ∂ui 0
2 ∂ui ∂ui 0
2
2of
P mean
DBP
(H|Y
)
<
α
≃
Φ
z
−
(1
+
σ
)βcurvaturem 6 (Rejection probability
doublebootstrap
probabilities).
α3 . For the1,σmean
curvature
of
the
mean
curvatureofthesurface%2
1 ∂ γ1 (h, u) %
β =
%
2

(32)
75

σ

0

−1
3
2
−1 ! (29)
−1/2
DBP
≃
1
−
Φ
β
τ
−
−
β
τ
−
β
τ
σ
.+
!
−2β1 τ error
0
2
3
"
!
−1
meaning(H|y)
thatLetthecomputationalisO(n) withrespect
tocomputed
the original
parame
DBP
(H|y)
=
DBP
(H|y)
O(n
),
2
2
DBP
(H|y)
denotethedoublebootstrap
probability
by replacing
!
1,1
1,1
DAU
=
PV
+
O(n
)
2
τ
,σ
DBP
(H|y)
=
DBP
(H|y)
+
O(n
1,1λ0 ∈ R, and σ >),0. Then,−132 (h, λ0 ) for h ∈ S,
Lemma 4. Let s = Lσ1,1
s is expressedas 2
τ 2 ,σ 2

2 ,σ 2 (H|y) ≃ 1 − Φ β0 τ
DBPτ!
β1 τ − β2 τdouble− β3 τ σbootst.
ＤＡＵは射影誤差に影響されにくい
say, η ins thespherical
example.
We(H|y)
forO(1),this−modified= M(h,λ) by specifying
λ(u)≃ denote
λ0 + λi uiDBP+ λij τu2i u,σj 2with
λ0 =
−1
q
!
DBP
=
DBP
+
O(n
)
projectionµ̂(H|y)
in
(4)
by
µ̃
=
(θ,
−h(θ))
∈
∂H
for
some
θ
∈
R
. Then,probability,
andderiveits
asymptotic
expansion
for
y
=
(0,
λ
−
h
).
!probability00
2 (H|y) denote the double bootstrap
computedbyreplacing
the
DBP
(H|y)
≃
DBP
(H|y).
!
2
2
1,−1
1,−1
DBP
DBP
(H|y).
λi terms= σ (−3h+ 6λhmli2(H|y)), 2 (H|y)
λij ≃
=
σ
(−6h
+ 2hmmhli hlj + probability
4h h hlj ).computed
(26)
1,−1
mmi
0 hDBP
ml
mmij!
HigherFirst
order
doublebootstrap
probabilities
+1,−1
+
Let
denotethedoublebootstrap
by r
τ
,σ
2
2
2 (h, λ0 ), andnotethat of
BP
(H|Y
)
≥
BP
(H|y)
⇔
Y
∈ R(s)
for s =mlLσmi
σ
σ
! ! q
"
−2
DAU
=
DAU
+
O(n
)
(H|y) in (4)µ̃ =(θ,−h(θ))
∈
∂H
for
some
θ
∈
R
.
Then,
−1
3
2
−1
2
2
Weby
have
sDBP
∈
S
with
coeﬃcients˜ q , (! ,σ2 (H|y) ≃ 1 − Φ β0 τ − β1 τ − β2 τ − β3 τ σ − τ (τ + σ )∆(θ)
projectionµ̂(H|y)
in
(4)
by
µ̃
=
(θ,
−h(θ))
∈
∂H
for
some
θ
∈
R . TheIf µ̂(H|y) =τ 2(0,isreplacedby−h(0))∈∂H!
DBPτ 2 ,σ2 (H|y)2= 1 − BPτ 2 (R(s)|µ̃).
s!0 = h0 − λ0 , si = hi − 2λ0 hm hmi + σ (3hmmi − 6λ0 hml" hmli − 6λ0 hmi hmll ),y22(6h!µ̂(H|y)
−1
2
3 4λ
2 h −12 ˜− 2h µ̃ h h − 4h h h ), (27)!
−1shijmli
=
hi−1
−
2λ
hβ
+
h
h
+
σ
DBP
=
DBP
+
O(n
H|y)
≃
1
−
Φ
β
τ
−
β
−
τ
−
β
τ
σ
−
τ
(τ
+
σ
))−1
∆(θ)
, following
(30)
!
ij
0τhmi
mj
ml
mi
lj
ijmm
mm
li lj 3ml mi2 lj −1
3h
−
6λ
h
)θ
+
(6h
−
2γ
h
h
−
4h
h
h
)θ
θ
0
0
1
2
3
˜
mmi
0
ml
mmij
1
mi
mj
ml
mi
lj
i
j
DBP
=
DBP
+
O(n
˜.
!
By applying where
TheoremBP(R(s)|µ̃),
thetheoremstrai∆(θ)2 =to(3h−6λ0 hml+0 τthe(6h
2 (H|y)
DBP
≃h1we
− )θ
Φget
−
β1 τ−
−2γ
β12 τhmi−hβ
σ 4h−ml
τ hvia(τh2alj+)θ
σi2θ)j∆ττ 22,σ
mmi
mli
iβ
mmij
mj
mi
3 τ−
sijk = hijk − 2λ0 (hmi hmjk + hmj hmik + hmk hmij ), sijkl = hijkl .
forward computation.
−222!
−2)the H
m 4 (Asymptotic
expansionof double
bootstrap
probabilities).
For
particular,
τ
=
1,
σ
=
±1,
gives
!
DAU
=
DAU
+
O(n
The
four
geometric
quantities
of
s
at
(0,
−s(0))
are
DAU
=
DAU
+
O(n
˜
= (3hmmi − 6λ0 hml hmli )θi + (6h
−
2γ
h
h
−
4h
h)mi
hlj)θ
)θii+
θj .(6hInmmij − 2γ1 hmi hmj − 4hml hm
mmij∆(θ) =1 (3hmi mmimj − 6λ0ml
where
h
hmli
ml
Theorem 45.4
(Scaling-law
ofthedoublebootstrap
probability).
For
the H and y = (0, λAsymptotic accuracyofbootstrap
methods2
20, λ0 − h0 ) given in Theorem 1,γthebootstrap(s,double0)≃ γmodified+
4λ0 γ3 +is
σexpressed(6γ4 − 2γprobability11,
1 − 2λ0 γ2probability
1 γ2 − 4γ3 ), with µ̃ = (θ, −h(θ)
inTheoremthedoublebootstrap
2 h0 ) given
2
2
2
(28)
= 1, σ = ±1, gives−1y4λµ̂(H|y)µ̃ ≃particular,τ =1,(H|y)σγ (s,
= 0)
±1,
gives
!
DBP
=
DBP
(H|y)
+
O(n
),
γ
(s,
0)
≃
γ
−
γ
,
γ
,
γ
(s,
0)
≃
γ
,
1,1
1,1
2 as
2
0 3
3
3
4
4
expressed
DBPasymptotically
becomes
The probability
is calculated by (1), and zα = Φ−1 (α), 0 < α < 1.
"
where γi = γi (h, 0), i = 1, . . !. , 4.
−1
−1
!
−1
3
2
−1
2
2
!
DBP
(H|y)
=
DBP
(H|y)
+
O(n
),
!
1,1
DBP
DBP≃
DBP
(H|y)
Φ̄ β0 τ+ O(n− β ),
τ −(H|y)β221,1
τ ≃−DBPβ3 τ σ −τ (τ + σ )κ(θ) ,1,1 (H|y)1,1 (H|y)!
τ 2 ,σ 2 =1,−1 (H|y).
We denotethe λ(u) of (26) as
λσ2DBP(u) =11,−1λ0 −σ κ(u)
with
Theorem 5 accuracy(Rejectionof bootstrap
probabilities). For the H given
5.4 Asymptotic
Asymptotic
accuracyprobabilityof bootstrap
bootstrap
methods
5.4
methods
error = O(n-1)
κ(u)
=
γ
(h,
u)
−
γ
(h,
0)
−
λ
(γ
(h,
u)
−
γ
(h,
0))
! 1,−1 (H|y) ≃ DBP1,−1 (H|y).
1
1
0 2
2 DBP
where κ(θ)!
is definedin (29).
DBP
(H|y)
≃
DBP
(H|y).
1,−11,−12 (H|Y ) is
Theorem
1,20and= (0, −h0 ) ∈ ∂H, therejection
NBPσ(29)≃ (3hmmi− 6λµ2hmm
hli hlj −probability
4hml hmi hlj )uof
0 hml hmli )ui + (6hmmij −−1
i uj .
−1The
probability
is calculated
calculatedby (1),and−hzα0 ).
= By
Φ letting(α),
Theprobabilityby(α), 00 <<α <<1. (31), we obtainWhenhi = 0, iswehave µ̂(H|y)
= (0,
θ α=
0 1.
in
$ DBP
! accurate (k=2),% −1but&
Corollary: DBP#becomes only second-order)! ! = DBP + O(n"4 2
22 (H|Y ) < α accuracy
P keepsNBP≃ Φ −1
zaccuracy+ (1 + σ )ofγ31bootstrap+ z α γ2 +γ3 − γ1 γ2DAU
2 3 zαmethodsα(k=4)
σfourth-order
5.4
Asymptotic
2 ,σ 2 (H|y) ≃ Φ̄ β0 τ
DBP
−
β
τ
−
β
τ
−
β
τ
σ
.
1
2
3
τ
Theorem55 (Rejection(Rejection probability
probability of bootstrap probabilities).
Theoremprobabilities). For
ForthetheHHgiven
giveninin%&"−2! = DBP + O(n−1 )! = DAU + O(n
DBP
DAU
)
2 2
24
4
+with
(1
+(1),
σ )and3γ
−O(n
γ3−1
−for
σ <3which
γα3<. 1.κ(θ) ≃
(
−1
4=
3 2−1
The
probability
iscalculated
by
z
Φ
(α),
0
When
h
=
O(n
),
we
have
µ̂(H|y)=(θ,−h(θ))
some
θ
=
)
α
i
2
Theorem
1,
and
µ
=
(0,
−h
)
∈
∂H,
the
rejection
probability
of
NBP
(H|Y
)
is
Theorem
1,
and
µ
=
(0,
−h
)
∈
∂H,
probability
of
NBP
(H|Y
)
is
00
σσ
76
Efron and Tibshirani (1998)の結果
Shimodaira (2014)
Therefore, (32)!
holds for any h ∈−2
S, and (21) follows. This argument also confirms that

Robustness to projection error

まとめ• ブートストラップ確率はベイズ事後確率• 頻度論のp-値との差は，仮説境界の「平均曲率」が原因• ダブルブートストラップ法でも修正したときの誤差は，「平均曲率の平均曲率」• マルチスケールブートストラップ法で誤差解
消• 検定以外にも使えるかも？ （境界までの距
離や境界の曲率を推定できる）
77

Summary and other issues
•
•
•
•
•

DAU = “DBP with m=-n” is proposed
The accuracy of BP is first order (k=1), AU is third-order (k=3), DBP is third-order(k=3)
DAU is fourth-order accurate (k=4)
DAU is robust to the projection error (surprisingly, k=4)
Geometry of surfaces played important roles•

Shimodaira (2008) showed another theory of AU using unusual asymptotic theoryof “nearly flat surfaces”•

Shimodaira (2004) discussed deviation from the multivariate normal model, andresults for exponential family distributions are given there for multistep-AU•

Future topics may be DAU for nearly flat surfaces, or for exponential family
distributions78