International Symposium on Statistical Theory andMethodology for Large Complex Data
November 16-18, 2018
Venue:
D509 Institute of Natural Sciences,

University of Tsukuba

1-1-1 Tennodai, Tsukuba, Ibaraki 305-8571, Japan
Organizers:
Makoto Aoshima (University of Tsukuba)
Mika Sato-Ilic(University of Tsukuba)

Kazuyoshi Yata

(University of Tsukuba)

Aki Ishii(Tokyo University of Science)

Supported by
Grant-in-Aid for Scientific Research (A) 15H01678 (Project Period: 2015-2019)“Theories and methodologies for large complex data”(Principal Investigator: Makoto Aoshima)
Grant-in-Aid for Challenging Research (Exploratory) 17K19956 (Project Period: 2017-2018)“New developments for big data by non-sparse modeling”(Principal Investigator: Makoto Aoshima)
Grant-in-Aid for Scientific Research (S) 18H05290 (Project Period: 2018-2022)“Introduction of general causality to various observations and the innovation for its optimal statistical inference”(Principal Investigator: Masanobu Taniguchi)

Program
November 16 (Friday)
14:00∼14:05

Opening

14:05∼14:45 Aki Ishii∗,a , Kazuyoshi Yatab and Makoto Aoshimabab(Department of Information Sciences, Tokyo University of Science)(Institute of Mathematics, University of Tsukuba)

Tests of high-dimensional mean vectors under the SSE model
14:55∼15:35 Hiroumi Misaki(Faculty of Engineering, Information and Systems, University of Tsukuba)

Comparison of ﬁnancial volatility estimators: RK, TS, PA and SIML
15:50∼16:30 Shota Katayama(Department of Industrial Engineering and Economics, Tokyo Institute of Technology)

Robust and sparse Gaussian graphical modelling under cell-wise contamination

16:40∼17:20 Mariko Yamamura(Graduate School of Education, Hiroshima University)

Estimation of spatiotemporal eﬀects by the fused lasso for densely sampled
spatial data using body condition data set from common minke whales

November 17 (Saturday)
9:25∼10:05 Kazuyoshi Yata∗ and Makoto Aoshima(Institute of Mathematics, University of Tsukuba)

A high-dimensional quadratic classiﬁer after feature selection
10:15∼10:55 Junichi Hirukawaa,∗ and Sangyeol Leebab

(Faculty of Science, Niigata University)(Department of Statistics, Seoul National University)

Asymptotic properties of mildly explosive processes
with locally stationary disturbance
11:05∼11:45 Hirokazu Yanagihara(Department of Mathematics, Graduate School of Science, Hiroshima University)

High-dimensionality adjusted asymptotically loss eﬃcient GCp criterion
in normal multivariate linear regression models
11:45∼13:15

Lunch

13:15∼17:30

Special Invited and Keynote Sessions

18:30∼

DinnerNovember 18 (Sunday)
9:40∼10:20 Kei Hirosea,∗ and Hiroki Masudabab(Institute of Mathematics for Industry, Kyushu University)(Faculty of Mathematics, Kyushu University)

Robust relative error estimation
10:30∼11:10 Tomonari Sei(Graduate School of Information Science and Technology, University of Tokyo)

Inconsistency of diagonal scaling under high-dimensional limit
11:20∼12:00 Takafumi Kanamori(Department of Mathematical and Computing Sciences, Tokyo Institute of Technology)

Statistical inference with unnormalized models
12:00∼ 12:10

Closing(∗ Speaker)

Special Invited Session
13:15∼14:05

Fixed support positive-deﬁnite modiﬁcation of covariance matrix estimators
via linear shrinkage

Speaker: Johan Lim(Department of Statistics, Seoul National University)

14:20∼15:10

Greedy active learning algorithm for logistic regression models

Speaker: Ray-Bing Chen
(Department of Statistics, National Cheng Kung University)

15:25∼16:15

Inference for LSHD time series models and applications to sensor monitoring
and ﬁnancial engineering

Speaker: Ansgar Steland(Institut fuer Statistik und Wirtschaftsmathematik, RWTH Aachen University)

Keynote Session
16:30∼17:30

Do we necessarily have more cost-eﬀective information with big data?

Speaker: Nitis Mukhopadhyay(Department of Statistics, University of Connecticut-Storrs)

Discussion Leader: Chikara Uno(Faculty of Education and Human Studies, Akita University)

Tests of high-dimensional mean vectors under the SSE model
Aki Ishiia , Kazuyoshi Yatab and Makoto Aoshimaba

Department of Information Sciences, Tokyo University of SciencebInstitute of Mathematics, University of Tsukuba

1 Introduction
We considered statistical inference on mean vectors in the high-dimension, low-sample-size (HDLSS) context. Let x1 , . . . , xn be a random sample of size n(≥ 4) from a p-variate distribution with an unknown
mean vector µ and unknown covariance matrix Σ. In the HDLSS context, the data dimension p is very
high and n is much smaller than p. We define the eigen-decomposition of Σ by Σ = HΛH T , whereΛ = diag(λ1 , ..., λp ) is a diagonal matrix of eigenvalues, λ1 ≥ · · · ≥ λp ≥ 0, and H = [h1 , ..., hp ] isan orthogonal matrix of theeigenvectors.
We write the sample mean vector and the sample∑n∑correspondingncovariance matrix as x = j=1 xj /n and S = j=1 (xj − x)(xj − x)T /(n − 1).
In this talk, we discussed the one-sample test:
H0 : µ = µ0

vs.

H1 : µ ̸= µ0 ,(1.1)

where µ0 is a candidate mean vector. We assume µ0 = 0 without loss of generality. One should note that
Hotelling’s T 2 -statistic is not available because S −1 does not exist in the HDLSS context. [6, 7] considered
the test when X is Gaussian. When X is non-Gaussian, [5] considered the test. Let us consider the following
eigenvalue condition:
λ21→ 0 as p → ∞.
(1.2)
tr(Σ2 )
Under (1.2), H0 and some regularity conditions, [5] and [1, 2] showed the asymptotic normality for their test
statistics.
[3] called the eigenvalue condition (1.2) the “non-strongly spiked eigenvalue (NSSE) model” and drewattention that high-dimensional data do not fit the NSSE model on several occasions. In order to overcome
this inconvenience, [3] proposed the “strongly spiked eigenvalue (SSE) model” defined by{ λ2 }1>0p→∞tr(Σ2 )

lim inf

and gave a data transformation technique from the SSE model to the NSSE model.

