「ベイズ統計の理論と方法」の補題4(2)の反例?

タイトルは釣り。
ベイズ統計の理論と方法」(渡辺澄夫)を読んでいて、2章でちょっと困っている。

ベイズ統計の理論と方法

ベイズ統計の理論と方法

他の本には書かれていないようなことが注意書きにたっぷり書かれていたりして、普通に通読するだけでもおもしろいが、やっぱり紙と鉛筆でじっくり楽しむのが本筋かな。

「本書を読むのに必要な予備知識は、大学初年度に習う線形代数微分積分だけで十分」「(それに含まれない基礎数学が必要な)場所においては重要概念について初等的に理解できるように導入部分を加えている」(まえがき)「本書ではルベーグ測度論を仮定しない」(p203)など、厳密さよりも前提知識のハードルを下げることを優先していることが伺える。
それはもちろん良いことなので全然構わないのだが、「これは書かれてない仮定がありそうだな~」ということがちょいちょいあるのが読んでいて困る。

一番明確なところを例に上げると、注意12(p36-37)には、分子分母にパラメータ集合 W 上で定義された関数を含む数式を受けて「本書では W としてコンパクト集合を考えていくので、その場合にはこの式は分母、分子ともに有限の値を取る」*1とある。この記述からその関数はおそらく連続なのだろうと推測できるが、この本では関数が連続かどうか全く触れられておらず、暗黙の仮定となっている。
まあこのくらいわかりやすかったら書いてなくても忖度できるのでなんとかなるんだが*2、2章の補題4(2)はそれがわからなくて困っている。

まずできるだけ簡単に補題4(2)に必要な記号を導入する。
\mathbb{R}^N 上で定義された真の確率分布 q(x) をパラメータ w\in W\subset\mathbb{R}^d を持つ確率モデル p(x|w) で表現する統計的推測を考える。
任意の x に対し q(x)=p(x|w_0) となる w_0\in W が存在するとき、q(x)p(x|w)実現可能という。
q(x)p(x|w) で実現可能なとき、対数尤度比関数 f(x,w) を以下のように定義する。*3

 \displaystyle f(x,w)=\log\frac{q(x)}{p(x|w)}

関数 g(x)q(x) で平均したものを \mathbb{E}_q[g(x)]=\int q(x)g(x)dx と書く。*4
ある定数 c_0>0 が存在して、任意の w\in W に対して、\mathbb{E}_q[f(x,w)]\geq c_0\mathbb{E}_q[f(x,w)^2] が成り立つとき、対数尤度比関数が相対的に有限な分散を持つという。
ここで補題4(2) は次のような命題である。

補題4(2)
q(x)p(x|w) で実現可能であれば、f(x,w) は相対的に有限な分散を持つ。

本では q(x) にも p(x|w) にも何の仮定もない。
W にもなにもないが、実は先の注意12はこの補題4(2)の証明を f(x,w)\approx0 の近傍で考えればよいということを言っていた。よって、命題には書かれていないものの W にはコンパクト性が仮定されていると思われる。
この状態で本に書かれている補題4(2) の証明のアウトラインを追いかけて、「ここ何の仮定もないと成り立たないよね」という部分を拾っていくと次のような反例が構成できてしまった。



真の分布: q(x)=1 \; (0\leq x\leq1)、パラメータ空間: W=\{w|0\leq w\leq1\}、確率モデル: p(x|w)=\frac1{z(w)}\exp(-x^{-w}), \;z(w)=\int_0^1\exp(-x^{-w})dx とする。0<\exp(-x^{-w})\leq \exp(-1)=1/e より z(w)有界である。
p(x|w=0)=q(x) かつ最適なパラメータの集合は W_0=\{0\} なので、q(x)p(x|w) によりユニークに実現可能である。

この確率モデルの対数尤度比関数は f(x,w)=\log \frac{q(x)}{p(x|w)}=\log z(w)+x^{-w}
その期待値は \mathbb{E}_q[f(x,w)]=\log z(w)+\int_0^1x^{-w}dx=\log z(w)+\frac1{1-w} となり、w<1 にて有限値を持つ。
ところが対数尤度比関数の2乗の期待値は
\displaystyle \begin{eqnarray}
\mathbb{E}_q[f(x,w)^2]&=&\int_0^1\{\log z(w)+x^{-w}\}^2dx\\
&=&\{\log z(w)\}^2+\log z(w)\int_0^1x^{-w}dx+\int_0^1x^{-2w}dx
\end{eqnarray}
となる。この第1項・第2項は有限値だが、第3項 \int_0^1x^{-2w}dxw\geq1/2 にて +∞ に発散することから、この対数尤度比関数は相対的に有限な分散を持たない。


この「反例」の q(x),p(x|w) はともに連続かつ W もコンパクトなので、即座に忖度できる範囲の仮定では足りないことがわかる。
\mathbb{E}_q[f(x,w)], \mathbb{E}_q[f(x,w)^2] が有限 or 有界とか? 要求される仮定としては使いにくすぎるし、それを満たす確率モデルというのも直感的にわかりにくい。

本に証明として書かれているのは、F(t)=t+e^t-1 という関数を天下りに考え、

 \displaystyle q(x)F\left(\log\frac{q(x)}{p(x|w)}\right)=q(x)\log\frac{q(x)}{p(x|w)}+p(x|w)-q(x)

から両辺積分して

 \displaystyle \int q(x)F\left(\log\frac{q(x)}{p(x|w)}\right)dx=\int q(x)\log\frac{q(x)}{p(x|w)}dx=\mathbb{E}_q[f(x,w)]

が得られる(p(x|w)-q(x)積分して消える)。
一方、平均値の定理から(と本には書いてあるが、F(t) をテイラー展開したときに0次と1次がともに消えるので、平均値の定理から派生したテイラーの定理を2次剰余項=F(t)に適用、といったほうが一般にはわかりやすいだろう)、 F(t)=\frac{t^2}2\exp(-t^*), |t^*|\leq|t| を満たす t^* が存在し、 f(x,w)\approx0 の近傍で

 \displaystyle \int q(x)F\left(\log\frac{q(x)}{p(x|w)}\right)dx\approx\int q(x)\cdot\frac12\left(\log\frac{q(x)}{p(x|w)}\right)^2dx=\frac12\mathbb{E}_q[f(x,w)^2]

より相対的に有限な分散を持つと言える、というもの。

自力ではとても思いつく自信がないトリッキーで見事な「証明」だが、最後の最後で突然えらく雑な評価に落とし込まれている。t^* とやらはこのままでは t=f(x,w) に依存してるので、そんな簡単に積分から消せない。
そこをちゃんと評価するには、 f(x,w)\approx0 の近傍という条件を使って |t^*|\leq|f(x,w)|<\epsilont^* を定数で押さえることで \exp(-t^*)>\exp(-\epsilon)積分の外に出し、近傍の外は有限開被覆有界性で処理かな。
でもそれをするには、上の積分が w を動かしたときにちゃんと収束しないといけないので、{}^\forall w_1\in W に対して \lim_{w\rightarrow w_1} p(x|w)p(x|w_1) に一様収束するって仮定とかあれば行けそう? この仮定があれば上の反例も弾けるはず(\lim_{w\rightarrow0}p(x|w) が p(x|w=0) に各点収束しない)。
そんなに変な仮定ではないとは思うけど、もっと弱い仮定でも通るんだろうか。統計の理論畑は全然見てきていないので、一般的な仮定がどれくらいなのかを全然知らないんだよなあ。

【追記】
その後つらつら考えてみたら、f(x,w) が発散せずに定義できるために、常に q(x)>0, p(x|w)>0 の仮定もあったほうがいい気がしてきた。
上に書いた一様収束の条件は、\mathbb{E}_q[f(x,w)] たちが W 上で連続を言うために使う。これでどうだ?
【/追記】

*1:この記述も問題があって、「分子・分母ともに有界」でないとわざわざコンパクト性を要求した意味がないし、その後に続く記述にも有界性が必要だったりする。

*2:数学的構造うんぬんを言うなら、当たり前の仮定だからこそちゃんと書いてほしいという気持ちもなくもない。

*3:本ではまず実現可能でない場合について定義するが、補題4(2) の説明には不要なので省略。

*4:本では確率変数 X を使って \mathbb{E}_X と表記しているが、notation を節約するために q を使った。こっちのほうが紛れも少ないし。

Randomized Response のベイズ推論(3): 変分ベイズ

Randomized Response はアンケートの回答をランダム化することで、個人の回答は伏せつつ平均などの統計量を得る手法の1つ。
前回までの記事で、ランダム化された回答から真の割合を最尤推定ベイズ推定(ギブスサンプリング)で推定する方法とそれらの実験を紹介した。今回は同じベイズ推定でも変分ベイズによる推論を行ってみる。
Randomized Response のギブスサンプリングによる推論までは普通に論文が見つかるが、変分ベイズ推論は探した範囲では見つけられなかったので、頑張って導出してみる。

Collapsed Variational Bayes 推論

記法は前回と同様。
D 通りの選択肢を持つ質問 Q に対する回答を確率変数 X\in\{1,\cdots,D\} で表す。
回答者は真の回答 X を知られたくないので、 X をランダム化した Y\in\{1,\cdots,D\} を集計者に返す。
このとき以下のモデルを考える。

 \begin{cases}
P(Y=j|X=i)=p_{ij} & i,j=1,\cdots,D, \; \sum_j p_{ij}=1 \\
P(X=i)=\pi_i & i=1,\cdots,D \\
P(\boldsymbol{\pi}) = \rm{Dir}( \boldsymbol{\pi}| \boldsymbol{\alpha} )
\end{cases}

ランダム化の確率 P(Y=j|X=i)=p_{ij} からなる計画行列を \boldsymbol{P}=\left(p_{ij}\right) と記す。
また \boldsymbol{\pi}=(\pi_i) にはパラメータ \boldsymbol{\alpha}=\left( \alpha_i \right) を持つディリクレ共役事前分布を入れる。

ユーザ n の真の回答を X_n、それをランダム化したものを Y_n で表す。
集計者は N 人の回答者のランダム化された回答の観測値 Y_1,\cdots,Y_N から \mathbb{E}[X] などを推定する。

f:id:n_shuyo:20210119180208p:plain:w320

変分ベイズ推論は事後分布に独立性の仮定(変分近似)を入れて、最終的に \log q(X_n)\log q(\boldsymbol{\pi}) などを計算して更新式を導出する

 q(X_1,\cdots,X_N,\boldsymbol{\pi}) \approx q(\boldsymbol{\pi})\prod_{n=1}^Nq(X_n)

しかしせっかくなので、より高い性能を期待して LDA(Latent Dirichlet Allocation) でおなじみの Collapsed Variational Bayes を使ってみよう(Teh+2007, Asuncion+2009)。
Collapsed Variational Bayes ではまず変分近似を以下のように弱める。

  q(\boldsymbol{X},\boldsymbol{\pi}) 
= q(\boldsymbol{\pi}|\boldsymbol{X})q(\boldsymbol{X})
\approx q(\boldsymbol{\pi}|\boldsymbol{X})\prod_{n=1}^Nq(X_n)

ただしこれ以降 X_1,\cdots,X_N\boldsymbol{X}Y_1,\cdots,Y_N\boldsymbol{Y}X_n 以外の X_* 全体を \boldsymbol{X}^{-n}=\{X_1,\cdots,X_{n-1},X_{n+1},\cdots,X_N\} で表す。
この仮定に基づき変分自由エネルギー \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right) を以下のように展開する。

 \begin{eqnarray}
\mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right)
&=& \mathbb{E}_{q(\boldsymbol{X},\boldsymbol{\pi})} [ -\log P(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{\pi}) ] - \mathcal{H}(q(\boldsymbol{X},\boldsymbol{\pi})) \\
&=& \mathbb{E}_{q(\boldsymbol{X})} [ \mathbb{E}_{q(\boldsymbol{\pi}|\boldsymbol{X})} [ -\log P(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{\pi}) ] - \mathcal{H}(q(\boldsymbol{\pi}|\boldsymbol{X})) ]
 - \mathcal{H}(q(\boldsymbol{X}))
\end{eqnarray}

ここで自由エネルギー \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right)q(\boldsymbol{X},\boldsymbol{\pi}) で最小化する代わりに、まず q(\boldsymbol{\pi}|\boldsymbol{X}) で最小化して、次に q(\boldsymbol{X}) で最小化する。q(\boldsymbol{\pi}|\boldsymbol{X}) で最小化したものを \mathcal{F}\left(q(\boldsymbol{X})\right) と書く。

 \begin{eqnarray}
\mathcal{F}\left(q(\boldsymbol{X})\right) 
&:=& \min_{q(\boldsymbol{\pi}|\boldsymbol{X})} \mathcal{F}\left(q(\boldsymbol{X},\boldsymbol{\pi})\right) \\
&=& \mathbb{E}_{q(\boldsymbol{X})} [  -\log P(\boldsymbol{Y},\boldsymbol{X}) ]
 - \mathcal{H}(q(\boldsymbol{X}))
\end{eqnarray}

この \mathcal{F}\left(q(\boldsymbol{X})\right) を、q(\boldsymbol{X})=\prod_{n=1}^Nq(X_n) に基づいて、各 q(X_n) で反復的に最小化するのが Collapsed Variational Bayes となる。
多項分布となる q(X_n) のパラメータを q(X_n=i)=\gamma_{ni} とし、c_k^{-n}=\sum_{m\neq n}I(X_m=k), \;\gamma_k^{-n}=\sum_{m\neq n} \gamma_{mk} とおくと(ただし I(\cdot) は指示関数)、その更新式は以下のように求められる。

 \begin{eqnarray}
\gamma_{ni} &\propto& \exp\left( \mathbb{E}_{q(\boldsymbol{X}^{-n})} [ \log P(\boldsymbol{Y},\boldsymbol{X}^{-n},X_n=i) ]\right)  \\
&\propto& \exp\left(  \log P(Y_n|X_n=i) + \mathbb{E}_{q(\boldsymbol{X}^{-n})} \left[ \log \int \pi_i \prod_k \pi_k^{\alpha_k+c_k^{-n}-1}d\boldsymbol{\pi}  \right]\right)  \\
&\propto& \exp\left(  \log p_{iY_n} + \mathbb{E}_{q(\boldsymbol{X}^{-n})} \left[ \log (\alpha_i+c_i^{-n})  \right]\right)  \\
&\approx& p_{iY_n}(\alpha_i+\gamma_i^{-n})
\end{eqnarray}

式の最後の近似は、LDA の CVB0 と同じくテイラー展開の 0 次の項のみに近似している。

実験

実装はこちら。

github.com

同じ反復法だが、データ件数 N でループする必要があるため、Gibbs Sampling より処理が10倍以上重い。そのため、推論の試行とサマリーを分離した上、試行の並列化までしている(苦笑)。
ここまでやっても、(N,\alpha) の組み合わせ 9通り×10000試行の推論には手元の環境で10時間かかった。

前回と同様に Randomized Response の計画行列 P(Y|X) は1個のパラメータ p で表現し、実験では p=1/5 を用いる。
 \displaystyle P(Y=j|X=i)=p_{ij}=\begin{cases}p + \frac1D(1-p)& (i=j) \\ \frac1D(1-p) & (i\neq j) \end{cases}

D=4 個の選択肢を持つ質問に対し、N=100,1000,10000 人のユーザが 1:2:3:4 の割合で回答(\boldsymbol{\pi}=(0.1\;0.2\;0.3\;0.4)^\top)p=1/5 の Randomized Response でランダム化された回答から推定した \hat\pi_1,\cdots,\hat\pi_4 の分布を見る。
事前分布のパラメータによる挙動の違いも見たいので、\alpha=1.0,\;0.1,\;0.01 を対称なディリクレ事前分布のパラメータとして用いる。

まずは手法ごとの違いを見るため、最尤推定、ギブスサンプリング、そして変分ベイズ(Collapsed Variational Bayes)に対して N=100,1000,10000 ごとの推定値の分布を見比べてみる。事前分布のパラメータは \alpha=1 のものを選んでいる。
今回は密度推定した線も書き加えてみたので、分布の形がよりわかりやすい。
 