2 A new test procedure under the SSE modelLet Ψr = tr(Σ2 ) −∑r−1

2
s=1 λs =∑p2
s=r λsfor r = 1, ..., p. We assumed the following model:(A-ii) There exists a fixed integer k (≥ 1) such that(i) When k ≥ 2, λ1 , ..., λk are distinct in the sense that
lim inf (λr /λs − 1) > 0 for 1 ≤ r < s ≤ k;p→∞(ii) λk and λk+1 satisfy

λ2k+1
λ2k
lim inf> 0 and→ 0 as p → ∞.
p→∞ ΨkΨk+1(1.3)

Note that (A-ii) is one of the SSE models. According to [3, 4], we considered transforming the data from the
SSE model to the NSSE model by using the projection matrix
A = Ip −k∑hj hTj =j=1

p∑

hj hTj .

j=k+1

We have that E(Axj ) = Aµ (= µ∗ , say) andVar(Axj ) = AΣA =p∑

λj hj hTj (= Σ∗ , say).

j=k+1

Note that tr(Σ2∗ ) = Ψk+1 and λmax (Σ∗ ) = λk+1 , where λmax (Σ∗ ) denotes the largest eigenvalue of Σ∗ .
Then, it holds that
λ2max (Σ∗ )/tr(Σ2∗ ) → 0 as p → ∞ under (A-ii).
Thus, the transformed data has the NSSE model.
By using the transformed data, we consider the following quantity:)∑n ( T∑k∑nTtr(AS)
l<l′ xl xl′ −j=1 xjl xjl′
l<l′ xl Axl′
2
TDT = ∥Ax∥ −
=2
=2
,
n
n(n − 1)n(n − 1)where
xjl = hTj xl for all j, l.
We discussed the asymptotic null distribution and the power of TDT . We applied the findings to the construction of confidence regions on the mean vector under the SSE model. We further discussed multi-sample
problems under the SSE models. Finally, we demonstrated the new test procedure by using actual microarraydata sets.

References[1] Aoshima, M., Yata, K. (2011). Two-stage procedures for high-dimensional data. Sequential Analysis(Editor’s special invited paper), 30, 356–399.
[2] Aoshima, M., Yata, K. (2015). Asymptotic normality for inference on multisample, high-dimensionalmean vectors under mild conditions. Methodology and Computing in Applied Probability 17, 419–439.
[3] Aoshima, M., Yata, K. (2018a). Two-sample tests for high-dimension, strongly spiked eigenvalue models. Statistica Sinica 28, 43–62.
[4] Aoshima, M., Yata, K. (2018b). Distance-based classifier by data transformation for high-dimension,strongly spiked eigenvalue models. Annals of the Institute of Statistical Mathematics, in press(doi:10.1007/s10463-018-0655-z).
[5] Bai, Z., Saranadasa, H. (1996). Effect of high dimension: By an example of a two sample problem.
Statistica Sinica, 6, 311–329.
[6] Dempster, A.P. (1958). A high dimensional two sample significance test. The Annals of Mathematical
Statistics, 29, 995–1010.
[7] Dempster, A.P. (1960). A significance test for the separation of two highly multivariate small samples.
Biometrics, 16, 41–50.

Comparison of financial volatility estimators:
RK, TS, PA and SIML
Hiroumi Misaki 1
In the recent decade, several estimation methods of volatility of financial asset prices
with high-frequency data have been developed to deal with the problem of market microstructure noise (MMN). Bandorﬀ-Nielsen et al. (2008) proposed the realized kernel(RK) method using autocovariances and kernels, which is widely used for empirical analysis in academic literatures. The two-scale (TS) estimator by Zhang et al. (2005) combines
sparsely and densely (sub)sampled realized volatility to eliminate the bias. Jacod et al.
(2009) developed pre-averaging (PA) approach for the estimation problem with MMN.
Kunitomo and Sato (2011, 2013) have proposed the separating information maximum
likelihood (SIML) method. Misaki and Kunitomo (2015) and Kunitomo, Misaki and Sato(2015) have further investigated the properties of the SIML estimation with the data thatis randomly sampled and includes MMN. The analytical investigations and the batch ofsimulations have shown that the SIML estimator has reasonable asymptotic properties aswell as finite sample properties.
Although some alternative methods have been proposed as described above, the relative
properties of these methods in actual market structures are yet to be assessed. Therefore,
in this presentation, we make a comparative study on RK, TS, PA and SIML.
We assume that the underlying continuous process X(t) (0 ≤ t ≤ 1) is that∫ t

X(t) = X(0) +0

σx (s)dB(s) (0 ≤ t ≤ 1),

where B(s) represents standard Brownian motion and σx (s) is the instantaneous volatility
function. The main statistical goal is to estimate the integrated volatility
σx2 =

∫ 1
0

σx2 (s)ds

of the underlying continuous process X(t) from the set of discretely observed prices y(tni )()

that are generated by y(tni ) = h X(tni ), y(tni−1 ), u(tni ) , where u(tni ) is market MMN.
We investigate the finite sample properties of the five estimators for the integrated
volatility based on a set of simulations. Generally, we can assume deterministic volatility][

function σx2 (s) = σ(0)2 a0 + a1 s + a2 s2 where ai (i = 0, 1, 2) are constants and σx (s)2 > 0
for s ∈ [0, 1]. In this case the integrated volatility is given by
σx2 =
1

∫ 1
0

[

σx (s)2 ds = σx (0)2 a0 +

a1 a2
+
2
3

]

.

Faculty of Engineering, Information and Systems, University of Tsukuba, Tennodai 1-1-1, Tsukuba
City, Ibaraki 305-8577, JAPAN, hmisaki@risk.tsukuba.ac.jp

We use several non-linear transformation models for the form of MMN in addition to
the standard ‘addtive plus noise’ model: h(x, y, u) = x + u. Each model corresponds toModel 1h(x, y, u) = gη (x + u),

Model 2h(x, y, u) = y + b(x − y),Model 3h(x, y, u) = y + gη (x − y + u),

Model 4h(x, y, u) = y + gη (x − y) + u,

Model 5 b1 (x − y) if x ≥ yh(x, y, u) = y + u + b (x − y) if y < x,2where gη (x) = η [x/η], and b, b1 , b2 are constants, respectively.
In the Tables we show the some of simulation results except for PA. Figures of MSE
for each estimator are attached. As we can see from the Figures, the RK is severely biased
to the round-oﬀ errors when the noise is small, whereas the TS is sensitive to the linear
price adjustment when the noise is relatively large. The SIML and LSIML are not thebest in most of the cases including the standard one, but they seem to have robustness tothe form of noise, compared to the other estimators. We have shown the more findings at
the session.

References[1] Barndorﬀ-Nielsen, O., P. Hansen, A. Lunde and N. Shephard (2008), “Designing realized kernels to measure the ex-post variation of equity prices in the presence of noise,” Econometrica,
Vol.76-6, 1481-1536.
[2] Jacod, J., Y. Li, P. A. Mykland, M. Podolskij and M. Vetter (2009), “Microstructure noise in
the continuous case: The pre-averagin approach,” Stochastic Processes and their Applications,Vol. 119, 2249-2276.
[3] Kunitomo, N., H. Misaki and S. Sato (2015), “The SIML Estimation of Integrated Covarianceand Hedging Coeﬃcients with Micro-market noises and Random Sampling,” Asia-PacificFinancial Markets, Vol. 22, 3, 333-368.
[4] Kunitomo, N. and S. Sato (2011), “The SIML Estimation the Integrated Volatility of Nikkei225 Futures and Hedging Coeﬃcients with Micro-Market Noise,” Mathematics and Computersin Simulations, Elsevier, 81, 1272-1289.
[5] Kunitomo, N. and S. Sato (2013), “Separating Information Maximum Likelihood Estimation of the Integrated Volatility and Covariance with Micro-Market Noise,” North American
Journal of Economics and Finance, Vol. 26, 282-309.
[6] Misaki, H. and N. Kunitomo (2015), “On robust properties of the SIML estimation of volatility under micro-market noise and random sampling,” International Review of Economics &Finance, 40, 265-281.
[7] Zhang, L. , P. Mykland and Ait-Sahalia, Y., (2005), “A tale for two time scale : determining
integrated volatility with noisy-high frequency data,” Journal of the American Statistical
Association, Vol. 100(472), 1394-1411.

Robust and sparse Gaussian graphical modelling undercell-wise contamination
Shota Katayama1 , Hironori Fujisawa2 and Mathias Drton3
1
2

Tokyo Institute of Technology, JapanThe Institute of Statistical Mathematics, Japan
3
University of Washington, USA

Let Y = (Y1 , . . . , Yp )T be a p-dimensional random vector representing a multivariate
observation. The conditional independence graph of Y is the undirected graph G = (V, E)
whose vertex set V = {1, . . . , p} indexes the individual variables and whose edge set E
indicates conditional dependences among them. More precisely, (i, j) ̸∈ E if and only if Yi
and Yj are conditionally independent given YV \{i,j} = {Yk : k ̸= i, j}. For a Gaussian vector,
the edge set E corresponds to the support of the precision matrix. Indeed, it is well known
that if Y follows a multivariate Gaussian distribution Np (µ, Σ) with mean vector µ andcovariance matrix Σ, then (i, j) ̸∈ E if and only if Ωij = 0, where Ω = Σ−1 .
Inference of the conditional independence graph sheds light on direct as opposed to indirect interactions and has received much recent attention (Drton and Maathuis, 2017). In
particular, for high-dimensional Gaussian problems, several techniques have been developedthat exploit available sparsity in inference of the support of the precision matrix Ω. Meinshausen and Bühlmann (2006) suggested fitting node-wise linear regression models with ℓ1
penalty to recover the support of each row. Yuan and Lin (2007), Benerjee et al. (2008) andFriedman et al. (2008) considered the graphical lasso (Glasso) that involves the ℓ1 penalized log-likelihood function. Cai et al. (2011) proposed the constrained ℓ1 minimization for
inverse matrix estimation (CLIME), which may be formulated as a linear program.
In fields such as bioinformatics and economics, data are often not only high-dimensional
but also subject to contamination. While suitable for high dimensionality, the above mentioned techniques are sensitive to contamination. Moreover, traditional robust methods may
not be appropriate when the number of variables is large. Indeed, they are based on themodel in which an observation vector is either without contamination or fully contaminated.
Hence, an observation vector is treated as an outlier even if only one of many variables

is contaminated. As a result these methods down-weight the entire vector regardless ofwhether it contains ‘clean’ values for some variables. Such information loss can become fatal
as the dimension increases. As a more realistic model in high dimensional data, Alqallaf etal. (2002) considered cell-wise contamination: the observations X1 , . . . , Xn with p variables
are generated by
Xi = (Ip − Ei )Yi + Ei Zi ,i = 1, . . . , n.

(1)Here, Ip is the p × p identity matrix and each Ei = diag(Ei1 , . . . , Eip ) is a diagonal randommatrix with the Eij ’s independent and Bernoulli distributed with P (Eij = 1) = εj . The
random vectors Yi and Zi are independent, and Yi ∼ Np (µ, Σ) corresponds to a clean
sample while Zi makes contaminations in some elements of Xi .
Our goal is to develop a robust estimation method for the conditional independence
graph G of Yi from the cell-wise contaminated observations Xi . Techniques such as nodewise regression, Glasso and CLIME process an estimate of the covariance matrix. Our
strategy is thus simply to apply these procedures using a covariance matrix estimator thatis robust against cell-wise contamination. However, while many researchers have considered
the traditional ‘whole-vector’ contamination framework (see, e.g., Maronna et al., 2006),
there are fewer existing methods for cell-wise contamination. Specifically, we are aware ofthree approaches, namely, use of alternative t-distributions Finegold and Drton (2011), useof rank correlations (Loh and Tan, 2015; Öllerer and Croux, 2015), and a pairwise covarianceestimation method by Tarr et al. (2016) who adopt an idea of Gnanadesikan and Kettenring(1972). In contrast, in this talk, we provide a robust covariance matrix estimator via γdivergence as proposed by Fujisawa and Eguchi (2008). The γ-divergence can automaticallyreduce the impact of contaminations, and it is known to be robust even when the numberof contaminations is large.

Estimation of spatiotemporal eﬀects by the fused lasso for densely sampled spatial data using body condition data set from common minke whales
Mariko Yamamura1 , Hirokazu Yanagihara2 , Keisuke Fukui3 ,Hiroko Solvang4 , Nils Øien4 , Tore Haug4
1

Graduate School of Education, Hiroshima University, 2 Graduate School of Science, Hiroshima
University, 3 Research & Development Center, Osaka Medical College,
4

Institute of Marine Research, Norway,

Samples evenly distributed all over the population are not always available for real data
analysis. As an example of a spatial data, a data from common mink whales in Norwegianwater provides values showing their body conditions with whaling locations such as longitudes
and latitudes. Though whales are distributed all over the Norwegian water, whaling locations
are almost the same every year, therefore samples are dense at particular locations. The space
to be analyzed is subdivided into several, and we estimate the spatial eﬀect by using fused lasso
with combining spatial eﬀects from subdivided space.
Let y and ε be n-dimensional vectors obtained by stacking the vectors of response variables′ )′ and ε = (ε′ , . . . , ε′ )′ ,and error variables of the j-th space, respectively, i.e., y = (y1′ , . . . , ymm1

and let X and B be n × k and n × b matrices obtained by stacking the matrices of explanatory′ )′ and B =variables and basis functions of the j-th space, respectively, i.e., X = (X1′ , . . . , Xm′ )′ . As a whole of space, the additive model is written such as(B1′ , . . . , Bm

y = Xβ + Bα + Rµ + ε,
where µ = (µ1 , . . . , µm )′ and R is an n × m matrix defined by
1n1 ⊗ e′1


..


R=
.
.


′
1nm ⊗ emHere, ej is the m-dimensional vector of which the j-th element is 1 while all the other elements
are 0, and ⊗ indicates the Kronecker product of the two matrices.
Yanagihara (2012) shows that choosing the smoothing parameters in the penalized smoothing
spline is equivalent to choosing the ridge parameters in the generalized ridge regression using
the matrix of transformed basis function values as the matrix of explanatory variables. And
then Yanagihara (2018) considers optimization of the ridge parameters in generalized ridge1regression by minimizing a model selection criterion, i.e., generalized cross-validation (GCV).
From Yanagihara (2012, 2018), estimates of α and β after optimizing smoothing parameters by
GCV are given byα̂γ = Qγ Vγ Dγ−1/2 zγ ,β̂γ = (X ′ X)−1 X ′ (y − Rµ̂ − B α̂γ ),respectively, where zγ = (z1 , . . . , zγ )′ , and Vγ is a γ × γ diagonal matrix as(

Vγ = diag(νγ,1 , . . . , νγ,γ ),νγ,j = I zj2 > s2γ,a∗γ

)

(
1−

s2γ,a∗γ)zj2(j = 1, . . . , γ),

where the γ is optimized by the GCV.
The penalized residual sum of squares (PRSSλ ) for the adaptive fused lasso is given by
PRSSλ (µ|fˆ) = ∥y − X β̂ − B α̂ − Rµ∥2 + λm ∑∑wjℓ |µj − µℓ |,j=1 ℓ∈Dj

where λ is the non-negative regularization parameter. The spatial eﬀect µ is estimated byminimizing PRSSλ asµ̂λ = arg minm PRSSλ (µ|fˆ).
µ∈R

The above minimization problem can be solved by the coordinate descent algorithm in Friedman
et al. (2007).
There is the geographic distribution of the five International Whaling Commission (IWC)
management areas. We subdivide each areas to have about 300 samples, and estimate µ̂ ofsubdivided areas. If the µ̂ of a subdivided area is equal to the one of its neighbor area, thesubdivided area and the neighbor are united. From the fused lasso estimation result, the subdivided areas are narrowed down the 11 spaces, and whales have the thickest blubber in thenorthernmost space.

References[1] Friedman, J., Hastie, T., Höfling, H. and Tibshirani, R. (2007). Pathwise coordinate optimization. Annals ofApplied Statistics, 1, 302–332.
[2] Yanagihara, H. (2012). A non-iterative optimization method for smoothness in penalized spline regression.
Statistics and Computing, 22, 527–544.
[3] Yanagihara, H. (2018). Explicit solution to the minimization problem of generalized cross-validation criterion
for selecting ridge parameters in generalized ridge regression. Hiroshima Mathematical Journal, 48, 203–222.

2

A high-dimensional quadratic classifier after feature selection
Kazuyoshi Yata and Makoto Aoshima
Institute of Mathematics, University of Tsukuba, Ibaraki, Japan

1 Introduction
A common feature of high-dimensional data is that the data dimension is high, however, thesample size is relatively low. This is the so-called “HDLSS” or “large p, small n” data situation,
here p is the data dimension and n is the sample size. In this paper, we mainly focus on the casewhen “n/p → 0”. Suppose we have independent and p-variate two populations, πi , i = 1, 2,
having an unknown mean vector µi = (µi1 , ..., µip )T and unknown positive-definite covariancematrix Σi for each i. Letµ12 = µ1 − µ2 = (µ121 , ..., µ12p )T .
We assume lim supp→∞ |µ12j | < ∞ for all j. Note that lim supp→∞ ||µ12 ||2 /p < ∞, where || · ||denotes the Euclidean norm. Let σi(j) be the j-th diagonal element of Σi for j = 1, ..., p (i =1, 2). We assume that σi(j) ∈ (0, ∞) as p → ∞ for all i, j. For a function, f (·), “f (p) ∈(0, ∞) as p → ∞” implies that lim inf p→∞ f (p) > 0 and lim supp→∞ f (p) < ∞. Here,“lim inf p→∞ f (p)” and “lim supp→∞ f (p)” are the limit inferior and the limit superior of f (p),respectively. Then, it holds that tr(Σi )/p ∈ (0, ∞) as p → ∞ for i = 1, 2. We do not assume Σ1 = Σ2 . The eigen-decomposition of Σi is given by Σi = H i Λi H Ti , where Λi =diag(λi1 , ..., λip ) is a diagonal matrix of eigenvalues, λi1 ≥ · · · ≥ λip > 0, and H i = [hi1 , ..., hip ]is an orthogonal matrix of the corresponding eigenvectors. We have independent and identically
distributed (i.i.d.) observations, xi1 , ..., xini , from each πi , where xik = (xi1k , ..., xipk )T , k =1, ..., ni . We assume ni ≥ 2, i = 1, 2. We estimate µi and Σi byT

xini = (xi1ni , ..., xipni ) =ni∑xik /nik=1∑ni

and S ini = k=1 (xik − xini )(xik − xini )T /(ni − 1). Let sini (j) be the j-th diagonal element of
S ini for j = 1, ..., p (i = 1, 2). Let x0 = (x01 , ..., x0p )T be an observation vector of an individual
belonging to one of the two populations. We assume x0 and xij s are independent. Let
nmin = min{n1 , n2 } and m = min{p, nmin }.
Note that the divergence condition “p → ∞, n1 → ∞ and n2 → ∞” is equivalent to “m → ∞.
Let |M | be the determinant of a square matrix M . When πi s are Gaussian, the Bayes
optimal rule (the minimum-error rate discriminant function) is given as follows: One classifies
the individual into π1 if−1T −1(x0 − µ1 )T Σ−1
1 (x0 − µ1 ) − log |Σ2 Σ1 | < (x0 − µ2 ) Σ2 (x0 − µ2 )

1

(1.1)

and into π2 otherwise. Since µi s and Σi s are unknown, one usually considers the following
typical classifier:−1
T −1(x0 − x1n1 )T S −1
1n1 (x0 − x1n1 ) − log |S 2n2 S 1n1 | < (x0 − x2n2 ) S 2n2 (x0 − x2n2 ).

The classifier usually converges to the Bayes optimal classifier when nmin → ∞ while p isfixed or nmin /p → ∞. However, in the HDLSS context, the inverse matrix of S ini does notexist. When Σ1 = Σ2 , Bickel and Levina (2004) considered an inverse matrix defined byonly diagonal elements of the pooled sample covariance matrix. When Σ1 ̸= Σ2 , Dudoit et al.
(2002) considered an inverse matrix defined by only diagonal elements of S ini . Aoshima andYata (2011) considered using {tr(S ini )/p}−1 I p instead of S −1
ini from a geometrical backgroundof HDLSS data and proposed geometric classifiers. Here, I p denotes the identity matrix ofdimension p. Chan and Hall (2009) and Aoshima and Yata (2014, 2018) considered distancebased classifiers and Aoshima and Yata (2014) gave the misclassification rate adjusted classifier
for multiclass, high-dimensional data whose misclassification rates are no more than specified
thresholds.
In this talk, we considered classifiers by the diagonal elements of S ini . We provided a
DQDA type classifier by feature selection and show that it has the consistency property evenwhen nmin /p → 0.

References[1] Aoshima M, Yata K (2011) Two-stage procedures for high-dimensional data. Seq Anal
(Editor’s special invited paper) 30:356–399[2] Aoshima M, Yata K (2014) A distance-based, misclassification rate adjusted classifier for
multiclass, high-dimensional data. Ann I Stat Math 66:983–1010[3] Aoshima M, Yata K (2018) Distance-based classifier by data transformation for
high-dimension, strongly spiked eigenvalue models. Ann I Stat Math, in press(doi:10.1007/s10463-018-0655-z)[4] Bickel PJ, Levina E (2004) Some theory for Fisher’s linear discriminant function, ‘naive
Bayes’, and some alternatives when there are many more variables than observations.
Bernoulli 10:989–1010[5] Chan YB, Hall P (2009) Scale adjustments for classifiers in high-dimensional, low sample
size settings. Biometrika 96:469–478[6] Dudoit S, Fridlyand J, Speed TP (2002) Comparison of discrimination methods for theclassification of tumors using gene expression data. J Am Stat Assoc 97:77–872                                         







 

 !    " #$%&'    #$%$'( ) 

#'

¾

    +          ,    

½+      

 





½





½ '#

½

¾

*  ½½' -#   '    





#

 

      !( ./         ,                          #'-    #$%$'

            ,                ,(  !( 0

Æ     

1 

 

   ,     

 #2334 '          # ½ '( 5      *     

#'( "   3      !                                (             0

1  #2334 '                    +   

#2'

2



.      0 1    #  '  

        ¼

 #2334 '      ¼

# '( . 

   (  




6  #2334'

                 ,         



 

( 1 





0 #233&'     

       





1 

 #232'         7   ( 8 ( #234'            ,       ,                  (              ( .                             (                     0    #0 #$9%''           7     
$$4'       

   

 :   #$$9  $$9  $$9     (                    ( :   #232'      
( 5          ( .       







     ;       ;        

 +   " #$&3'    " #$$$'(

            

 #232'( <  





                 ;                 #233$'  0  ( #23'( <         

 0 =         : 7< 



7     ! 

     (         5            













     



 ¼ 










# '







¼

 





 

½































 - 























 

 #' 

 # '  
¼½
 ¼












   # '
¼½
 ¼







¼





¼

½




½

      7    >3?
  >3 ?
   1 Æ         
¼

   3 ¾ 

  3



 



        1 #'  - 

 

½





½









    7

-







,      (Æ-"      *     
#



 -

#@'



   











     





½

¾


½
'

½


#



'

½


 


 

 ½



¾

  





#A'





-




 



½





        



2

 

½

¾

-







½

 





 ´½ ½µ #3' #'  ´¼ ½µ  







     

















     


 






½







   ( "        (        ,    # ;'   B#' 

. 

 ½  ½ ¾
 ½  ½½

 ½ 

 


 

 

 

2









Æ-High-dimensionality Adjusted Asymptotically Loss Efficient GCp Criterion
in Normal Multivariate Linear Regression Models
Hirokazu Yanagihara
Department of Mathematics, Graduate School of Science, Hiroshima University
1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8626, JapanThe multivariate linear regression model is one of basic models of multivariate analysis.
This model is introduced in many multivariate statistical textbooks (see e.g., Srivastava, 2002,
chap. 9; Timm, 2002, chap. 4), and even now is widely used in chemometrics, engineering,
econometrics, psychometrics, and many other fields, for the predication of multiple responses toa set of explanatory variables. Let Y = (y1 , . . . .yn )′ be an n × p matrix of p response variables,
and let X = (x1 , . . . , xn )′ be an n × k matrix of non-stochastic k explanatory variables, where nis the sample size. In order to ensure the possibility of estimating the model and the existenceof a variable selection criterion, we assume that rank(X) = k (< n) and n − p − k − 1 > 0.
Suppose that j denotes a subset of ω = {1, . . . , k} containing kj elements, and Xj denotes then × kj matrix consisting of the columns of X indexed by the elements of j, where kA denotes thenumber of elements of a set A, i.e., kA = #(A). For example, if j = {1, 2, 4}, then Xj consists ofthe first, second, and fourth columns of X. We then consider the following multivariate linear
regression model with kj explanatory variables as the candidate model:
Y ∼ Nn×p (Xj Θj , Σj ⊗ In ),(1)

where Θj is a kj × p unknown matrix of regression coefficients, and Σj is a p × p unknown
covariance matrix. We identify the candidate model by the the set j and call the candidate
model in (1) the model j. Especially, the model with Xω (namely X) is called the full model.
We will assume that the data are generated from the following true model:
Y ∼ Nn×p (Xj∗ Θ∗ , Σ∗ ⊗ In ),(2)

where j∗ is a set of integers indicating the subset of explanatory variables in the true model,Θ∗ is a kj∗ × p matrix of true regression coefficients, and Σ∗ is a p × p true covariance matrix
with rank(Σ∗ ) = p. We call the model in (2) the true model j∗ . Henceforth, for simplicity, we
represent Xj∗ and kj∗ as X∗ and k∗ , respectively.
In this paper, we focus on a variable selection method by minimizing a generalized Cp (GCp )criterion, which is called a GCp -minimization method, when p may be large but still smallerthan n. The GCp criterion is defined by adding a positive constant value α times the number ofparameters in the mean structure to the minimum value of the residual sum of squares (RSS).
The GCp criterion in an univariate linear regression model was proposed by Atkinson (1980)
and that in a multivariate linear regression model was proposed by Nagai et al. (2012). The
family of GCp criteria contains many widely known variable selection criteria, e.g., Mallows
Cp criterion proposed by Sparks et al. (1983) (the original Cp was proposed by Mallows, 1973,
under the univariate linear regression model), the modified Cp (M Cp ) criterion proposed by
Fujikoshi and Satoh (1997), which is the completely bias-corrected version of Cp criterion. Since
we deal with the multivariate linear regression model of which the dimension p may be large,
the following asymptotic framework is used for assessing an asymptotic property of a variableselection method:• Large-sample and high-dimensional (LSHD) asymptotic framework: n and p/n approach ∞and c0 ∈ [0, 1), respectively. For simplicity, we will write it “n → ∞, p/n → c0 ”.
It should be emphasized that we do not care about whether p goes to ∞ or not in the above
asymptotic framework. Suppose that regression coefficients of the ath explanatory variable(a ∈ j∗ ) and diagonal elements of the true covariance matrix are corresponding to sequences
∗ }
∗
{βa,i
i=1,2,... and {ψi }i=1,2,... , respectively. These mean that new elements are added to current

regression coefficients and diagonal elements of the true covariance matrix when p increases.
There are two important properties of a variable selection method by minimizing a variableselection criterion. One is a consistency property that a selection probability of the true modelby a variable selection criterion goes to 1 asymptotically, the other is an efficiency property
that a ratio of a loss function of a selected model by a variable selection criterion and theminimum loss function goes to 1 asymptotically (see e.g., Shibata, 1980, 1981; Shao, 1997).
Recently, Yanagihara (2016) clarified a sufficient condition of α of GCp criterion in (1) to satisfy
a consistency property under the LSHD asymptotic framework. An aim of this paper is to derivea sufficient condition of α of GCp criterion in (1) to satisfy an efficiency property under the LSHD
asymptotic framework. Then, we can propose asymptotically loss efficient GCp criterion evenunder high-dimensionality of p. The asymptotically loss efficient GCp criterion is defined by thefollowing α:α(β) =

2n
+ β,
n−p

β > 0 s.t.

limn→∞,p/n→c0√

pβ = ∞

andp
β = 0.
n→∞,p/n→c0 n
lim

(3)

Fixed support positive-definite modification ofcovariance matrix estimators via linear shrinkage
Young-Geun Choi
Data R&D Center, SK Telecom, Seoul, Korea
Johan Lim
Department of Statistics, Seoul National University, Seoul, Korea
Anindya Roy and Junyong Park
Department of Mathematics and Statistics, University of Maryland,
Baltimore County, MD, USA

Abstract
In this work, we study the positive definiteness (PDness) problem in covariancematrix estimation. For high dimensional data, many regularized estimators areproposed under structural assumptions on the true covariance matrix including
sparsity. They are shown to be asymptotically consistent and rate-optimal in estimating the true covariance matrix and its structure. However, many of them do nottake into account the PDness of the estimator and produce a non-PD estimate. To
achieve the PDness, researchers consider additional regularizations (or constraints)
on eigenvalues, which make both the asymptotic analysis and computation muchharder. In this paper, we propose a simple modification of the regularized covariancematrix estimator to make it PD while preserving the support. We revisit the idea of
linear shrinkage and propose to take a convex combination between the first-stage
estimator (the regularized covariance matrix without PDness) and a given formof diagonal matrix. The proposed modification, which we denote as FSPD (Fixed
Support and Positive Definiteness) estimator, is shown to preserve the asymptotic
properties of the first-stage estimator, if the shrinkage parameters are carefully selected. It has a closed form expression and its computation is optimization-free,unlike existing PD sparse estimators. In addition, the FSPD is generic in the sensethat it can be applied to any non-PD matrix including the precision matrix. The
FSPD estimator is numerically compared with other sparse PD estimators to understand its finite sample properties as well as its computational gain. It is also
applied to two multivariate procedures relying on the covariance matrix estimator1

— the linear minimax classification problem and the Markowitz portfolio optimization problem — and is shown to substantially improve the performance of bothprocedures.
Key words: Covariance matrix; fixed support; high dimensional estimation; linear
minimax classification problem; linear shrinkage; mean-variance portfolio optimization; precision matrix; positive definiteness.

2

Greedy Active Learning Algorithm for Logistic Regression Models

Ray‐Bing Chen
Department of Statistics,
National Cheng Kung University

We study a logistic model‐based active learning procedure for binary classification
problems, in which we adopt a batch subject selection strategy with a modified
sequential experimental design method. Moreover, accompanying the proposed
subject selection scheme, we simultaneously conduct a greedy variable selection
procedure such that we can update the classification model with all labeled training
subjects. The proposed algorithm repeatedly performs both subject and variableselection steps until a prefixed stopping criterion is reached. Our numerical results
show that the proposed procedure has competitive performance, with smaller
training size and a more compact model compared with that of the classifier trainedwith all variables and a full data set. We also apply the proposed procedure to awell‐known wave data set (Breiman et al., 1984) and a MAGIC gamma telescope data
set to confirm the performance of our method.

Inference for LSHD Time Series Models and Applications to Sensor Monitoringand Financial Engineering
Ansgar Steland
Institut fuer Statistik und Wirtschaftsmathematik, RWTH Aachen University

In data science applications such as sensor monitoring or financial portfolio
optimization the available data for modeling, computations and analysis of thecovariance matrix has often be done in a low-sample-size-high-dimensional (LSHD) regime.
Especially, if the dimension is larger than the sample size, classic methods fail andneed to be replaced by procedures which are designed for high-dimensional data.

Methods for testing and change detection in the covariance matrix can be based on recent
results on LSHD asymptotics of bilinear forms of the sample covariance matrix. Thisapproach allows us to detect and infer changes in averages of covariances of variables
or in the variance of projections. The theoretical results hold without the need toconstraint the dimension relative to the sample size. For the statistical estimationof unknowns one often needs a (large) learning sample. To circumvent this, we propose
in-sample estimators not requiring a learning sample.

Simulations show that the proposed methods work reliable for realistic mathematicalmodels. As a real world application the method is applied to analyze monitoring datafrom ozone sensors. The sensor data is compressed by projecting it onto sparse principal
directions obtained by a sparse principal component analysis (SPCA). It turns out that
the SPCA automatically learns the spatial locations of the sensors and leads to a spatial
segmentation. Analyzing the projections for a change-point provides a mean to detectchanges in the spatial dependence structure of the sensor network measuring ozone.

Do We Necessarily Have More Cost-Effective
Information with Big Data?*
Nitis Mukhopadhyay
Department of Statistics
University of Connecticut-Storrs

Abstract
It is a normally held belief that more data provide more information without mentioning what kind of
information one may be looking for. Such a feeling is widespread especially in the face of big data movement.
In order to hold a reasonable discourse, pros and cons both, I begin working under a certain stochastic models
of one kind or another and then try to grasp what additional useful information may entail as more data come in.
I will illustrate situations where a sentiment that more data may amount to more information is actually
nearly valid. But, I will also show that the same belief may not be entirely justified in other situations where therelative gain in information is not very substantial as one utilizes more available data. Indeed, the relative gain
in information may go down with more data accrued.
In this presentation, I will share preliminary ideas and show some analysis to validate my thoughts. It is
my earnest belief that there must be more to it than just using more and more available data simply because they
are out there!

Robust relative error estimation
Kei Hirose 1 and Hiroki Masuda 2
1

Institute of Mathematics for Industry, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

2

Faculty of Mathematics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

E-mail: hirose@imi.kyushu-u.ac.jp, hiroki@math.kyushu-u.ac.jp

In regression analysis, many analysts use the (penalized) least squares estimation, whichaims at minimizing the mean squared prediction error [5]. On the other hand, the relative(percentage) error is often more useful and/or adequate than the mean squared error. For
example, in econometrics, the comparison of prediction performance between diﬀerent stock
prices with diﬀerent units should be made by relative error; we refer to [9] and [11] among others.
Additionally, the prediction error of photovoltaic power production or electricity consumption
is evaluated by not only mean squared error but also relative error (see, e.g., [10]).
In relative error estimation, we minimize a loss function based on the relative error. An
advantage of using such a loss function is that it is scale free or unit free. Recently, several
researchers have proposed various loss functions based on relative error [9, 11, 1, 8, 2, 3]. In
practice, a response variable y(> 0) can turn out to be extremely large or close to zero. For
example, the electricity consumption of a company may be low during holidays and high on
exceptionally hot days. These responses may often be considered to be outliers, to which therelative error estimator is sensitive because the loss function diverges when y → ∞ or y → 0.
Therefore, a relative error estimation that is robust against outliers must be considered. Recently, Chen et al. [2] discussed the robustness of various relative error estimation procedures byinvestigating the corresponding distributions, and concluded that the distribution of least product relative error estimation (LPRE) proposed by [2] has heavier tails than others, implying that
the LPRE might be more robust than others in practical applications. However, our numericalexperiments show that the LPRE is not as robust as expected, so that the robustification of the
LPRE is yet to be investigated from the both theoretical and practical viewpoints.
To achieve a relative error estimation that is robust against outliers, this paper employs the
γ-likelihood function for regression analysis by Kawashima and Fujisawa [7], which is constructedby the γ-cross entropy [4]. An analysis of electricity consumption data is presented to illustrate
the usefulness of our procedure. For detail of the theoretical properties, algorithms, and Monte
Carlo simulations, please refer to Hirose and Masuda [6].

1

References[1] K. Chen, S. Guo, Y. Lin, and Z. Ying. Least Absolute Relative Error Estimation. Journalof the American Statistical Association, 105(491):1104–1112, Sept. 2010.
[2] K. Chen, Y. Lin, Z. Wang, and Z. Ying. Least product relative error estimation. Journalof Multivariate Analysis, 144:91–98, Feb. 2016.
[3] H. Ding, Z. Wang, and Y. Wu. A relative error-based estimation with an increasing numberof parameters. Communications in Statistics—Theory and Methods, 47(1):196–209, Nov.
2017.
[4] H. Fujisawa and S. Eguchi. Robust parameter estimation with a small bias against heavy
contamination. Journal of Multivariate Analysis, 99(9):2053–2081, Oct. 2008.
[5] T. Hastie, R. Tibshirani, and J. Friedman. The elements of statistical learning. Springer
Series in Statistics. Springer, New York, New York, NY, second edition, 2009.
[6] K. Hirose and H. Masuda. Robust relative error estimation. Entropy, 20(9):632, 2018.
[7] T. Kawashima and H. Fujisawa. Robust and Sparse Regression via γ-Divergence. Entropy,
19(12):608–21, Dec. 2017.
[8] Z. Li, Y. Lin, G. Zhou, and W. Zhou. Empirical likelihood for least absolute relative error
regression. TEST, 23(1):86–99, Sept. 2013.
[9] H. Park and L. A. Stefanski. Relative-error prediction. Statistics & Probability Letters, 40(3):227–236, 1998.
[10] D. W. van der Meer, J. Widén, and J. Munkhammar. Review on probabilistic forecastingof photovoltaic power production and electricity consumption. Renewable and Sustainable
Energy Reviews, 81(Part 1):1484–1512, Jan. 2018.
[11] J. Ye. Price Models and the Value Relevance of Accounting Information. SSRN Electronic
Journal, 2007.

2

Inconsistency of diagonal scaling under high-dimensional limit1

Tomonari Sei
Graduate School of Information Science and Technology,
The University of Tokyo.

1

Summary

We claim that diagonal scaling of a sample covariance matrix is asymptotically inconsistent if the ratio of the dimension to the sample size converges toa positive constant, where the population is assumed to be Gaussian with aspike covariance model. Our non-rigorous proof relies on the replica method
developed in statistical physics. In contrast to similar results known in theliterature on principal component analysis, strong inconsistency is not observed. Numerical experiments support the derived formulas.

2

Main results

Let x(1) , . . . , x(n) be independent and identically distributed according to
the p-dimensional Gaussian distribution with mean vector 0 and covariancematrix Σ ∈ Rp×p . Denote the (uncentered) sample covariance matrix by∑
S = (1/n) nt=1 x(t) x⊤(t) . We assume n ≥ p, which implies that S is positive
definite with probability one, unless otherwise stated.
Let R+ be the set of positive numbers. By a diagonal scaling theorem,
there exists a unique vector ŵ ∈ Rp+ such thatŴ S Ŵ 1 = 1,

(1)where Ŵ = diag(ŵ) and 1 = 1p = (1, . . . , 1)⊤ . In other words, all row sumsof the scaled matrix Ŵ S Ŵ are unity.
1

This study was supported by Kakenhi Grant Numbers JP17K00044 and JP26108003.

1

Let w0 be the population counterpart of ŵ, which means W0 ΣW0 1 = 1,
W0 = diag(w0 ). If p is fixed and n → ∞, a standard argument of asymptotic
statistics shows that ŵ converges almost surely to the true parameter w0
because S converges to Σ. However, if p is getting large as well as n, thenthe limiting behavior of ŵ is not obvious. We are interested in the behaviorof ŵ if αp := n/p converges to some α ∈ [1, ∞) as p → ∞.
In principal component analysis, this type of high-dimensional asymptotics has been deeply investigated. In particular, the angle between the first
eigenvectors of S and Σ converges to a non-zero value. Furthermore, thelimit becomes π/2 if α is less than a threshold. We call these phenomena
inconsistency and strong inconsistency, respectively.
We obtained similar conclusions for the diagonal scaling problem, at least
numerically. The following formula is derived with the help of the replica
method from statistical physics. See [1] for details.
Claim 1. Let Σ = I. Suppose that αp = n/p converges to some α ∈ [1, ∞)as p → ∞. Then we have
1− 3ŵ⊤ w0= √ 8α .
p→∞ ∥ŵ∥∥w0 ∥
1− 1
lim

(2)

2α

√
The right-hand side falls within (5 2/8, 1).
We also established formulas for a class of covariance matrices
Σ=Ω

11⊤
+ I,
p

(3)

where Ω is a positive constant expressing the signal-to-noise ratio.

References[1] Sei, T. (2018). Inconsistency of diagonal scaling under high-dimensional limit:
a replica approach, Preprint, arXiv:1808.05781.

2

Statistical Inference with Unnormalized Models
Takafumi Kanamori1
1

Tokyo Institute of Technology/RIKEN AIP

Abstract
In this talk, we propose a class of estimators for unnormalized models, and show that our estimators have the property of the eﬃciency in the parametric inference. The key concept of our method
is the density-ratio matching under Bregman divergences. This is the joint work with T. Takenouchi
of Future University Hakodate/RIKEN AIP, and M. Uehara and X.-L. Meng of Harvard University.

Statistical inference with parametric statistical models is an important issue in the fields of machinelearning and statistics. The maximum likelihood estimation (MLE) or its variants is often used for theparameter estimation. However, the direct application of the MLE to, say, the Boltzmann machine[6] is
not tractable. A diﬃculty comes from the calculation of the normalization constant. Its computational
cost is of the exponential order in the variable dimension. The problem of computational cost is common
to many probabilistic models, and several solutions for the estimation of unnormalized models have beensuggested; see Markov Networks[3], the Boltzmann machine (with hidden variables)[6, 1, 2], models in
independent component analysis[7], truncated distribution[8], exponential-polynomial distribution[4], andreferences therein.
When the normalization constant in a parametric model is not computationally tractable, we needtodealwith unnormalized statistical model, p(x; θ), θ ∈ Θ. Here, the integral of the probability density∫p(x; θ)dx is not necessarily equal to one, i.e., it may depend on the parameter θ. An approach for thestatistical inference with unnormalized models is to approximate the unnormalized model by a tractablemodel by the mean-field approximation, which considers a model assuming independence of variables [9].
Another approach such as the contrastive divergence [5] avoids the exponential time calculation by theMarkov Chain Monte Carlo (MCMC) sampling.
In this talk, we propose a class of estimators for unnormalized models. Our method does not requirecalculation of the normalization constant. The proposed estimator is defined by minimization of Bregman
divergence in which density-ratio matching is used. Our approach works for statistical models with both
discrete and continuous random variables. We show that in many case our estimator achieves the Fishereﬃciency and also possesses the convex property in the parameter for the unnormalized exponential family.
So far, estimators for unnormalized model rarely achieves the fisher eﬃciency except the estimator basedon the pseudo-spherical divergence on discrete sample space[10]. Our estimator is regarded as an extensionof [10] to the model on the continuous sample space. We show that not only pseudo-spherical divergencebut other Bregman divergences yield eﬃcient estimators for unnormalized models. Some theoreticalbackground and numerical experimenters will be shown in the presentation. The paper including thecontent of the presentation is now in preparation[11].

References[1] D. H. Ackley, G. E. Hinton, and T. J. Sejnowski. A learning algorithm for boltzmann machines.
Cognitive Science, 9(1):147–169, 1985.
[2] S. Amari, K. Kurata, and H. Nagaoka. Information geometry of Boltzmann machines. In IEEE
Transactions on Neural Networks, volume 3, pages 260–271, 1992.
1[3] J. Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal
Statistical Society. Series B (Methodological), 36(2):192–236, 1974.
[4] J. Hayakawa and A. Takemura. Estimation of exponential-polynomial distribution by holonomic
gradient descent. Communications in Statistics -Theory and Methods, 45(23):6860–6882, 2016.
[5] G.E. Hinton. Training Products of Experts by Minimizing Contrastive Divergence. Neural Computation, 14(8):1771–1800, 2002.
[6] Geoﬀrey E Hinton and Terrance J Sejnowski. Learning and relearning in boltzmann machines. MIT
Press, Cambridge, Mass, 1:282–317, 1986.
[7] Aapo Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of
Machine Learning Research, 6:695–708, 2005.
[8] N. L. Johnson, S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions. Wiley, New
York, 1995.
[9] M. Opper and D. Saad, editors. Advanced Mean Field Methods: Theory and Practice. MIT Press,
Cambridge, MA, 2001.
[10] Takashi Takenouchi and Takafumi Kanamori. Empirical localization of homogeneous divergences on
discrete sample spaces. In C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett,
editors, Advances in Neural Information Processing Systems 28, pages 820–828. Curran Associates,
Inc., 2015.
[11] M. Uehara, T. Kanamori, T. Takenouchi, and X. L. Meng. Unified eﬃcient estimation framework
for unnormalized models. in preparation.

2