N=100 N=1000 N=10000
MLE f:id:n_shuyo:20210204133030p:plain:w200 f:id:n_shuyo:20210204133034p:plain:w200 f:id:n_shuyo:20210204133037p:plain:w200
Gibbs f:id:n_shuyo:20210204130853p:plain:w200 f:id:n_shuyo:20210204130856p:plain:w200 f:id:n_shuyo:20210204130900p:plain:w200
VB f:id:n_shuyo:20210204130903p:plain:w200 f:id:n_shuyo:20210204130907p:plain:w200 f:id:n_shuyo:20210204130911p:plain:w200

最尤推定解は \pi_i によらず均等な散らばり具合で [0,1] からはみ出すが、ベイズ化によって \pi_i を確率値としてモデリングすると解消できることは前回までの記事で見たとおり。変分ベイズももちろんベイズモデルなので、同様に [0,1] に収まるように推定される。
上の図ではギブスサンプリングと変分ベイズの結果はそっくり同じに見えるが、実は縦軸・横軸の数値をよく見ると分かるように、変分ベイズの方が分散が小さい(山も高い)。
それがもっと明確に読み取れるように、\hat\pi_1\hat\pi_4 について箱ひげ図を書いて、特にそれらの散らばり具合の違いを見比べてみる。

N=100 N=1000 N=10000
f:id:n_shuyo:20210204152557p:plain:w200 f:id:n_shuyo:20210204152600p:plain:w200 f:id:n_shuyo:20210204152603p:plain:w200

すべての場合で 変分ベイズ(左) のほうがギブスサンプリング(右) より分散が小さい。横軸の下に書かれた推定値の標準偏差からも、グラフの見た目だけではなく実際に分散が小さくなっていることがわかる。これは推定値の精度が良くなっていることを意味している。*1

次にディリクレ事前分布のハイパーパラメータによる挙動の違いを見る。
N=10000 ではハイパーパラメータによる差異がほとんどない(データが多いと事前分布の影響が小さくなる)のは前回ギブスサンプリングでも見たとおりなので、N=100,1000 の図を掲載する。

N=100 α=0.01 α=0.1 α=1.0
Gibbs f:id:n_shuyo:20210204161937p:plain:w200 f:id:n_shuyo:20210204161941p:plain:w200 f:id:n_shuyo:20210204161944p:plain:w200
VB f:id:n_shuyo:20210204161959p:plain:w200 f:id:n_shuyo:20210204162002p:plain:w200 f:id:n_shuyo:20210204162005p:plain:w200
N=1000 α=0.01 α=0.1 α=1.0
Gibbs f:id:n_shuyo:20210204162023p:plain:w200 f:id:n_shuyo:20210204162026p:plain:w200 f:id:n_shuyo:20210204162030p:plain:w200
VB f:id:n_shuyo:20210204162047p:plain:w200 f:id:n_shuyo:20210204162051p:plain:w200 f:id:n_shuyo:20210204162054p:plain:w200

ギブスサンプリングでは \alpha を 1 より小さくすると事前分布の影響が強く出た推定になっていた。特に N=100,\alpha=0.01 の両端が上がった分布(つまり \hat\pi_i のどれか1つが 1 で残りが 0)が顕著だろう。*2
それに対し、変分ベイズではそうした極端さはほぼなく、直感に即した真値の周りの山形の分布となっている。

\alpha を小さくすることによる平滑化の低減(トレードオフとして分散の増加)を箱ひげ図で確認しよう。平均が理屈通りに振る舞うのは当然なので、今回は中央値に着目してみる。

α=0.01 α=0.1 α=1.0
N=100 f:id:n_shuyo:20210204165425p:plain:w200 f:id:n_shuyo:20210204165428p:plain:w200 f:id:n_shuyo:20210204152557p:plain:w200
N=1000 f:id:n_shuyo:20210204165531p:plain:w200 f:id:n_shuyo:20210204165534p:plain:w200 f:id:n_shuyo:20210204152600p:plain:w200

まず前回の記事で指摘しそびれていたこととして、ギブスサンプリングの \alpha<1 での極端な分布では中央値も真値から大きく離れていることが箱ひげ図から読み取れる。推定値の分散も非常に大きい。N=100,\alpha=0.01\hat\pi_4 では台が [0,1] しかないのに標準偏差が 0.430 と、一様分布(標準偏差 0.289)のほうがマシなレベル。
ギブスサンプリングで \alpha<1.0 のハイパーパラメータを使うのは難しそうだ。

変分ベイズでは、真値の周りの山型分布から想定できるとおり、中央値がちゃんと真値に近い値として推定できている。そして \alpha を小さくすると中央値はより真値に近づく。

このように変分ベイズ(Collapsed Variational Bayes)による推論は、最尤推定やギブスサンプリングより良い推定値が得られる。\alpha<1 でも直感的に良い推定値が得られる可能性が十分高く、欲しい精度に合わせてハイパーパラメータを厳選する甲斐がありそう。
実際の観測値に対する推定を行う場合、平均より中央値が外れてないことのほうが嬉しいので(確率 1/2 で中央値以上/以下になる)、その点でも変分ベイズによる推論は筋が良いとわかる。

同じ(真の事後)分布の推定を目的としているギブスサンプリングとの意外なほどの大きな差異が何によって生じているのかは興味深いところ。思い当たる要因は変分近似と CVB0(テーラー展開の0次近似)しかなく、どちらも原理的には精度を下げる効果しかないはずなのに、変分ベイズのほうが自明に良い結果を叩き出すのは、LDA-CVB0 の不思議な性能の良さと同じ根っこにつながっていそう(そう言えば CVB0 の理論的解釈とか誰かやってたような……)。

References

  • Teh, Yee W., David Newman, and Max Welling. A collapsed variational Bayesian inference algorithm for latent Dirichlet allocation. In Advances in Neural Information Processing Systems 19, 2007.
  • Asuncion, A., Welling, M., Smyth, P., and Teh, Y. W. On smoothing and inference for topic models. In Proceedings of the International Conference on Uncertainty in Artificial Intelligence, 2009.

*1:ギブスサンプリングはサンプルの平均を推定値とするが、その長さが小さいために散らばりが大きくなっているという可能性を確認するため、元のソースでは 200サンプルを burn-out して、続く 200サンプルの平均をとっていたのを 500サンプル burn-out & 500サンプル平均 にしてみたが、傾向は全く変わらなかった。

*2:それでいて \mathbb{E}[\hat\pi_i] はちゃんと真値に近くなるのがよくできているというかなんというか

Randomized Response のベイズ推論(2): ギブスサンプリング

Randomized Response はアンケートの回答をランダム化することで、個人の回答は伏せつつ平均などの統計量を得る手法の1つ。
前回記事では、回答の割合の推定量最尤推定で得る手順を紹介したが、割合の推定値が負になる可能性があることを示した。
shuyo.hatenablog.com
その問題を解消する方法の1つに、統計モデルをベイズ化するアプローチがある。今回はその中でも一番シンプルなギブスサンプリングによる推論を紹介する。

2値アンケートだとベイズ化の効果が見えにくいので、ベイズ化の前にまずは多値化する。
なお、i 番目の要素が a_i であるベクトルを \left(a_i\right)(i,j)要素が a_{ij} である行列を \left(a_{ij}\right) のように簡易的に表記している。

質問 Q に対する回答者の答えを確率変数 X で表す。
回答者は真の回答 X を知られたくないので、 X をランダム化した Y を集計者に返す。
集計者は N 人の回答者のランダム化された回答 Y_1,\cdots,Y_N から \mathbb{E}[X] などを推定する。

QD 通りの回答を持つとし、X, Y ともに 1,\cdots,D の値を取るとする。
X をランダム化した Y の確率を P(Y=j|X=i)=p_{ij} とし\left(\sum_j p_{ij}=1\right)p_{ij} を要素とする計画行列 \boldsymbol{P}=\left(p_{ij}\right) を考える。
X のモデルを P(X=i)=\pi_i としたとき、\boldsymbol{\pi}=\left(\pi_i\right) の推定値が欲しい回答の割合となる。

Y の観測値 \mathcal{D}=\left\{y_1,\cdots,y_N\right\} に対し、過程は省略するが、最尤推定\hat{\boldsymbol{\pi}} は以下のように求められる。ただし c_j=\left|\left\{ n | y_n=j \right\}\right|, \boldsymbol{c}=\left(c_j\right) とする。

\displaystyle \hat{\boldsymbol{\pi}}=\boldsymbol{P}^{-1}\boldsymbol{c}

計画行列 \boldsymbol{P}=\left(p_{ij}\right)D(D-1) 個のパラメータを持つが、自由度が多すぎてそのままでは使いにくい。
そこで真値が X=i であるとき、確率 pY=i、確率 1-p\{1,\cdots,D\} から一様乱数で選んだ値を Y とするランダム化メカニズムがよく採用される。
この場合、1つのパラメータ p に対し、計画行列は以下のように定められる。*1

\displaystyle p_{ij}=\begin{cases}p + \frac1D(1-p)& (i=j) \\ \frac1D(1-p) & (i\neq j) \end{cases} (式1)

最尤推定\hat{\boldsymbol{\pi}} は 0 から 1 の範囲に収まらないことがある。前回と同様にシミュレーションで確認しておこう。
実験に使ったコードはこちら:
github.com

4個の選択肢を持つ質問に対し、N(=100,1000,10000) 人のユーザが 1:2:3:4 の割合で回答するとき、ランダム化された回答から \hat{\boldsymbol{\pi}} を推定する。
p=1/5 の Randomized Response でランダム化された回答から推定した \hat\pi_1,\cdots,\hat\pi_4 をプロットした。

f:id:n_shuyo:20210119133717p:plain:w320f:id:n_shuyo:20210119133720p:plain:w320
f:id:n_shuyo:20210119133723p:plain:w320

N=10000 ではいい感じに推定できているが、N=100 では混然一体となってしまっており、\pi_i 間の大小さえうまく推定できるかどうかわからないことが読み取れる。
Randomized Response のパラメータ p を大きくすれば精度は改善できるので、集めたい N の大きさや最低限望む安全性などと相談して p を決める必要がある。

N=100 における各 \pi_i の平均・標準偏差、および 95%信頼区間も見ておこう。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.1 0.2 0.3 0.4
\rm{std}[\hat\pi_i] 0.21 0.21 0.22 0.22
95%区間 -0.3~0.5 -0.2~0.6 -0.1~0.75 0.0~0.85

これほど重なっていても、\mathbb{E}[\hat\pi_i]=\pi_i はちゃんと成立している。
分散(標準偏差)は i によらず一定だ。ここはベイズ版との比較ポイントになるので覚えておこう。
95%信頼区間は思いっきり負の領域にはみだしている。ベイズ化によってこれの改善を図る。

ベイズ版 Randomized Response

\boldsymbol{\pi} に共役事前分布としてディリクレ分布を入れて、その事後分布を推定する。

\displaystyle P(\boldsymbol{\pi}) = \rm{Dir}( \boldsymbol{\pi}| \boldsymbol{\alpha} ) = \frac{\Gamma\left(\sum_{i=1}^D\alpha_i\right)}{\prod_{i=1}^D\Gamma(\alpha_i)}\prod_{i=1}^D\pi_i^{\alpha_i-1}

\boldsymbol{\alpha}=\left( \alpha_i \right) はディリクレ事前分布のパラメータである。

事後分布の推定には方法はいくつかあるが、ここではまずギブスサンプリングを使った推論を導出しよう。
ユーザ n の真の回答を X_n、それをランダム化したものを Y_n という確率変数で表す。
また X_n 以外の X_* 全体を \boldsymbol{X}_{-n}=\{X_1,\cdots,X_{n-1},X_{n+1},\cdots,X_N\} で表す。

f:id:n_shuyo:20210119180208p:plain:w490
Randomized Response のグラフィカルモデル

このとき観測値 \mathcal{D}=\left\{Y_1=y_1,\cdots,Y_N=y_N\right\} に対する事後分布 P(\boldsymbol{\pi}|\mathcal{D}) をギブスサンプリングで求めるのに必要な、潜在変数 X_1, \cdots, X_N および \boldsymbol{\pi} の全条件付き事後分布を計算する。

\displaystyle \begin{eqnarray}
&& P(X_n=i|\boldsymbol{X}_{-n},\boldsymbol{\pi},\mathcal{D}) \\
&\propto& P(X_n=i,Y_n=y_n|\boldsymbol{\pi}) \\
&=& P(X_n=i|\boldsymbol{\pi}) P(Y_n=y_n|X_n=i) = \pi_i p_{iy_n}
\end{eqnarray} (式2)

\displaystyle \begin{eqnarray}
&& P(\boldsymbol{\pi}|X_1,\cdots,X_N,\mathcal{D}) \\
&\propto& P(\boldsymbol{\pi},X_1,\cdots,X_N,\mathcal{D}) \\
&=& P(\boldsymbol{\pi})\prod_{n=1}^N \left\{ P(X_n|\boldsymbol{\pi})P(Y_n|X_n) \right\}
\propto \prod_{i=1}^D \pi_i^{\alpha_i-1+d_i}
\end{eqnarray} (式3)

ただし d_iX_n の値が i である個数 d_i = \left|\{n|X_n=i\}\right| とする。
P(\boldsymbol{\pi}|X_1,\cdots,X_N,\mathcal{D}) はその式の形から  \rm{Dir}\left( \boldsymbol{\pi}| (\alpha_i+d_i) \right) とわかる。

X_1, \cdots, X_N\boldsymbol{\pi} の初期値を適当に決めて、(式2) と (式3) からのサンプリングを繰り返せばギブスサンプリングとなる。
(式2) は、X_n のサンプリングが他の \boldsymbol{X}_{-n} に依存しないことを表している。実際グラフィカルモデルを見ると、\boldsymbol{\pi} を縛れば X_1,\cdots,X_N は条件付き独立になる。
つまり X_1,\cdots,X_N 全体を一度にサンプリングするようなコードを書くことができる*2

実験に使ったコードはこちら:
github.com

実装ではサンプリングを400回繰り返し、最初の200回を捨てて、残る 200回のサンプル \boldsymbol{\pi} の系列の平均をとって割合の推定値としている(回数は適当w)。
事前分布のパラメータは {}^\forall \alpha_i=\alpha と対称とし、まずは \alpha=1 に対してデータの件数 N=100,1000,10000 に対して、Randomized Response の試行とギブスサンプリングによる推論を繰り返して得られた分布が以下のようになった。

f:id:n_shuyo:20210120155525p:plain:w240f:id:n_shuyo:20210120155534p:plain:w240f:id:n_shuyo:20210120155546p:plain:w240

まず N=10000 ではきれいな結果になっている。これは問題ないだろう。

N が小さくなると分散が大きくなり、最尤推定解が負になることがある問題を解決するためにベイズ化したわけだが、N=100 や 1000 のヒストグラムの横軸を見れば分かる通り、\hat\pi_1 などの分布もすべて 0 以上の範囲に収まっている。
わかりやすい N=1000 で見ると、\hat\pi_2,\hat\pi_3,\hat\pi_4 は左右対称な釣鐘型であるのに対し、\hat\pi_1 の分布は左の裾が 0 のところにストンと落ちる非対称な分布となっている。
N=100 になると分布の重なりは著しくなってしまうが、それでも分布の台が 0 から 1 の範囲に収まることはちゃんと守られている。

この通りベイズ化は期待通りの効果を発揮してくれた様子だが、別の問題も生じている。
N=100, \alpha=1.0 における推定値の平均 \mathbb{E}[\hat\pi_i] を見てみよう。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.18 0.22 0.27 0.32
\rm{std}[\hat\pi_i] 0.09 0.11 0.12 0.13
95%区間 0.06~0.41 0.07~0.48 0.09~0.54 0.11~0.59

分散が \pi_i ごとに異なり、また最尤推定解のそれより小さいなど嬉しいところもあるが、今見てほしいのは最尤推定では成立していた \mathbb{E}[\hat\pi_i]=\pi_i| が成立しなくなっている点である。
これはディリクレ事前分布によって正則化の効果が生じて、解が平準化されているためだ(0.1,0.2 の推定解は大きく、0.3,0.4 のは小さくなっている)。
N=10000 ではまたちゃんと \mathbb{E}[\hat\pi_i]=\pi_i| が成立しており、データが少ないと事前分布に引っ張られるという正しくベイズらしい振る舞いになっている。

では事前分布による正則化の効果を下げるために \alpha を小さくしてみたらどうだろう。N=100,1000,10000 それぞれについて事前分布のパラメータを \alpha=0.1 にしてみたときの推定解の分布が以下の図だ。

f:id:n_shuyo:20210120155522p:plain:w240f:id:n_shuyo:20210120155532p:plain:w240f:id:n_shuyo:20210120155542p:plain:w240

N=10000 では事前分布の違いを感じさせないが、N=100,1000 では \pi_1 が高い確率で 0 に縮退することがわかる。
N=100, \alpha=0.1 における推定値の平均や分散も確認する。

\pi_i 0.1 0.2 0.3 0.4
\mathbb{E}[\hat\pi_i] 0.13 0.19 0.28 0.4
\rm{std}[\hat\pi_i] 0.18 0.23 0.27 0.3
95%区間 0~0.69 0~0.81 0~0.89 0~0.93

\alpha を小さくしたことで、たしかに平準化の効果は押さえられ、\mathbb{E}[\hat\pi_i] は真の割合に近づいている。しかし 95% 区間を見れば分かる通り、一番大きい \pi_4 でさえも 0 に縮退する可能性が十分ある( PRML に書いてあった「関連度自動決定(ARD)」と同じ効果)。
0 に潰れる可能性があるということは、広い範囲の値を取りうるということで分散も大きくなっている。

このようにベイズ推論を用いる場合は事前分布をうまく決める必要がある、という至極当たり前な結論がでたところで*3、Randomized Response のグラフィカルモデルをもう一度よく見てみると、LDA の簡易版(ドキュメントが1個で、単語の生起確率が given)であることがわかる。
となるとせっかくだから LDA ではポピュラーな Collapsed Variational Bayes 推論も試してみたいよね! ということで次回記事に続く。

shuyo.hatenablog.com

Reference

  • Oh, Man-Suk. "Bayesian analysis of randomized response models: a Gibbs sampling approach." Journal of the Korean Statistical Society 23.2 (1994): 463-482.

*1:本稿では触れないが、local differential privacy における privacy budget が一定の範囲で推定値の分散が一番小さくなるのは(式1)の形で表現できることが示せる。

*2:\boldsymbol{X}_{-n} に依存しないので、観測データのループを書かなくていい。今回の Python 実装では、numpy.random.multinomial を選択肢の種類数である M 回呼ぶだけで済んでおり、データ数の N を増やしても実行時間は変わらない。

*3:今回はまだそこまで手を伸ばしてないが、経験ベイズ的アプローチでエビデンスが最大となるようにハイパーパラメータを決定してもおもしろい結果が出るかも?

Randomized Response のベイズ推論(1)

Randomized Response は、センシティブな質問を含むアンケートなどを実施するとき、個々人の回答を知らずにその統計量(yes の割合など)を得る手法の1つ。

例えばキセル乗車や浮気などの発生割合を調べたいとき、「あなたはキセル乗車をしたことがありますか?」とアンケートを取っても正直に答えてはもらえないだろう。
そうしたセンシティブな質問に対する個人の回答は隠しながら、調査者がほしい全体における割合を統計的に推定する手法がいくつかあり、Randomized Response はその元祖とも言える手法になる。

基本的な枠組みは次の通り。
質問 Q に対する回答者の答えを確率変数 X で表す。
回答者は真の回答 X を知られたくないので、 X をランダム化した Y を集計者に返す。
集計者は N 人の回答者のランダム化された回答 Y_1,\cdots,Y_N から \mathbb{E}[X] などを推定する。
この枠組みは untrusted curator model と呼ばれる。

この X, Y をどのように定めるかでランダム化のモデルが決まる。
一番基本的な Randomized Response では、yes/no で回答する質問 Q を想定し、XX=1 なら yes、X=0 なら no という2値の確率変数とする。
Y も同じ2値の確率変数で、X=1(yes) のとき、確率 p_{11}Y=1、確率 1-p_{11}Y=0 をとる。X=0(no) のときは、確率 p_{00}Y=0、確率 1-p_{00}Y=1 をとる。
パラメータ p_{11}, p_{00} は「回答を正直に送信する確率」と考えればわかりやすいかも。


\begin{cases}
P(Y=1|X=1)=p_{11} & 
P(Y=0|X=1)=1-p_{11} \\
P(Y=0|X=0)=p_{00} & 
P(Y=1|X=0)=1-p_{00}
\end{cases}

X のモデルを P(X=1)=\pi とおき、パラメータ \piY の観測値 \left\{y_1,\cdots,y_N\right\} から最尤推定することで、真の解答における X=1(yes) の割合を推定することができる。
y_n=1 の個数を a とおくと、尤度 L

\begin{eqnarray}
L&=&\prod_{n=1}^N P(Y=y_n)=P(Y=1)^a \cdot P(Y=0)^{N-a}\\
&=&
\left(P(Y=1|X=1)P(X=1)+P(Y=1|X=0)P(X=0)\right)^a \\
&&\cdot \left(P(Y=0|X=1)P(X=1)+P(Y=0|X=0)P(X=0)\right)^{N-a}\\
&=&
\left(1-p_{00}+(p_{11}+p_{00}-1)\pi\right)^a
\left(p_{00}-(p_{11}+p_{00}-1)\pi\right)^{N-a}
\end{eqnarray}

となり、これを最大にする \pi=\hat\pi

\displaystyle\hat\pi=\frac{p_{00}-1+a/N}{p_{11}+p_{00}-1}

と求められる。
この推定値 \hat\pi の平均と分散も計算できる。過程は省略すると、

\displaystyle\mathbb{E}[\hat\pi]=\pi
\displaystyle\rm{Var}[\hat\pi]=\frac{1/4-\{p_{00}-1/2-(p_{00}+p_{11}-1)\pi\}^2}{(p_{00}+p_{11}-1)^2N}

となり、\hat\pi は不偏推定量であることがわかる。

分散の式より、パラメータ p_{11}, p_{00} で推定値の精度が変わることがわかるが、同様に安全性(プライバシー保護の度合い)も変わる。例えば極端な場合を考えるとわかりやすい。
p_{11}=p_{00}=1 では、常に Y=X であり、集計者は真の回答を知ることができる。したがって推定値はそのまま真の解(標本平均)となるが、安全性は皆無である。
一方、p_{11}=p_{00}=1/2 では、Y は常に 1/2 の確率で 1 となり、X によらず独立である。つまり Y から X を知ることは不可能で完璧に安全だが、欲しい統計量も当然得られない。

このように推定の精度と安全性はトレード・オフの関係にあり、許容できる安全性の範囲で精度が高くなるように p_{11}, p_{00} を決める必要がある。
そのために必要な「安全性」を定量的に測る方法が局所差分プライバシー(local differential privacy)である。*1
untrusted curator model に基づくランダム化手法は Randomized Response 以外にもいくつかあり、そういった手法を同じ安全性(とみなす)で揃えたときの精度の差異などを検証するときにも差分プライバシーの考え方は役に立つ。
が、この記事ではそこは本題ではないので紹介だけにとどめておく。

シミュレーション

Randomized Response をシミュレーションで実験してみよう。
N 人のユーザのうち 3割が yes と答えるアンケート(\pi=0.3)について、p_{00}=p_{11}=2/3 の Randomized Response でランダム化された回答(1/3 の確率で嘘を返す)から \pi を推定してみる。
これを N=100,1000,10000 でそれぞれ推定してみた結果をヒストグラムにしたのが次の図である。
回答者数 N を増やすごとに理論値通りに分散が小さくなり、推定の精度が上がることわかる。

f:id:n_shuyo:20210114182952p:plain
f:id:n_shuyo:20210114182959p:plain
f:id:n_shuyo:20210114182955p:plain

また最尤推定解は実のところ線形方程式の解なので、\pi は割合であり、0 から 1 の間の値であることなど斟酌してくれない。実際、N=100 では推定値が負になることもあることがヒストグラムからもわかる。
これを解決する方法の1つが、次回の記事で紹介する推論のベイズ化である。

shuyo.hatenablog.com

References

  • Warner, Stanley L. "Randomized response: A survey technique for eliminating evasive answer bias." Journal of the American Statistical Association 60.309 (1965): 63-69.
  • N. Holohan, D. J. Leith, and O. Mason. “Optimal differentially private mechanisms for randomised response.” IEEE Trans. Inf. Forensics Security, 12(11):2726–2735, Nov 2017.

*1:差分プライバシーは上述の通り「ランダム化手法の安全度(プライバシー保護の度合い)を定量的に測る枠組み」のことなので、最近巷でまかり通っている「(局所)差分プライバシーを使ってユーザーのプライバシーを守る」という謳い文句はおかしい。「差分プライバシーによって安全性を保証されたプライバシー保護技術を採用している」とか言ってほしい。

Nreal Light のイメージトラッキング(画像タグ)実装方法

日本でも 12月1日から一般販売が始まるとアナウンスされた AR グラス Nreal Light。

www.moguravr.com

大昔に注文した開発キットもまだ届いてない人がいっぱいいるのに……とか、8月に先行して一般販売が始まっているはずの韓国で買えたという話がまだ聞こえてこないのに……とか、いろいろ懸念はあるけど、楽しみに待つことにしよう。

気を取り直して。
Nreal Light では画像タグを扱うことができる。
あらかじめアプリケーションに登録された画像が目の前にあるとき、その画像の位置・向き・大きさが位置とクォターニオンのベクトルで得られるので、それらを使ってオブジェクトを表示すると、画像タグにオブジェクトが乗っているように見える。

f:id:n_shuyo:20201110220855j:plain

Nreal の SDK に画像タグを扱う Demo シーン ImageTracking が含まれており、これをビルドすれば簡単に試すことができる。
Nreal Developer の Image Tracking のドキュメントはそのサンプルの解説であり、イメージトラッキングのプログラミング方法は書かれていないため、ドキュメントを読んでも自分のアプリケーションに画像タグ機能を組み込むことはできない。
この記事では Nreal の Demo を参考に、画像タグ認識機能を組み込む方法を説明する。


まず最初に Nreal Developer の Image Tracking のドキュメント にしたがって、認識したい画像タグを格納したデータベースを静的に作成し、コンフィグファイルに登録する。
簡単に説明すると、 Unity プロジェクトのアセットの適当なディレクトリに認識させたい画像たちを置き(複数ディレクトリに分散させることはできない)、それらの画像を選択してから右クリック > Create > NRSDK > TrackingImageDatabase でデータベースファイルが生成される。
アセットの Assets/NRSDK にNRKernelSessionConfig というコンフィグファイルがあるので、その項目 Tracking Image Database に先ほど生成したデータベースファイルをセットすればOK。

f:id:n_shuyo:20201110164353p:plain

データベースの Inspector にて、Width という項目が変更できるようになっているが、これがどういう意味を持っているのかはわからなかった。読める範囲のソースからは参照されておらず、Native DLL はソース非公開……。
Quality は画像に対して計算される値で、平坦な画像だと Quality が小さくなるので、画質というより画像に含まれる特徴量の割合を表しているように思われる。
ドキュメントには Quality が 65 以上となる画像をお勧めとあるが、65 より小さい画像だと認識しにくいという体感はないので、まあまあ大丈夫かな。

作成したデータベースはコンフィグに設定する。Assets > NRSDK の中にある NRKernelSessionConfig を開いて、Inspector にて項目 Tracking Image Database に先ほど作ったデータベースを指定する。
また Image Tracking Mode が ENABLE になっていることも確認しておこう(動的に ENABLE/DISABLE を切り替えることもできる。デモのソースを参照)。


次に、認識した画像にオブジェクトを追随させるために以下のようなコンポーネントを用意する。区別しやすいように MyTrackableImage という名前にしておいた。

using UnityEngine;
using NRKernal;

public class MyTrackableImage : MonoBehaviour
{
    public NRTrackableImage Image;

    public void Update()
    {
        if (Image != null && Image.GetTrackingState() == TrackingState.Tracking)
        {
            Pose center = Image.GetCenterPose();
            transform.position = center.position;
            transform.rotation = center.rotation;
        }
    }
}

画像にかぶせるオブジェクトを Instantiate して、この MyTrackableImage コンポーネントを紐づければ、視界の中で画像が動いたときにオブジェクトが自動的に追随してくれるようになる。
そのあたりを司るコードが以下の MyTrackingImageController だ。適当な空 GameObject でも作って、このコンポーネントを Add しておく。

using System.Collections.Generic;
using UnityEngine;
using NRKernal;

public class MyTrackingImageController : MonoBehaviour
{
    public GameObject TrackingImageVisualizerPrefab;

    private Dictionary<int, MyTrackableImage> m_Visualizers = new Dictionary<int, MyTrackableImage>();
    private List<NRTrackableImage> m_TempTrackingImages = new List<NRTrackableImage>();

    public void Update()
    {
        if (NRFrame.SessionStatus != SessionState.Running) return;
        NRFrame.GetTrackables<NRTrackableImage>(m_TempTrackingImages, NRTrackableQueryFilter.New);
        foreach (NRTrackableImage image in m_TempTrackingImages)
        {
            int idx = image.GetDataBaseIndex();
            m_Visualizers.TryGetValue(idx, out MyTrackableImage visualizer);
            TrackingState state = image.GetTrackingState();
            if (state == TrackingState.Tracking && visualizer == null)
            {
                Pose pose = image.GetCenterPose();
                visualizer = Instantiate(TrackingImageVisualizerPrefab, pose.position, pose.rotation, transform)
                    .AddComponent<MyTrackableImage>();
                visualizer.Image = image;
                m_Visualizers.Add(idx, visualizer);
            }
            else if (state == TrackingState.Stopped && visualizer != null)
            {
                m_Visualizers.Remove(idx);
                Destroy(visualizer.gameObject);
            }
        }
    }
}

あとはこの MyTrackingImageController のインスペクタで、適当に作ってプレハブ化しておいたオブジェクトを Tracking Image Visualizer Prefab にセットすれば完成。
ビルド& Nreal にインストールすれば、データベースに登録しておいた画像を認識した場所に先ほどセットしたオブジェクトが表示される。

f:id:n_shuyo:20201110225604p:plain

このサンプルコードでは、すべての画像に対して Tracking Image Visualizer Prefab に設定された一種類のオブジェクトを表示するが、ピカチュウの絵にはピカチュウの 3D モデルを、マリオの絵にはマリオを、みたいに画像の種類ごとに異なるオブジェクトを表示したいこともあるだろう。
NRTrackableImage.GetDataBaseIndex() によって認識された画像のインデックス番号が得られるので、番号によって分岐することで実現できる。
ただし、データベースの番号は作り直すたびに変わる可能性があるため(おそらくファイル名順)、画像のファイル名をキーにしたいところ。
以下のコードで認識した画像のファイル名を得ることができる。NRTrackingImageDatabase のインスタンスは Start() あたりで取得して private 参照できるようにしておけばいいだろう。

NRSessionConfig config = NRSessionManager.Instance.NRSessionBehaviour.SessionConfig;
NRTrackingImageDatabase imgDB = config.TrackingImageDatabase;

// NRTrackableImage image;
int idx = image.GetDataBaseIndex();
string name = imgDB[idx].Name;


あとは Nreal Light の画像トラッキングについての特徴や要望(文句?)をつらつら書いてみる。

Nreal が認識できる画像は、ある程度大きく、かつ自分の正面に正対していないといけない。壁に貼ったポスターの正面に立って見るイメージだ。
ニンテンドー 3DS の AR カード(53mm×85mm、いわゆる名刺サイズ)や一般的なトレーディングカード(63mmx88mm)などの大きさだと顔の前 10cm くらいまで持ってこないと認識しないし、机に置いた状態のように斜めの角度も認識しにくい。
一度正面から見て認識させると、少し離して、少し斜めになってもある程度追随して認識し続けることはできるが、簡単に見失うのであまり期待しないほうがいい。
ポケモンカードで AR デュエルだ! みたいのを実現するには、やっぱり普通に机に置いてあるカードを認識できるようになってほしいところ。

また、別々の画像なら同時に複数枚認識させることはできるが、同じ画像についてはどれか1つだけしか認識されない。
これもトレーディングカードの類を認識させたい人には悲報だろう。

Nreal がタグとして認識する画像は、アプリケーションをビルドする前に画像タグ用の静的なデータベースに登録しておく必要があり、ビルド後に追加することはできない(少なくとも NRSDK 1.4.8 以前は)。
認識できる画像の種類を後から動的に追加したい要件はとても多くありそうな気もするが、NRTrackingImageDatabase に Add メソッドは生えているものの、#if UNITY_EDITOR でくるまれており、エディタ上で対話的に作成することしか想定されていない様子。

QR コードの認識と読み取りができたら嬉しかったのだが、それもできない*1

Nreal のイメージトラッキングを体験して、「オブジェクトが画像にちゃんと重なってないなー」と感じたなら、効き目が右目なのかもしれない。右目を閉じて左目だけで見てみよう。
Nreal のトラッキングの座標はどうも左目が基準になっているようで、実は左目ではぴったり重なっているが、右目ではものすごくずれている。
このズレは画像が小さいほど顕著で、カードサイズだと以下の写真くらいズレてしまう。

f:id:n_shuyo:20201110223232j:plain

画像が十分大きいとズレは気にならないレベルに収まるので、その点でもやはり「正面のポスター画像」を想定した設計になっているのだろうな、と思わされる。

*1:QR コードを静的な画像として登録してみたら、位置の認識だけはできたが……それじゃないよねw

NRSDK(Nreal SDK) for Unity の新バージョン 1.4.8 の新機能を試す

AR グラス Nreal Light のアプリを開発するための NRSDK(Nreal SDK) for Unity のバージョン 1.4.8 がこの9月に出た。1.3.0 の次が 1.4.8 なのがちょっと謎。
Release Note のうち、機能追加にあたるのは以下の通り。

  • Added running state tips (temperature, battery, lost tracking notifications)
  • Added dynamic switch of 6DoF/3DoF/0DoF
  • Supported Nreal Dimmer
  • Added "IsTouching" API in NRInput
  • Adapted to smartphone screen resolution change when NRSDK is running
  • Adapted to Unity 2019.4 LTS

Unity 2019.4 LTS への対応を明示してくれるようになったのが地味に嬉しい。
Unity の特定のバージョン以外で NRSDK を使うと、アプリを起動してしばらくしたらメモリリークして落ちたりするので、Unity ガチャで NRSDK がちゃんと使えるバージョンを引くまで結構面倒だった。
2019.4 LTS と NRSDK 1.4.8 の組み合わせでは今のところかなり安定して使えている。

それ以外も地味な機能ばかり……。
Nreal のサイトに、まるですでに対応しているかのように書かれているハンドトラッキングは、いったいいつになったらサポートされるんだ。
AR グラスのインターフェースの本命はやっぱりハンドトラッキングだと思うので、今の画像タグや平面認識と同レベルでもいいから早いところ公開してほしい。

気を取り直して。
地味とはいえ、せっかくの新機能なので使い方をチェックしておきたいところ。
しかし Nreal のドキュメントは1日で読んで試し終わるくらいの分量で、新機能についての記述はなく。API Reference はカラッポ。
せめて NRSDK の github があればよかったのだが、いつまでたっても Coming Soon。

というわけで 1.3.0 と 1.4.8 の diff 見比べて、ソースを読み、新機能の使い方を洗い出してみた。
Nreal Dimmer はなんのことかわからず、IsTouching in NRInput はさすがに想像がつく(そして想像通り)。smartphone screen resolution change は純正 computing unit には多分関係ない?*1
そこで残りの running state tips と dynamic switch of 6DoF/3DoF/0DoF についてまとめる。

running state tips (temperature, battery, lost tracking notifications)

バッテリーの残量やグラスの温度について2段階の警告(Middle/High Level)を表示できる。トラッキング(SLAM State)は失っているかどうかだけで段階はない。

これを使うには、まず UI Canvas をシーンに追加して、そこに UI 部品の Text を2個(title, message), Image, Button を適当に配置する。Button は無くてもいいが、追加する場合は Nreal のドキュメント Controller の Building a Project with User Input にしたがって Canvas の構成をいじる。
この Canvas に "NR Notification Window" というコンポーネントを追加し、Icon/Title/Message/Confirm Btn にさきほどの Image, Text×2, Button をひもづける。*2

f:id:n_shuyo:20200928160829p:plain

NR Notification Window の Middle Level Info と High Level Info を開いて、Sprite/Title/Message に適当なアイコン画像やテキストを設定する。アイコン画像は NRSDK にそれっぽいのが含まれているので、そのまま使ってもいい。

f:id:n_shuyo:20200928161751p:plain

High と Middle の違いは、ソースやアイコン画像名(苦笑)を信じれば以下の通り。tracking では Middle Level Info しかつかわない。

Middle High
電池残量 40% 以下 30%以下
グラスの温度 45度以上 55度以上
ラッキング Lost ---

次に空のゲームオブジェクトをシーンに追加して(紛れなければ既存のオブジェクトでもOK)、Component に "NR Notification Listener" を追加する。
Low Power と SLAM State と High Temp のうち、使いたいものを選んで Enable にチェック、Notification Prefab に先ほどの NR Notification Window を追加したオブジェクト(今の場合は Canvas)を設定する。
ここで注意しないといけないのは、使わない通知項目も含めて3か所ともすべてに NR Notification Window を設定しないといけないこと。これを怠ると、NRNotificationListener.Awake() が NullReferenceException を吐いて通知機能が働かない(Nreal ぇ……)。

f:id:n_shuyo:20200928192318p:plain

これで完成。上述の条件を満たすと UI Canvas が表示され、各レベルに設定したアイコンやテキストが表示される。ボタンは High Level 時のみ表示され、これをクリックするとアプリケーションが終了する。
電池や温度は確認に時間かかるが、トラッキングはグラス前面のカメラをふさぐだけでロストしてくれるので、動作確認が楽。

f:id:n_shuyo:20200928193719j:plain

dynamic switch of 6DoF/3DoF/0DoF

サンプルの NRCameraRig をインスペクタで開いて、 "NRHMD Pose Tracker" コンポーネントの項目 "Tracking Type" から "Tracking 6/3/0 Dof" を選ぶことで、もともと静的にトラッキングモードを変更することはできていた。
バージョン 1.4.8 では この NRHMDPoseTracker に ChangeTo6Dof() などが生えて、これをアプリ実行中に叩くことでトラッキングモードを動的に変更できるようになった。使い道は……思いつかない(笑)。

この機能を使う最小サンプル。
以下のコンポーネントを作って、NRCameraRig に追加。

 public class ChangeTrackingMode : MonoBehaviour {
     private NRHMDPoseTracker tracker;
     void Start() {
         tracker = GetComponent<NRHMDPoseTracker>();
     }
     public void OnChangeTo3Dof() {
         tracker.ChangeTo3Dof();
     }
     public void OnChangeTo6Dof() {
         tracker.ChangeTo6Dof();
     }
 }

あとは UI Canvas に Change 6 DoF などのボタンを配置し、上で作った OnChangeTo6Dof() などに紐づければ OK。
3DoF にしたときの position は (0,0,0) になるみたいで、6DoF の状態で移動した後切り替えても初期位置に戻される。

*1:もし関係あっても、scrcpy の解像度が変わって嬉しいこと何もないだろうから、パスで。

*2:"NR Notification Window" を追加するのは別の GameObject でもよい。構成上、Canvas と "NR Notification Window" が一対一対応するので、Canvas に追加しておけばいいだろう。

AR グラス Nreal Light ファーストインプレッション

AR グラス Nreal Light の開発キットを1月に注文。
3月に届くはずが、コロナ禍のせい(だけかどうかわからないが)で遅れに遅れて8月にようやく到着。

f:id:n_shuyo:20200908163801j:plain

AR と言えば Microsoft の HoloLens が代表的な製品で*1、他に Magic Leap One などもあるけど、「電脳コイル*2を実現するには、メガネ型で、日常的にずっとかけていられて、子供が学校に持っていける価格でないといけない。

後述するが、Nreal Light(以降 "Nreal") もまだまだ全然ずっとかけていられる代物ではないけど、現行機(開発版含む)で唯一のメガネ型であり、グラス単体なら $500 程度と、ぎりぎり小学校に持っていけなくもない(例:牛革のランドセルは4~5万円)、という理想に一番近い(一応)、夢を見させてくれるデバイスである。
まあまだ夢なんだけど。

というわけで、まずはハードウェア周りのファーストインプレッションから。

かけてみて真っ先に気づくのが、ツルの幅が広すぎてずり落ちる。調整機能も見当たらない。
しょうがないので、ダイソーで「メガネストッパー」なる、メガネのツルに通すシリコン製の柔らかい爪を買ってきて Nreal のツルに頑張って付けた。
これが大正解。ずり落ちなくなって、しかも Nreal の重さを鼻だけではなく両耳の後ろにも分散できるようになって、掛け心地もかなり良くなった。電源を入れてなければ、1時間以上かけていてもしんどくない。
ただし、Nreal の左のツルからは本体であるコンピューティングユニットと接続するための約 1メートルの USB ケーブルが伸びていて、そのコネクタもふくめてメガネストッパーを通す必要がある。かなり時間のかかる作業で、ケーブルやその接続部分に負荷をかける可能性もあるので、ご興味のある方は自己責任で。

f:id:n_shuyo:20200908171014j:plain

映像の見え具合は、高精細で十分きれい。立体感もちゃんとある。
Nreal の映像の視野角は 52度あり、実際の見える範囲としては60センチ先の27インチディスプレイといったところ。それをはみ出した映像はぶった切られる。
現行機最大クラスの視野角だが、「拡張現実」と言い張るにはまだまだやっぱり狭い。現行 VR ヘッドセットで標準の 100~110度にはなんとか早いところ追いついてほしい。
またハーフミラー方式(実視界と CG をハーフミラーで合成)なためしかたないのだが、仮想の映像を透かしてその向こうの「現実」がうっすら見えちゃっている。いわゆるホログラム(特に映画やゲームにおける表現*3)を思い浮かべてもらうとわかりやすいだろう。

電脳コイル」では、電脳物質(仮想の物体やペット)かどうかをメガネをずらして確かめるというシーンがちょいちょいあるのだが、Nreal では仮想の映像と現実の区別がつかない心配はない。
視野角はともかく、「ホログラムっぽさ」はハーフミラー式ではどうしようもないので、それ以外の表示方式が安価かつ軽量で実現するのを期待するしかないだろう(ビデオパススルーか、網膜投影か、はたまたブレインマシンインターフェースか……って、そこまでいったらメガネいらん)。

f:id:n_shuyo:20200908173100j:plain

さらに 6DoF トラッキングの精度がまだ低いのか、体を動かすと映像がちょっとついてきてしまう。それどころか、自分が静止していても表示が微妙に動く。
Oculus Quest など現行の VR ヘッドセットでは、静止している仮想オブジェクトは自分がいくら移動してもその場に静止して見えるおかげで、物理的に質量を持って存在しているかのように感じさせてくれる。XR におけるトラッキングの重要性を改めて実感。
このあたりはソフトウェアアップデートによる向上を期待したいところだけど、Nreal のカメラは前面の前方向にしかついていないので、果たして。
とはいえ Qualcommチップセットにトラッキング技術を積んでくる気満々なので、こうした心配は今だけになりそう。

クアルコムVR/ARデバイスの新型リファレンスモデル スマホ接続と一体型モードを切り替え | Mogura VR
https://www.moguravr.com/qualcomm-vr-ar-reference-model/

多くの VR ヘッドセットと同様に近視の人が裸眼で見ることはできないので、視力矯正の必要がある。
形状的にリアルメガネと重ねるわけにもいかないので、視力矯正用レンズを付けられるようになっている。
標準で矯正用フレームのサンプルが付いていて、これをメガネ屋さんに持って行って自分用の矯正レンズを入れてもらうか、開発キットと一緒に売っている Lens Box($499) を購入する必要がある。Lens Box は度数 1.0 から 8.0 まで 0.5 刻みの 15×2 枚のレンズが入ってて、Nreal の内側にマグネットで簡単につけたり外したりできる。
周りのいろんな人に AR グラスを体験してもらうべく Lens Box もゲットしたのだが、コロナ禍のせいでその機会は当分難しそうだ。無念。

f:id:n_shuyo:20200908183535j:plain

この写真でレンズの上に少し平たい部分があるが、ここが主な放熱を担っているようで、稼働しているとかなり熱くなる。温度を測る機器を持っていないので正確な数字を出すことができないのだが、とても触ってはいられないほどである。
その状態でこのメガネを掛けると額のすぐそばに高熱源がある状態になり、連続利用は我慢して30分が限界だった(個人の感想です)。
メガネ側の仕事は映像のデコードとカメラのエンコードあたりか。そのため、ある程度熱を持つのは仕方ないが、顔に近い側で触れないほどの熱を持つのはさすがに困る。せめて顔から遠い側ならよかったのだが(Oculus Go のように)、前面にはカメラとレンズがあるから今の場所しかなかったということだろう。

連続使用を阻むのはこの熱問題だが、日常使用を妨げる問題がもう1つある。
Nreal を前から見ると、レンズの上半分にはカメラがある。その後ろには実は液晶パネルが下向きに入っており、ハーフミラーで実視界と映像が合成されるとてもシンプルな構造をしている。そのおかげで Nreal は他製品に比べて安価だが、レンズ部分にどうしても厚みが出てしまっている。

f:id:n_shuyo:20200820151923j:plain
問題は、このように上半分にいろんなモノが入っているため、実視界は下半分しかないという点である。Nreal をかけて道を歩くと、信号を見るために顔を上に向ける必要がある。車の運転はありえないだろう。
またハーフミラー式ならある程度は仕方ないものの、おそらくコストダウンか発熱を少しでも抑えるためあたりの事情で、実視界がかなり暗くなっている(その方が液晶パネルの輝度を上げなくても良い)。Nreal を掛けながら PC を操作するには、ディスプレイの輝度を最大近くしないと文字が読めないほどの暗さである。
視界の暗さは最初からサングラスかけてるつもりになればなんとかなる(かなあ……)としても、実視界の狭さは絶対に越えないといけないハードルなので、頑張ってほしいところ。
とはいえ、HoloLens がまさにその視界の広さ・明るさ(そしてレンズ部の薄さ)に技術をぶっこむことであの価格になっている様子を見ると、すぐに解決できる問題ではないんだろうなあ。

HoloLens光学系の謎に迫る
https://www.slideshare.net/AmadeusSVX/hololens-85758620

あとは、ここまでに書いたことと比べると些細な話だが、Nreal は公称 88g の「軽さ」をめちゃめちゃウリにしている。しかし実際に持って、掛けてみた感じはもうちょっと重く感じる。
また本記事ではほとんど触れなかった Nreal の本体であるコンピューティングユニット(こいつもめっちゃ熱くなる)は公称 140g なはずが、持った瞬間にわかる、それより絶対重い。
というわけで測ってみた。

f:id:n_shuyo:20200908191842j:plain

グラスの重さはケーブルを手で支えながら計測している。メガネストッパーも付けちゃった後なので、もう少し数字を割り引かれる可能性はあるが(鼻当ても外せるが、めんどくさくてやってない)、それでも 88g まで減るだろうか。
コンピューティングユニットはこれ以上外せるものは何もないので、140g は何かの間違いだろう(電池を除外した重量? 外せないが)。
それにしてもディスプレイもカメラもないのに重すぎ。Galaxy Note20 Ultra でやっと 208g なのに。

とここまで期待の裏返しで文句ばかり書いてしまったが、いろいろ割り引いた後でも十分楽しいデバイスである。
ソフトウェア面でもハードウェア面でも、VR における Oculus DK1 に相当すると言えば雰囲気は伝わるだろうか。Oculus Go/Quest のように人を選ばないレベルになるにはまだ数年かかるだろう。
そして「電脳コイル」のように日常的に使えるようになるのは何年後か。そういう夢を実物のデバイスを手にしながら語れるようになったのがすごいよね。

Nreal のソフトウェアや開発周りはまた別の記事で。

*1:単に AR というとポケモン GO なども含まれてしまうので、この場合 MR と言った方がいいのかも。

*2:電脳コイル」は、メガネ型のウェアラブルコンピュータが普及して、仮想の電脳世界が現実世界に交じり合った拡張現実が当たり前になった近未来(202X年)の地方都市を舞台に、主人公の小学生たちが日常の遊びの延長で電脳世界をハッキングしたり、破棄された電脳世界に隠された秘密を探したり……という名作アニメ。一言で説明するの、難しいな。https://www6.nhk.or.jp/anime/program/detail.html?i=coil

*3:スターウォーズのホログラムや、「ホライゾン ゼロ ドーン」のフォーカスでみるホログラフィックなどなど