ディスプレイの片眼用の画素数$d$と左右眼の全視野角$\theta$（図中の青）のみが判明している場合の角画素密度PPDの計算法をメモしておく． 一般に両眼視野角（左右片眼視野の論理積）は約120度と言われており，単眼視野角40度を加えたものが片眼の視野角100度（図中の緑と赤）となる． ビデオ透過型のHMDでは，この片眼の視野の一部をカバーしており，両眼視野については完全にカバーされていると仮定する． つまり片眼用のディスプレイがカバーする範囲は，
\[ 120/2+\theta/2 \]であり，全視野の角画素密度が一様であるとするとそのディスプレイの角画素密度は
\[ d/(120/2+\theta/2)=2d/(120+\theta) \]として算出できる．
EdelridのMega Julは駆動部分が無いにもかかわらず懸垂下降もできるオートロック式の軽量ビレイデバイス。Mega Jul→Mega Jul Sport→Mega Jul2と微妙な進化を遂げたけど、どこがどう変わったのかあまり情報がない、というか日本語ではMega Jul2の情報が殆どWebで見つけられない（2024年7月9日現在）ので分かる範囲内でまとめてみた。
Pythonスクリプト「報告書_20240523.pdf」 のようなファイル名を 「20240523_報告書.pdf」 に変換するスクリプトをChatGPT-4oに相談して作り理解したのでメモ． シェルスクリプトよりも理解し易い気がする．
import os
import re
# 変更したいファイルが存在するディレクトリを指定
directory = './' # カレントディレクトリを指定
# ファイル名のパターンを定義
pattern = re.compile(r'(.+)\_([0-9]{8})(.*)(\.pdf)')
# ディレクトリ内のすべてのファイルをループ
for filename in os.listdir(directory):
# ファイル名がパターンに一致するか確認
print(filename)
match = pattern.match(filename)
if match:
# 新しいファイル名を構築
new_name = f"{match.group(2)}_{match.group(1)}{match.group(3)}.pdf"
# ファイルのパスを構築
old_file = os.path.join(directory, filename)
new_file = os.path.join(directory, new_name)
# ファイル名を変更
os.rename(old_file, new_file)
print(f"-> {new_name}")
ファイル名は次のように分解できる。
（部分A）_（部分B）.（拡張子）
これを
（部分B）_（部分A）.（拡張子）
に変換したいだけだけど、（部分B）が「20240525」のような数字8文字だということを利用する。 （部分Aは）「(.+)」、（部分B）は「([0-9]{8})」でマッチさせる。どちらも丸括弧()で囲うことで、キャプチャグループという単位にすることができて、後でmatch.group(1)のように再利用できる。
Matplotlibのカラーマップ(Colormap）は便利で連続的な数値を滑らかに変化する色（疑似カラー）で可視化することができる． ただ，例外的な値（NaNや-1）の時だけ別の色で描きたいときはどうすれば良いのか？実はMatplotlibの標準的な機能で実現できる．実際に使ってみたのでソースコードをメモしておく．
その方法のポイントは，既存のカラーマップ（ここではjet）を少し変更することと，numpyのndarrayではなくてMaskedArrayというクラスに変換してplt.imshowに渡すこと． Colormapクラスを操作するのでそれに必要なモジュールやクラスをimportしておく．
import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap
次に描画する対象の配列を準備する． ここではnpzで保存された271$\times$271の画像で，各画素には正の実数が格納されているとする． 確認のためshapeとtypeをprintしておく．
array = np.load("positive.npz")["arr_0"]
print(array.shape)
print(type(array))
(271, 271, 1) <class 'numpy.ndarray'>
ただし，この正方形の画像の右上半分には有効な値は入っておらず，ここでは便宜的に0が格納されている． これを実数値を虹色で表すカラーマップjetを使って描画するとこうなる．
plt.imshow(array, cmap = "jet")
ただ，これだと0付近の値（青）と値が無いことを意味する数値ラベル0が殆ど同じ色になってしまう．そこで値が無い箇所だけグレーで表すことを考える． そのためにMaskedArrayというクラスを使ってarrayの値が0のときだけ無効であるという情報を保持させる．
marray = np.ma.array(array, mask=array == 0) type(marray)
(271, 271, 1) <class 'numpy.ma.core.MaskedArray'>
mask=array == 0の部分が解りにくいので詳しく説明する．arrayが27$\times$27の2次元配列なのでarray == 0とするとbroadcastの効果によりスカラーの0が全要素0の27$\times$27配列に変換される． ==の演算子によってarrayの要素が0のときだけTrue，そうでないときがFalseの値を要素に持つ27$\times$27の配列が生成される．maskにはそれが代入されている．
カラーマップjetを少し改造して，無効な値のときの色をset_bad関数でグレーに設定した上で， ndarrayの代わりにこのMaskedArray型の画像marrayをimshowに代入すると自動的に色分けしてくれる．
jet = plt.get_cmap("jet").copy()
jet.set_bad("gray")
plt.imshow(marray, cmap = jet)
GoProのようなアクションカムで動画撮影した場合，そのひとコマを写真として保存して残しておきたいことが良くある． 大抵両手が塞がってるのでカメラのシャッターを押せないことが多いから． でも，殆どのアプリでは動画から静止画を保存していくと撮影日時を埋め込んでくれないので，ファイル名順になり順番はめちゃくちゃになる．
動画の撮影日時にそのフレームまでの再生時間を足して静止画の撮影日時とするプログラムがあれば，順番がめちゃくちゃにならずに済む． MP4動画には撮影日時が埋め込まれており，JPEG画像にもExifという撮影に関する情報を埋め込むことができる． そこでPythonで動画を再生しながら適当な箇所で一時停止して撮影日時付きの画像を保存できるツールを作ったので紹介する．
このプログラムのポイントは，MP4の撮影日時情報の取得にffmpegのprobeという機能を使い， JPEG保存後にExif情報だけをpiexifというライブラリで上書きするところ． それ以外は基本的なライブラリしか使わない．ここでは動画ファイルをtest.mp4としておく．
import sys import cv2 import piexif import datetime as dt import ffmpeg fname = "test.mp4"
次は動画キャプチャの骨格部分．main関数になってるので最下部のmain関数の実行部分がこのソースの本体． 基本的にOpenCVで動画再生して必要な箇所をキーボード入力で一時停止，キャプチャなどするための機能だけ． もう少しシンプルにしたかったけど一応動作してるのであきらめた． 重要な箇所は「# MP4の撮影日時取得」，「# 動画再生時間の取得」，「# JPEG保存とExif情報の上書き」の部分だけ．ここだけ簡単に解説．
def main(argv):
# MP4の撮影日時取得
video_info = ffmpeg.probe(fname)
date_time = dt.datetime.strptime(video_info['format']['tags']['creation_time'],"%Y-%m-%dT%H:%M:%S.%fZ")
print("creation_time: " + date_time.strftime("%Y:%m:%d %H:%M:%S"))
cap = cv2.VideoCapture(fname)
n_frames = cap.get(cv2.CAP_PROP_FRAME_COUNT)
ret, frame = cap.read()
print("number of frames: "+ str(n_frames))
is_playing = True
cv2.namedWindow("win")
while True:
# 動画再生時間の取得
play_time = int(cap.get(cv2.CAP_PROP_POS_MSEC))
frame_id = cap.get(cv2.CAP_PROP_POS_FRAMES)
play_datetime = date_time + dt.timedelta(microseconds=play_time*1000)
# key input
key = cv2.waitKey(1)
# ESC
if key == 27:
break
elif key == 13: # Enter
is_playing = not is_playing
elif key == 32: # space
frame_id += 1
is_playing = False
elif key == 8: # back space
frame_id -= 1
is_playing = False
elif key == ord('c'): # capture !!
# JPEG保存とExif情報の上書き
cv2.imwrite("%06d.jpg" % frame_id, frame)
insertTime("%06d.jpg" % frame_id, play_datetime)
if is_playing:
if 0 <= frame_id and frame_id < n_frames:
ret, frame = cap.read()
else:
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)
ret = True
else:
if 0 <= frame_id-1 and frame_id < n_frames:
cap.set(cv2.CAP_PROP_POS_FRAMES, frame_id-1)
ret, frame = cap.read()
ret = True
if ret:
showFrame("win", frame, is_playing, frame_id, play_time, play_datetime)
else:
break
cap.release()
cv2.destroyWindow("win")
まず「# MP4の撮影日時取得」の部分． ffmpeg.probe関数でファイル名を直接していすると必要な情報が取得できる． ただ，撮影日時以外の様々な情報が階層構造をもつ辞書形式になっているのでvideo_info['format']['tags']['creation_time']で取得すると中身が文字列として得られる．それをdatetime.datetime.strptime関数で指定した書式として解釈して日時情報に変換している．これで動画の撮影日時が得られた．
# MP4の撮影日時取得 video_info = ffmpeg.probe(fname) date_time = dt.datetime.strptime(video_info['format']['tags']['creation_time'],"%Y-%m-%dT%H:%M:%S.%fZ")
次に「# 動画再生時間の取得」の部分． ここではOpenCVのcv2.VideoCapture.get関数で現在のフレームの再生時間を取得している． この時間をさっき取得した動画撮影日時date_timeに加えることで動画の各フレームの撮影日時を計算している． ミリ秒とマイクロ秒の変換を間違えないように．
# 動画再生時間の取得 play_time = int(cap.get(cv2.CAP_PROP_POS_MSEC)) frame_id = cap.get(cv2.CAP_PROP_POS_FRAMES) play_datetime = date_time + dt.timedelta(microseconds=play_time*1000)
最後は「# JPEG保存とExif情報の上書き」の部分． 写真の保存は単純に現在のフレームをcv2.imwriteで保存するだけ． その直後にその画像ファイルに対して自作のinsertTime関数を使ってExif情報を加えている．
# JPEG保存とExif情報の上書き
cv2.imwrite("%06d.jpg" % frame_id, frame)
insertTime("%06d.jpg" % frame_id, play_datetime)
insertTime関数は次のとおり． piexif.load関数で直接JPEGファイル名を指定して一旦Exif情報を取り出して，その必要箇所だけ変更して再度保存する方法をとっている．piexif.ExifIFD.DateTimeOriginalには，Exif情報の中のDateTimeOriginalという撮影日時を表す項目のIDで36867が入ってる． piexif.ExifIFD.SubSecTimeOriginalには，SubSecTimeOriginalという撮影日時の秒単位よりも小さい時間を表す項目のIDで37521が入っている． この2つの情報をdate_time.strftimeで指定された書式にして追加している．
def insertTime(fname, date_time):
exif_dict = piexif.load(fname)
exif_ifd = {piexif.ExifIFD.DateTimeOriginal: date_time.strftime("%Y:%m:%d %H:%M:%S"),
piexif.ExifIFD.SubSecTimeOriginal: date_time.strftime("%f"),
}
exif_dict['Exif']=exif_ifd
exif_bytes = piexif.dump(exif_dict)
piexif.insert(exif_bytes, fname)
showFrame関数は現在のフレームを表示する関数． cv2.putTextで各種情報を同時に表示している．
def showFrame(wname, frame, is_playing, frame_id, play_time, date_time):
dsp = cv2.resize(frame, (frame.shape[1]//2, frame.shape[0]//2))
cv2.putText(dsp,
text=str("%s: %06d-th frame, %06d msec, %s"
% ("play" if is_playing else "pause", frame_id, play_time, date_time.strftime("%Y:%m:%d %H:%M:%S.%f"))),
org=(20, 20),
fontFace=cv2.FONT_HERSHEY_SIMPLEX,
fontScale=0.5,
color=(255, 255, 255),
thickness=1,
lineType=cv2.LINE_4)
cv2.imshow(wname, dsp)
プログラム本体をmain関数化しておくための行． これのおかげで関数の定義順序を自由にできる．
if __name__ == "__main__":
sys.exit(main(sys.argv))
意外とサンプルが存在しないmatplotlibを使った数直線の描き方．
import numpy as np import matplotlib.pyplot as plt
必要なライブラリはこれだけ．
b = np.array(["$B_1$", "$B_2$", "$B_3$", "$B_4$", "$B_5$", ]) x = np.array([-0.86, 0.90, -0.88, -0.27, 1.10]) y = np.zeros_like(x) x_ano = x + 15 x_ano[0] += 10 # B1 x_ano[2] -= 10 # B3
B_1，B_2はデータ点につけるラベル．x_anoはラベルの水平位置で$B_1$と$B_3$の位置（b[0]とb[2]）だけ重なってしまうので手動で左右にずらしている． こういうのを自動でやってくれるとありがたいけどそういう機能は見当たらないのであきらめた．
#数直線 figsize_px = np.array([1024, 102]) dpi = 200 figsize_inch = figsize_px / dpi fig,ax=plt.subplots(figsize=figsize_inch,dpi=dpi)
ここは普通にplt.subplotsを読んで問題ない．画素数をしていしたかったのでこんな感じに．
# メモリと枠線を消す ax.tick_params(labelbottom=True, bottom=False) ax.tick_params(labelleft=False, left=False) plt.gca().spines['right'].set_visible(False) plt.gca().spines['top'].set_visible(False) plt.gca().spines['bottom'].set_visible(False) plt.gca().spines['left'].set_visible(False)
数直線を描くためにmatplotlibデフォルトの枠線とメモリを全部消去する．
# 自前でメモリと水平軸を描画 xmin, xmax= -1.2,1.2 ys = 0.005 yl = 0.01 xs = np.linspace(int(xmin),int(xmax), int(xmax)-int(xmin)+1) xl = np.linspace(int(xmin),int(xmax), (int(xmax)-int(xmin))*10+1) plt.hlines(y=0, xmin=xmin, xmax=xmax, zorder=1) #横軸 plt.vlines(x=xs, ymin=-yl,ymax=yl, zorder=1) #目盛り大 plt.vlines(x=xl, ymin=-ys, ymax=ys, zorder=1) #目盛り小
ここが数直線のほぼすべて．水平線と垂直線の組み合わせ手で自前で軸を描いている．
# データ点と注釈を描画
plt.scatter(x,y,c='r') #散布図
for i in range(x.shape[0]):
ax.annotate('%s' % b[i],
size=12,
xy=(x[i],y[i]),
xytext=(x_ano[i], y[i]+25.0),
textcoords='offset points',
arrowprops=dict(arrowstyle="->",connectionstyle="arc3,rad=.3")
)
plt.savefig('plot_out.png', bbox_inches="tight", transparent=True)
データ点はscatterで注釈はannotateで描くだけ．
「サーストンの一対比較法を理解する」シリーズ． 前回「サーストンの一対比較法って何？」で，サーストンの一対比較法がどんな手法なのか，他の方法との違いとか欠点とか説明したので， 今回は具体的にサーストンの一対比較法でどのようにデータを収集して，解析するか説明する．
5種類のビール（ビール1～ビール5）の「キレ」を数値化するために100人の被験者にアンケートを取ってみたとする． アンケート方法は，2種類のビールを飲み比べて『ビール1よりビール2の方がキレが高いですか？』と問うもの． どちらも同じとか，回答拒否とかの選択肢はない．必ずどちらが「キレ」が高いかを回答しないといけない． だから『ビール1よりビール2の方がキレが高いですか？』と聞かれたらYesかNoのどちらかを回答する． Yesは「「ビール1よりビール2の方がキレが高い」という意味になり， Noは「ビール2よりビール1の方がキレが高い」という意味という意味になる．
\[ \def\S#1#2{\small B_{#1}\gt B_{#2}} \begin{array}{ccccc} \hline \hline & \S{1}{2} & \S{1}{3}& \S{1}{4}& \S{1}{5}& \S{2}{3}& \S{2}{4}& \S{2}{5}& \S{3}{4}& \S{3}{5}& \S{4}{5}\\ \hline \rm 被験者1& \rm No & \rm Yes & \rm Yes & \rm No & \rm No & \rm Yes & \rm Yes & \rm Yes & \rm No & \rm No\\ \rm 被験者2& \rm Yes & \rm No & \rm Yes & \rm No & \rm Yes & \rm Yes & \rm No & \rm Yes & \rm No & \rm Yes\\ \rm 被験者3& \rm No & \rm Yes & \rm No & \rm No & \rm Yes & \rm Yes & \rm Yes & \rm Yes & \rm No & \rm No\\ \rm \vdots& \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ \rm 被験者100& \rm No & \rm Yes & \rm No & \rm No & \rm No & \rm Yes & \rm No & \rm No & \rm No & \rm No\\ \hline \rm Yesの数& 6 & 46 & 28 & 2 & 99 & 87 & 31 & 27 & 6 & 8\\ \rm Yesの割合&0.06 & 0.46 & 0.28 & 0.02 & 0.99 & 0.87 & 0.31 & 0.27 & 0.06 & 0.08\\ \hline \end{array} \]この表は，集めたデータ例．「$\S{1}{2}$」は「ビール1はビール2よりキレが良い」という意味．それに同意すれば「Yes」となる． 「ビール2はビール1よりキレが良い」は「ビール1はビール2よりキレが良い」という問の逆なので省略している． だから，ビール5種類から2種類を抽出する組み合わせ数が$_{5}C_{2}=10$なので10列のデータになっている． 行数は，被験者数分だけある．
サーストンの一対比較法で直接用いるのは，各問いに対する「Yesの割合」だけ．誰がどう回答したかは関係ない． だから順序効果を考慮して問いの順序を無作為にするとか，あえて逆の質問をして逆にカウントするとか，工夫をする余地がある．
得られたデータの解析法は3ステップ． まず上で計算した「Yesの割合」を使って次の表を作る． この表は「ビール$j$はビール$k$よりキレが良い」と回答する割りあい（確率）を表している． ちょうどリーグ戦の勝率表に似ている．
\[ \begin{array}{rrrrr} \hline \hline j \backslash k & 1& 2& 3& 4&5\\ \hline 1& 0.50& 0.06& 0.46& 0.28& 0.02\\ 2& 0.94& 0.50& 0.99& 0.87& 0.31\\ 3& 0.54& 0.01& 0.50& 0.27& 0.06\\ 4& 0.72& 0.13& 0.73& 0.50& 0.08\\ 5& 0.98& 0.69& 0.94& 0.92& 0.50\\ \hline \hline \end{array} \]この表を作るときには，$j$と$k$には1から5の整数を一通り入れて考えて，値は「Yesの割合」からもってくればいい． そうやって数値を入れていくと対角成分よりも上がわの三角の範囲しか値が入らないことが分る． 実は上の三角形部分は，$j\lt k$の場合だけで$j\gt k$の場合はない．つまり「$\S{1}{2}$」は有っても「$\S{2}{1}$」はない． $j\gt k$は下の三角形部分にいれるべきだけどこの部分のアンケートを取らなかったのは$\S{j}{k}$がYesとなる確率が分れば， $\S{k}{j}$がYesとなる確率は次の式で自動的に分かるから．
\[ P(\S{k}{j}) = 1 - P(\S{j}{k}) \]例えば，「ビール1はビール2よりキレが良い」がYesである確率は$P(\S{1}{2})=0.29$だけど， 「ビール2はビール1よりキレが良い」がYesである確率は自動的に$P(\S{2}{1})=1 - 0.29=0.71$になる． これはアンケートで「Yes」か「No」以外の回答を許さないから必ずそうなる．
もう一つデータが無い箇所があり，それは対角成分である$j=k$の場合． 例えば「ビール1よりビール1はキレが良い」と言われてもYesでもNoでもないからそもそもそういうアンケートに含めていない． この場合の確率は，無理やりにでも「Yes」か「No」どちらか答えろと言われたら「Yes」の割合は半分になるはずなので0.50とする． だから対角成分には全部0.50が埋まっている．
ここまではアンケート結果から確率を求めたりして何となく何がしたいのか想像できたかもしれないけど，次の手順は理論を知らないと意味不明かもしれない． でもとりあえずここでは理論は無視してどう計算するかだけ説明する．手順自体はエクセル等を使えば簡単． 何をするか一言で言うと標準正規分布の累積確率$\Phi(z)$が$\Phi(z)=P(\S{j}{k})$となる$z=z_{jk}$を求めるということ． つまり標準正規分布の累積確率密度関数を$\Phi(z)$は
\[ \Phi(z)=\int_{-\infty}^{z} \frac{1}{\sqrt{2\pi}} \exp^{-\frac{x^2}{2}}dx \] と表せる．この積分の中が標準正規分布．標準正規分布とは正規分布の関数に$\mu=0$と$\sigma=1$を代入したもの． 累積確率とは確率密度関数を積分したものだけど「下側」とは$-\infty$から$z$まで積分したもの．下のグラフで言えば青い部分の面積にあたる．
確率密度関数を積分した面積は確率そのものになる．ここではこの確率が$P(\S{j}{k})$となるような$z$を$z_{jk}$として，この$z_{jk}$を求める．
\[ \Phi(z_{jk})=P(\S{j}{k}) \]この計算はエクセルやGoogleスプレッドシートのnorm.s.inv()関数を使うと簡単に計算できる． 例えばnorm.s.inv(1-0.06)=1.55と計算できる． セル全部を計算した結果がこれ．
\[ \begin{array}{rrrrr} \hline \hline j \backslash k & 1& 2& 3& 4&5\\ \hline 1& 0.00& 1.55& 0.10& 0.58& 2.05\\ 2& -1.55& 0.00& -2.33& -1.13& 0.50\\ 3& 0.10& 2.33& 0.00& 0.61& 1.55\\ 4& -0.58& -1.13& -0.61& 0.00& 1.41\\ 5& 2.05& -0.50& -1.55& -1.41& 0.00\\ \hline \hline \end{array} \]計算結果をよく見ると対角成分（$j=k$）は0で，対角成分を挟んで反対側は正負反対になっている． これは標準正規分布が原点に対して左右対称だから．
最後はさらに理論を知らなければ意味不明だけど，計算は簡単． 上の表の各行を平均するだけ． ちなみに添字がいつのまにか$i$になっているのは理論の方と整合とってて多少の理由がある．
\[ \hat{\mu}_i=\displaystyle\frac{1}{5}\sum_{j=1}^5 \hat{z}_{ji} \]この計算によってこの数直線のような解析結果が得られる． ここで注意が必要なのはこのデータ点の重心が原点，分散が1になるように数値が決まっているのでデータ点1点の数値そのものには意味がない． データ点どうしの大小比較や2組以上のデータ点間の距離の比については意味がある数値になっている（グラフの描き方は「Pythonでラベル付きの数直線を描く」を参照）．
「サーストンの一対比較法を理解する」シリーズとして今回はサーストンの一対比較法がどんな手法なのか，他の方法との違いや欠点はどんなものなのか解説する． サーストンの一対比較法は，感覚を数値化する方法． 感覚とは，例えば，「イケメン度合い」とか「美味さ」，「肌触り」とか一見，数値化が難しそうな感覚のこと． そんな曖昧で抽象的な感覚の優劣を数値化できる方法がある．
感覚の数値化と聞いてすぐに思いつくのは5段階評価． 例えば，10種類のビールの「キレ」の良さを評価したいなら「このビールはキレが良い」のような説明文を設けて， その文に同意できるか否かを次の1～5の選択肢から回答するというもの．
この数値をリッカート尺度という．こうすると複数の異なる説明文に対して同じ基準で感覚を数値化できる． 多数の被験者の回答を平均すれば1.0～5.0の実数値で各ビールの「キレ度」が得られる． ただ，尺度を求めるべく実際に被験者にアンケートしようとすると結構難しい． というのも一本目を試飲した段階でキレが良いか回答するのは難しく，他のものを飲んで比べてみたくなる． 何らかの基準となるビールが用意されているか，全てのビールの「キレ」のおおよその分布を知らない段階で，最初に試したビールがどのくらい「キレ」が良いのか判断しずらい． というのも，こうした感覚量は相対的なものだから． でも，2つビールを比べてどちらがよりキレが良いかなら答えやすい． 全体に対して今試飲しているビールがどの位置にあるのかなどとは考えなくても良いから．
そんなときに使える方法がサーストンの一対比較． 一対比較法は「AよりBの方が良い」などと優劣を答えるだけでAのキレ度は1.5などの感覚量を数値化できる手法． ただし，得られる数値の原点（0をどこにするか）や単位には意味がない（どのようにも設定できる）ので注意が必要． 意味があるのは「AよりBの方が良い」とか「AとBの差よりもBとCの差の方が2倍ある」など相対的なもの．
サーストンの一対比較法にも欠点がある．致命的な欠点は，全員が同じ回答する程2つの優劣が明らかな場合計算できないこと．ちょっとした対処法で計算できなくなる事態は回避できるけど理論的な正しさはなくなってしまう．
もう一つの欠点も結構致命的．比べるものが多くなると被験者の回答時間が急激に増加してしまうってこと． 例えば，ビール10種類の場合は，$_{10}C_{2}=45$回も2杯のビールを飲まないといけない． リッカート尺度を使う場合は10杯で済む．
ところで，なぜ2つを比較するだけで数値化できるのか？不思議じゃないですか？ それを理解するためには理論を理解しないといけない． 幸い，この手法の理論は，大学の教養で習う統計の知識があれば理解できる． このシリーズでは「統計」も忘れちゃってる（授業中寝てた）人にも分かるように説明するつもり． そのために，まず「サーストンの一対比較法でのデータ収集方法と解析方法」で，どうやってサーストンの一対比較を使うのかデータの集め方や解析方法を説明してみる． その後，「サーストンの一対比較法の理論を理解する」で，なんでその方法で数値化できるのか理論の部分を説明してみる．
クライミングは危険を伴うスポーツです．殆どの事故は広く知られた安全な手順を間違えたり忘れていたり伝えられていなかったりすることで起きています． 特に，初心者に伝えるべき事は多いので，経験者が伝え忘れに注意しなければなりません． 最近まで初心者だった自分が初心者の気持ちを忘れないうちに伝えるべき事，自分が疑問に思ったことをメモしています． チェックリストは伝えるべき順に配置しています． 特に始めて一緒に登るけれど，何度か経験があるという人が要注意です． このチェックリストにより事故が減ることを願います．
光学透過型ディスプレイで提示したCGの仮想物体をUnity内で再現するには，単にオブジェクトを半透明するだけではダメで， 背景像に仮想物体の輝度が加算されないといけない． それを簡単に再現するUnity用シェーダを作ったので忘れないように重要な部分だけメモしておく．
Unity上では全部CGだけど，不透明なオブジェクトを実物体，光学透過型ディスプレイで提示されたつもりの半透明なオブジェクトを仮想物体と呼んでおく．
今回紹介する方法では，上の特徴全部を再現できる．
このシェーダを光学透過型ディスプレイで表示されたようにしたい物体にアタッチするだけ． このシェーダはUnlitシェーダを数行改変しただけ．
Shader "Unlit/OstDisplay"
{
Properties
{
_MainTex ("Texture", 2D) = "black" {}
_Luminance ("Luminance", Range(0.0, 1.0)) = 1.0
}
SubShader
{
Cull Back ZWrite On ZTest Less
Blend One One
Tags {
"Queue"="AlphaTest"
"IgnoreProjector"="True"
"RenderType"="Transparent"
"PreviewType"="Plane"
}
Pass
{
CGPROGRAM
#pragma vertex vert
#pragma fragment frag
#include "UnityCG.cginc"
struct appdata
{
float4 vertex : POSITION;
float2 uv : TEXCOORD0;
};
struct v2f
{
float2 uv : TEXCOORD0;
float4 vertex : SV_POSITION;
};
v2f vert (appdata v)
{
v2f o;
o.vertex = UnityObjectToClipPos(v.vertex);
o.uv = v.uv;
return o;
}
sampler2D _MainTex;
float _Luminance;
fixed4 frag (v2f i) : SV_Target
{
fixed4 col = tex2D(_MainTex, i.uv);
return col*_Luminance;
}
ENDCG
}
}
}
左が実際に使用したテクスチャで右がレンダリング結果． 結果をみると分かるようにチェスボード模様の黒い部分は完全に透明になっていて背景をそのまま透過させている（特徴1）． 黒じゃない部分は背景色に白や緑が加算されている（特徴1）． 青い球と後ろにある仮想物体を見てみるとちゃんと仮想物体が実物体に遮蔽されている（特徴2）． 一方，仮想物体どうしの重なりを見ると後ろの仮想物体の「あ」の一部は描画されていないので仮想物体どうしでも遮蔽が実現できている（特徴3）． 一つの仮想物体だけに注目すると立方体の3面しか見えていないので単一の仮想物体の背面が前面に遮蔽されているのも分かる（特徴4）． ということで必要な全ての特徴を満たすことができた． だからこのシェーダは提示する仮想物体の幾何学的整合性が全部満たされたときに，光学透過型ディスプレイでどのように見えるかを再現していると言える．
まずこの部分．この部分はディスプレイに表示するテクスチャ（_MainTex）とディスプレイの輝度（_Luminance）の初期化に関する部分． Properties内に変数と初期値を定義するとインスペクタ上で選択できる．初期値はデフォルトの値となる．
Properties
{
_MainTex ("Texture", 2D) = "black" {}
_Luminance ("Luminance", Range(0.0, 1.0)) = 1.0
}
インスペクタ上で設定した値はこの変数に代入される． これを定義していないとコンパイルエラーになる．
sampler2D _MainTex;
float _Luminance;
次にこの部分．透明な物体を描く上で遮蔽関係は重要．ディスプレイにどのように3D CGを表示したいかによるけど， ここでは3Dオブジェクトでも一旦は普通に描画してその結果が光学透過型ディスプレイに表示されている，という感じにしたい． だから，背面は普通にカリングされて表示されないようにした上で，後ろの物体とブレンドされないといけえない．
Cull Back ZWrite On ZTest Less
Cull Backはカリング（Cull）される対象が物体の背面（Back）という意味．左がCull Back，中央がCull Off，右がCull Frontの結果． 各オブジェクト内のポリゴンで背面を描かないという設定で特徴4に対応．
左：Cull Back，中央：Cull Off．右：Cull Front．
ZWrite Onはこの仮想物体のデプスをzバッファに書き込むかどうかを指定している． 左がZWrite On，右がOff．透明な仮想物体（QueueがAlphaTestやTransparent）は，不透明な物体の描画が終わってから描かれる． だから仮想物体自体のデプスが影響するのは別の仮想物体．右図では，後ろの仮想物体のZバッファないため，どちらの仮想物体を描くときも実物体の奥行きと仮想物体の奥行きが比較されて，実物体より手前にある部分が描画されている．そのため仮想物体どうし重なる部分の画素が2重に描かれてしまっている．
左：ZWrite On，右：ZWrite Off．
ZTest Lessは，既に描かれたzバッファに対してこの仮想物体のポリゴンの奥行きが手前（Less）にあるときのみ描画することを意味している．下図の通りZTest Less（左）のときは正常に描画されているけど，ZTest Greater（右）のときは描画された実物体のzバッファの値よりも奥（Greater）にある透明物体のポリゴンだけが描かれる．
左：ZTest Less，右：ZTest Greater．
透明な物体を描くということは前景と背景が混ぜるということなので，ここが一番重要． 光学透過型ディスプレイでは，実シーンからの光に，ディスプレイの光が加算される．
Blend One One
左のOneが手前の物体のアルファ値で，後ろのOneが背景のアルファ値と考えればいい． Unityのマニュアルでは左がsource，右がdestinationとなっているけどそのネーミングが分かりにくい． 数式で表せば分かりやすい．透明な物体は，不透明な物体が全部描き終えてから，一番奥の透明物体から手前の透明物体へ順に描画される． $i$番目の透明物体（もしくは最後の不透明物体）が描かれた時点での画素値を$C_{i}$，今から描画しようとする透明物体の画素値を$C_{\rm obj}とすると，その物体の描画結果$C_{i+i}$は次式で表すことができる．
\[ C_{i+1} = w_{\rm obj} C_{\rm obj} + w_{i} C_{i} \]$w_{\rm obj}$と$w_{i}$は重みでこれが左のOne（$w_{\rm obj}=1$）と右のOne（$w_{i}=1$）に対応している． この結果$C_{i+1}$は次の物体を描画する際には右辺に入ることになる． この重みを変えてみた結果が下の図．
左：Blend One One．右：Blend Zero Zero．
左：Blend Zero One．右：Blend One Zero．
Tagはシェーダの動作を切り替える役割を持っている． デフォルトのUnlitシェーダからの変更点は，1行目と3行目．
Tags {
"Queue"="AlphaTest"
"IgnoreProjector"="True"
"RenderType"="Transparent"
"PreviewType"="Plane"
}
Queueは，描画順序を決めるタグ． 透明物体に関するタグはAlphaTestとTransparentがあり，不透明物体はGeometry．他にもあるけどこの3つなら Geometry→Transparent→AlphaTestの順に描画される．
次の図はこの3つのタグを試してみた結果． AlphaTestとGeometryが全く同じ結果に見えるし，TransparentはZWriteをOffにしたときの結果と結果と全く同じに見える． 描画順序だけの問題ならこうはならないように思えるのでこのタグによって何が変化するのか理解できていないかも． マニュアルを読んだだけでは理解できなかった．
左：AlphaTest．中央：Transparent．右：Geometry．
RenderTypeは，特定のオブジェクトだけシェーダを切り替えるためのタグ． これはOpaqueに変えようがTransparentに変えようが結果に変化なかった．今回は関係なさそう．
最後は仮想物体の輝度を反映させる部分． インスペクターで設定した輝度（Luminance）を描画結果の輝度に掛け算する．
fixed4 col = tex2D(_MainTex, i.uv);
return col*_Luminance;
沢登りをしていて滝の上下に分かれたときにコミュニケーション手段がなくて困ることが良くある． 例えば，リードで滝を登る人（リーダ）とそれに続く人（フォロア），ビレイをする人（ビレイヤ）との間の会話を
みたいなやりとりをする．このやりとりはフリークライミングなどでは成り立っても，滝の音がごうごうと鳴り響く沢では，殆どまともにできないことが多い． これらのコールの多くはそのタイミングが伝わればいいだけなので，作業が順調なら笛で代用できるけれども，ちょっとしたトラブルがあるともう何の笛なのか分からなくなる．
そこで沢登り用に無線機を持っていけば良いと考えるけど沢登りの特殊な環境を考えると要求される機能がある． デジタル簡易無線なら出力も大きく免許不要でクリアな音になるということで比較してみたい．
まず当たり前だけど防水性能が高くないといけない． 防水袋に入れておけば良いと考えるかもしれないけど，袋の上からでは操作しにくく袋から出したりすると濡れた手で触らないといけないかもしれない． それに防水袋でも浸水することは良くあること．
普通無線は自分が発話するときだけボタンを押すPPT (Push to Talk)という使い方をする． ただ，滝を登っている最中にボタンを押すことは難しい場合もあるので，発話しているときだけボタンを押したことにしてくれるVOX (Voice Operated Transmit)の機能があると良い．実際に滝の音がごうごうと鳴り響く状況下で機能するかどうかは不明．ノイズキャンセル機能のあるインカムも あるのでちょっと期待．
無線機本体は厳重に防水したいのでBluetoothでインカムと本体とを分離できたら浸水したときにも被害がインカムだけになる． マイク付きのイヤホンは無線用途に限らずスマホでの会話用にも売られているので比較的安いというのもある．
市販のBluetoothマイクが完全防水ならいいんだけど，滝の水しぶきを浴びながらVOXが機能するのかあやしい． そこで有線マイクも装備しておきたいのでマイクも防水であるべき． 本体は防水しておいてマイクだけ外にだして別途防水して使う．もちろん水が侵入する可能性があるのでもともと防水であることが望ましい．
代表的な無線機メーカーは，ICOM，ALINCO，KENWOOD，八重洲無線． この中から上の条件を満たすものをピックアップして比較してみる． 重量は全て付属のバッテリを装着したときのもので，使用時間は各社のバッテーリーセーブ機能ONの時の値． 全て取り扱い説明書を確認の上，Bluetooth接続時にVOX機能が使えることを確認済み．
| 重量 | 使用時間 | 防塵防水 | マイク | |
|---|---|---|---|---|
| ICOM IC-DPR7SBT | 245g | 13時間 | IP67 | HM-159FS/SJ*** |
| ALINCO DJ-DPS71* | 244g | 15時間 | IP68 | EME-80BMA87B (IP67)** |
| KENWOOD TPZ-D563BT | 247g | 15時間 | IP68 | KMC-55 (IP67) |
| 八重洲無線 SR740 | 237g | 16時間 | IP68 | SSM-10C (IP57) |
| 八重洲無線 SR741 | 248g | 16時間 | IP68 | SSM-10C (IP57) |
* 取り扱い説明書にBluetooth使用時にVOX機能が動作するか，明確な記載がなかった．「VOX 動作レベル」の設定の説明項目にBluetooth非対応のマークが無いのでBluetooth対応と読み取れる．
** Bluetooth接続（他は有線）
***防水性能は本体に準じる
結局どれがいいのだろう？ 現状，沢登りに十分耐えうる性能を取説で保証している製品はない． 市販BluetoothヘッドセットでもVOX機能がノイズに反応せずに機能するのかはかなり怪しい． そうかと思って純正マイクを取り付けようとしてもマイクが本体よりも防水性能が低いものばかり．なんで？ 純正のBluetoothマイクヘッドホンはそもそも防水でもない． どうして無線という極めてアウトドア向きな製品なのに防水性能を軽視するのだろう？ 重さや動作時間は僅かな差なので，この表ではKENWOOD TPZ-D563BTが一番ましなのかもしれない．
この式で表される多次元ガウス関数をscikit-learnなど使わずに自分で実装したい． 下の通りだけど実装するにあたりWeb上のいろんな解説を見ていると ガウス関数のPythonコードみたいに何故かexp関数に代入する前にdiag関数を呼んでいる． なんでdiag関数が必要なのか？これが本当にいい方法なのか？実はこのdiag関数を使うと異常にメモリを消費して計算速度が遅くなる． 今回はそれらを確認してみた．
\begin{equation} \newcommand{\vec}{\boldsymbol} f(\vec{x}) = \frac{1}{\sqrt{(2\pi)^n \left|\Sigma\right|}}\exp \left \{-\frac{1}{2}{} \vec{d}^T\,{\Sigma}^{-1} \vec{d} \right\},\ \ \vec{d}=\vec{x}-\vec{\mu} \end{equation}
import numpy as np
def gaussian(x,mu,sigma):
det = np.linalg.det(sigma)
inv = np.linalg.inv(sigma)
d = x - mu
n = x.ndim
norm = (np.sqrt((2 * np.pi) ** n * det))
return np.exp(-np.diag(d@inv@d.T)/2.0) / norm
そもそも何が疑問なのかというと，exp関数の入力はスカラーなのだからnp.diagなんて要らないのでは？というもの． diag関数は行列が入力されたときに対角成分を出力する関数で例えば次のような行列Aに対してはベクトルが出力になる．
A = np.array([
[1, 2],
[3, 4]])
print(np.diag(A))
[1 4]
そこでdiagにどのような次元の行列が入力されることになるのかを見てみる． 多次元ガウス関数は入力が多次元なのであって，出力は1次元． なのでexp関数は，いわゆる高校数学で習う指数関数なので入力も出力も1次元． 指数部分を見ると
diff@inv@diff.T
は数式では$\vec{d}^T\,{\Sigma}^{-1} \vec{d}$． $\vec{d}=\vec{x}-\vec{\mu}$は$n\times 1$のベクトルどうしの引き算なのでその結果も$n\times 1$のベクトル． ${\Sigma}^{-1}$は分散共分散行列であり$n\times n$の対称行列． $n\times 1$ベクトルと$n\times n$行列の積は$n\times 1$ベクトル．つまり${\Sigma}^{-1}\vec{d}$は$n\times 1$ベクトル． それと$\vec{d}^T$の積は，内積そのものだから$\vec{d}^T\,{\Sigma}^{-1} \vec{d}$はスカラーになる． つまり，行列を入力とするdiag関数なんて不要じゃねーか？というのが疑問．
実際にdiag関数はスカラーでは次のようなエラーがでる．
np.diag(100)
ValueError: Input must be 1- or 2-d.
そのため上のガウス関数gaussianはデータ点1つでは正しく動作しない． やはりdiag関数にスカラーを入力したときと同じエラーがでる．
mu = np.array([10,10])
sigma = np.array([
[4**2, 2],
[2, 3**2]])
x = np.array([1,1])
gaussian(x, mu, sigma)
ValueError: Input must be 1- or 2-d.
要するにgaussian関数は，入力のベクトルが多次元というだけでなく， そのベクトルが多数ある場合のbroadcastを前提とした設計だった． そこで2次元んベクトルが3つの場合を確認してみる．
x = np.array([
[1,1],
[2,2],
[3,3])
diff@inv@dif.T
[[12.15 10.8 9.45] [10.8 9.6 8.4 ] [ 9.45 8.4 7.35]]
なんと$3\times 3$の行列になってしまう． gaussian関数の入力が2次元ベクトルのとき結果がスカラーなら， 入力が2次元bベクトル3個分の$3\times 2$行列なら結果は，スカラー3個分の$3\times 1$などになりそうなものだけどそうではなかった．
一体どうなっているか考えてみた．まず$\vec{d}^T\Sigma^{-1}\vec{d}$の部分を成分表示してみる．
\begin{eqnarray} \vec{d}^T\Sigma^{-1}\vec{d} &=& \left[ \begin{array}{cc} d_{11} & d_{12}\\ \end{array}\right] \left[ \begin{array}{ccc} \sigma_{11} & \sigma_{21}\\ \sigma_{21} & \sigma_{22}\\ \end{array}\right] \left[ \begin{array}{c} d_{11} \\ d_{12} \\ \end{array}\right] \end{eqnarray}
変数の中身はともかくこんな形になっているはず．データ点が3つあれば変数d（つまり数式上の$d^T$）は3行になるので↓のようになるはず．
\begin{eqnarray} \left[ \begin{array}{cc} d_{11} & d_{12}\\ d_{21} & d_{22}\\ d_{31} & d_{32}\\ \end{array}\right] \left[ \begin{array}{cc} \sigma_{11} & \sigma_{21}\\ \sigma_{21} & \sigma_{22}\\ \end{array}\right] \left[ \begin{array}{c} d_{11} & d_{21} & d_{31}\\ d_{12} & d_{22} & d_{32}\\ \end{array}\right] \end{eqnarray}この行列のまま計算すると数式が複雑になりすぎるので成分のみの計算にする． 左2つのベクトルの積を$\vec{b}=d^T\Sigma^{-1}$とすると，その成分は
\begin{eqnarray} b_{ij} &=&\sum_l d_{il}\,\sigma_{lj} \end{eqnarray}と表せる．これは行列の積の公式そのもの．これを用いて残りも含めてた全体を成分で表示してみる． $\vec{a}=d^T\Sigma^{-1}\vec{d}=\vec{b}\vec{d}$とすると，
\begin{eqnarray} a_{ij} &=&\sum_k b_{ik}\,d_{jk}\\ &=&\sum_k \sum_l d_{il}\,\sigma_{lk}\,d_{jk}\label{tensor}\\ \end{eqnarray}と表せる．この成分をよく見ると$d_{ik}$の$i$がデータ数に対応する． 最後の式で$d_{il}$と$d_{jk}$とあるけど，もしデータが1つしか無かったら$i=1$の一行だけになりこのとき$d_{1x}$の形の成分しか存在しない． ということは$i=j=1$に限られることになる．2つ目のデータしかないときは$i=j=2$に限られ，以下同様． 結局$a_{ij}$の対角成分が$d^T\Sigma^{-1}\vec{d}$の結果を表していることになる．
\begin{eqnarray} \left[ \begin{array}{ccc} a_{11} & & \\ & a_{22} & \\ & & a_{33}\\ \end{array}\right] \end{eqnarray}だからdiag関数が必要だった．ブロードキャストが原因だった．
ここまで見てすぐに分かるのは行列$a_{ij}$の対角成分以外は必要ない． なのにデータが$N$個ある場合，一旦$N\times N$行列が生成されていることになる． これはデータ数が多いとすぐに破綻するメモリ的にも計算時間的にも無駄が多い．
def gaussian2(x, mu, sigma):
det = np.linalg.det(sigma)
inv = np.linalg.inv(sigma)
d = x - mu
n = x.ndim
norm = np.sqrt((2 * np.pi) ** n * det)
power = - np.einsum('il,lk,ik->i', d, inv, d) /2.0
return np.exp(power) / norm
einsum関数は，添え字を指定するだけでの行列（テンソル）どおしの必要な積が定義できる関数． 式$\ref{tensor}$の添え字をそのまま記入すればいい．$i=j$なので$j$は$i$に置き換えられている．
x = np.linspace(0, 40,40)
y = np.linspace(0, 40,40)
X, Y = np.meshgrid(x, y)
XY = np.c_[X.ravel(), Y.ravel()]
mu = np.array([10,10])
sigma = np.array([
[4**2, 2],
[2, 3**2]])
これらのパラメータでgaussianとgaussian2の計算時間をテストしてみた．
timeit gaussian(XY, mu, sigma)
16.5 ms ± 207 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
timeit gaussian2(XY, mu, sigma)
198 µs ± 1.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
0.198/16.5=0.012とおよそ100倍のスピード． loopの回数は自動的に決まるみたい． とにかく圧倒的にdiag関数は無駄無駄無駄．
共分散行列に全部値が詰まった一般の二次元ガウス関数の最小二乗法の実装メモ． scikit-learnのsklearn.mixture.GaussianMixtureを用いてEMアルゴリズムにより混合ガウス分布を推定するプログラムのサンプルはあちこちで見かけるけど最小二乗法が意外と無い． EMアルゴリズムと最小二乗法の違いは，EMアルゴリズムはデータ点$(x_i, y_i)$からガウス分布のパラメータを推定するのに対して，最小二乗法はデータ点$(x_i, y_i, z_i)$からガウス分布のパラメータを推定する部分． 最小二乗法の関数scipy.optimize.leastsqの使い方とか独特なのでちょっと手こずったから忘れないように．
まずは肝心のガウス関数． ベクトル・行列の表記だけを見ると二次元でも多次元でも同じなので，とりあえず多次元版を載せておく．
\begin{equation} \newcommand{\vec}{\boldsymbol} f(\vec{x}) = \frac{1}{\sqrt{(2\pi)^n \left|\Sigma\right|}}\exp \left \{-\frac{1}{2}{} (\vec{x}-\vec{\mu})^T\,{\Sigma}^{-1} (\vec{x}-\vec{\mu}) \right\} \end{equation}でもやっぱり成分表記しないとイメージしにくいので二次元版（$n=2$）のために成分も載せておく， ちなみに$\sigma_{xy}=0$のときは，$x$軸か$y$軸方向にガウス関数が広がりを持っていて，$\sigma_{xy}\neq 0$の時はそれが回転したような形になっている．
\begin{equation} \vec{\mu}= \left[ \begin{array}{c} \mu_x \\ \mu_y \end{array} \right] ,\ \ \ \Sigma = \left[ \begin{array}{cc} \sigma_x^2 & \sigma_{xy} \\ \sigma_{xy} & \sigma_y^2 \end{array}\right] \end{equation}これをPythonコード化したものがこれ↓↓． leastsqをが最小二乗法の関数． 「＠」は数学で習う行列とベクトルの積．横ベクトルと縦ベクトルを内積とって左辺の要素を計算するあれ． A @ Bを別の表記にすればnp.dot(A, B)です．
import numpy as np import matplotlib.pyplot as plt from scipy.optimize import leastsq
def gaussian(x,mu,sigma):
det = np.linalg.det(sigma)
inv = np.linalg.inv(sigma)
n = x.ndim
return np.exp(-np.diag((x - mu)@inv@(x - mu).T)/2.0) / (np.sqrt((2 * np.pi) ** n * det))
この関数が正しく実装できているか確認するために可視化する． 可視化するためにmeshgridによって$(x,y)$の密な点配列を生成する． XYは$40\times40=1600\,$個の二次元座標列を$1600\times2$の行列として成形されている．
x = np.linspace(0, 40,40) y = np.linspace(0, 40,40) X, Y = np.meshgrid(x, y) XY = np.c_[X.ravel(), Y.ravel()]
これが関数本体の計算部分． 最初に適当なパラメータ$\vec{\mu}$や$\Sigma$の値を設定している． gaussian関数の戻り値Zは，$1600\times1$のベクトルなので，描画用に元の$40\times40$の行列に戻している．
# パラメータ設定
mu = np.array([25,15])
sigma_x = 8.0
sigma_xy = 3.0
sigma_y = 3.0
sigma = np.array([
[sigma_x**2, sigma_xy],
[sigma_xy, sigma_y**2]])
# ガウス関数の計算
Z = gaussian(XY, mu, sigma)
z = Z.reshape(X.shape)
# ガウス関数の計算
max = np.max(Z)
min = np.min(Z)
print("min: %f, max: %f" % (min, max))
最大値，最小値の結果は次のとおり．
min: 0.000000, max: 0.006615
ここは普通のグラフの描画と同じなので特に説明の必要はないかも． x軸，y軸の範囲を[0, 39]しているのは40を最大にすると空白の1画素ができてしまったから．
fig = plt.figure(figsize = (6, 6))
plt.xlim([0, 39])
plt.ylim([0, 39])
plt.imshow(z, cmap="jet", vmin=min, vmax=max)
plt.subplots_adjust(left=0.04, right=1, bottom=0.04, top=1)
plt.savefig("gauss.png")
plt.show()
データ点は乱数で生成してみる．素なグリッドでも良い． $(x,y)$だけ乱数でテータ点$(x_i,y_i)$を生成して，各点に対して$f(x_i,y_i)$を計算してそこに正規分布の誤差を加えてみる．正規分布の分散はかなり小さい値にしないといけない．そもそも，二次元ガウス関数の最大値もかなり小さいので．
n = 1024*768//2**8 x_rand = np.random.random_sample(n)*40 y_rand = np.random.random_sample(n)*40 xy_rand = np.c_[x_rand.ravel(), y_rand.ravel()] z_rand = gaussian(xy_rand, mu, sigma)+ np.random.normal(0.0, 0.0001, xy_rand.shape[0])
ここが可視化部分．分かりやすいように上の画像をα値0.5の半透明な背景としている．
fig = plt.figure(figsize = (6, 6))
plt.xlim([0, 39])
plt.ylim([0, 39])
plt.imshow(z, cmap="jet", alpha=0.5)
plt.scatter(x_rand, y_rand, c = z_rand, cmap="jet", edgecolors="black")
plt.subplots_adjust(left=0.04, right=1, bottom=0.04, top=1)
plt.savefig("gauss_data.png")
最小二乗法には評価関数を定義しないといけない． 評価関数はデータ$z_i$とデータ点$(x_i,y_i)$でのガウス関数の値$f(x_i, y_i)$の差の二乗和． ただし，leastsqではLevenberg-Marquardt法が使われているので与える必要があるのは$\epsilon_i = z_i - f(x_i, y_i)$の部分だけ． これがn個あるのでxdataにはn$\times$2の行列が，zdataにはn$\times$1の行列が与えられるはず． paramsは5次元ベクトルで，前の2要素に[$\mu_x, \mu_y$]，後の3要素[$\sigma_x, \sigma_xy, \sigma_y$]が格納されているとする．
#評価関数定義
def error(params, xdata, zdata):
mu = params[0:2]
sigma = np.array([
[params[2], params[3]],
[params[3], params[4]]
])
return zdata - gaussian(xdata, mu, sigma)
最小二乗法のleastsq関数に与えるのは評価関数errorだけでなく初期パラメータx0，データargsを与えないといけない． 初期パラメータは上でテストしたときの値（真値）に適当な誤差2.0と1.5を加えて代入． データはn$\times$2のxy_randとn$\times$1のz_randを与えてそれぞれxdata，zdataに代入させる．
# パラメータを1次元ベクトル化 mu_init = mu+2 sigma_init = sigma.reshape(4)+1.5 params = np.array(empty(5) params[0:2] = mu_init params[2:] = np.array([sigma[0,0], sigma[0, 1], sigma[1, 1]]) # 最小二乗法 results = leastsq(error, x0=params, args = (xy_rand, z_rand)) print(results[0])
これが結果．
[25.04568978 14.93655215 61.81134792 -1.25051091 5.7730626 9.1594342 ]
分かりにくいので$\vec{\mu}$と$\sigma_x, \sigma_y$に分けて表示してみる．
mu_est = params_est[0:2]
sigma_est = params_est[2:].reshape(2,2)
print("mu:")
print(mu_est)
print("sigma:")
print(sigma_est)
print("sigma_x: %f" % np.sqrt(sigma_est[0,0]))
print("sigma_y; %f" % np.sqrt(sigma_est[1,1]))
mu: [25.11389534 14.9731751 ] sigma: [[66.22148794 3.23015128] [ 3.23015128 9.05173089]] sigma_x: 8.137659 sigma_y: 3.008609 sigma_xy: 3.230151
真値は$\vec{\mu}=[25, 15]^T$なの中心位置は大体あってる． $(\sigma_x, \sigma_{xy}, \sigma_y)=(8.0, 3.0, 3.0)$についてもおおよそ正しい．
ここからは最初の繰り返し． 求まったパラメータを使って再度ガウス関数を描画するだけ． XYだけは使い回せる．ぱっと見近い結果だし問題ないかな．
param_est = results[0]
# 推定されたパラメータをセットする
mu_est = params_est[0:2].reshape(2)
sigma_est = params_est[2:6].reshape(2,2)
# ガウス関数計算
Z_est = gaussian(XY, mu_est, sigma_est)
z_est = Z_est.reshape(X.shape)
# 描画
fig = plt.figure(figsize = (6, 6))
plt.xlim([0, 39])
plt.ylim([0, 39])
plt.imshow(z_est, cmap="jet", vmin=min, vmax=max)
plt.subplots_adjust(left=0.04, right=1, bottom=0.04, top=1)
plt.savefig("gauss_est.png")
plt.show()
真の分布zと推定されたガウス分布z_estとの差を調べてみる． これには単純に二つの差の画像を作って二乗和を計算してみる．
E = np.abs(z_est - z) fig = plt.figure(figsize = (6, 6)) plt.xlim([0, 39]) plt.ylim([0, 39]) plt.subplots_adjust(left=0.04, right=1, bottom=0.04, top=1) plt.imshow(E, cmap="jet")
誤差の分布が変な形になっているけど，平均値は1.109483e-05であり，加えた誤差の標準偏差0.0001=(1.0e-04)よりも小さい．なので上手くいったということにしておく．
e = np.sum(E)/E.size print(e)
1.1094833018535245e-05
二次元の確率密度分布$f(x,y)$が与えられているとする．三次元空間$(x,y,z)$内に一様分布の乱数点$(x_i,y_i,z_i)$を生成し，次式を満たす箇所のみ点$(x_i,y_i)$を残すことで目的の乱数が得られる．一様分布を生成する$z$の範囲は$f(x,y)$の最大値$z_{\rm max}$を用いて$0 \leq z < z_{\rm max}$とすれば，$0\leq z < 1$としたときと等価な乱数が得られる．
\begin{eqnarray} z_i < f(x_i,y_i) \end{eqnarray}必要なパッケージのimportとよく使う関数を定義しておく． plt.imshowとplt.savefigで画像を保存しても良かったけど画像サイズが変わってしまうのでcv2.imwriteを使うことにした．
import cv2
import numpy as np
def jet(img):
img_norm = np.zeros_like(img)
cv2.normalize(img, img_norm, 0, 255, cv2.NORM_MINMAX);
img_jet = cv2.applyColorMap(img_norm.astype(np.uint8), cv2.COLORMAP_JET)
return img_jet
まずは，二次元の確率密度分布$f(x,y)$を画像として表す． 単純な発想では顕著性マップは注視位置の確率分布と考えられるので，適当な画像から顕著性マップ作ってみる． ここでは結果だけ載せておくけど使用したのはDeepGaze II．入力画像はDeepGaze IIの入力画像サイズに合わせて1024$\times$576を1024$\times$768にしてあるので少し縦に伸びている． ここでは入力画像も顕著性マップも何でもいいのでこれで良しとする．
sal = np.load("keiba_sal.npy")
norm = np.linalg.norm(sal, axis=None)
sal = sal/norm
cv2.imwrite("keiba_sal.png", jet(sal))
$(x,y,z)$の各変数で独立に一様乱数を発生させれば$(x,y,z)$の一様乱数ができるけど，$(x,y)$に関しては画像の座標なので整数値に限定した． ここでは10000個の乱数点を発生させる．といっても，x,y,zにそれぞれ座標（乱数列）が代入されるだけ．
max = np.max(sal) n = 10000 x = np.random.randint(0,sal.shape[1],n) y = np.random.randint(0,sal.shape[0],n) z = np.random.rand(n)*max
これは結果の確認のためのコード無くても問題ない．画像uniは後で使わない．点かなり多め．．．
uni = np.zeros(shape = sal.shape, dtype="int8")
uni[y,x] = 1
plt.imshow(uni, cmap="gray")
plt.axis('off');
plt.subplots_adjust(left=0, right=1, bottom=0, top=1)
cv2.imwrite("uniform.png", rnd)
この記事の一番肝のコードがこれ． 1行目でrndには画像と同じサイズの中身が0の画像が作製される．int8型にしたのは後でor演算して画像に収めるため． 2行目では，$z < f(x,y) $を満たすとき1，満たさないとき0を，x,yのデータ数（ベクトルの次元数）分だけpntに代入している． なのでpntはxやyと同じshapeになる． 3行目では，画像rndの各画素$(x,y)$にpntの1か0が代入される．
rnd = np.zeros(shape = sal.shape, dtype="int8") pnt = np.where(z < sal[y,x], 1, 0) rnd[y,x] = rnd[y,x] | pnt
左が結果．なんか点が少ない．確かに$z_i < f(x_i,y_i) $を満たす点だけに絞られるので数がかなり減る．$z_{\rm max}$付近の値をとるのは$(x,y)$空間中のごく一部だけ．その他の点は殆ど無駄になる．これがこのアルゴリズムの難点． nを100倍のn = 10000000にして再実行すると次が得られる．それっぽい結果（右）が得られた．もっと効率の良い方法があると思う．
パイプは，プロセス間通信の一番簡単な手段だと思う． 「|」という記号でコマンドを接続するだけで異なるプログラム通しがデータをやり取りできるようになる． これが便利なのは異なる言語で書かれた既存のプログラムをできるだけ変更せずに接続したい場合に， C言語ならprintf/scanf，Pythonならprint/input文を送信側/受信側に追加するだけでプロセス間通信が実現できる． この便利なパイプがWindows上で妙な挙動をすることを発見したので記録しておく．
送信側（sender.py）のプログラムから1秒おきにインクリメントされた数字を標準出力に出力し， 受信側（receiver.py）ではinput関数により受け取った標準入力を即座に標準出力に出す． この2つのプログラムを次のようにパイプでつなぐ．
> python sender.py | python receiver.py
import sys
import time
for i in range(5):
print(i, flush=True)
print("sender: "+ str(i), file=sys.stderr, flush=True)
time.sleep(1)
for i in range(5):
a = input()
print("reciever: "+ a)
結果は，どちらもパイプを通してデータの受け渡しが出来ているけど，そのタイミングが全然違う． コマンドプロンプトの場合は送信側が送信し終えたらすぐにパイプを伝って受信側に行き即座に受信側から出力されているけど， PowerShellの場合は，送信側が終了するまで受信側が動作していない．
sender: 0 reciever: 0 sender: 1 reciever: 1 sender: 2 reciever: 2 sender: 3 reciever: 3 sender: 4 reciever: 4
sender: 0 sender: 1 sender: 2 sender: 3 sender: 4 reciever: 0 reciever: 1 reciever: 2 reciever: 3 reciever: 4
要するにPowerShell上だと送信側が終了しないと受信側にデータが受け渡されない． これってコマンドプロンプトからPowerShellになって退化したってことでしょうか？ 不具合ではなくそういう仕様？今更PowerShellを勉強する気もないのでそれ以上調べずにここで終わることにした． そういう仕様だとしてもややこしいので同じ記号にするなよって思う． 変なシェルを作らずにUnix系のシェルで統一して欲しい．
バックカントリー（BC）用のザック（バックパック）どれにしようかと，散々迷って比較用のページまで作った挙げ句に，一番格好いいNorrona lofoten 30Lを購入した． 迷った理由は，格好いいのはいいけど，重さと外観くらいでどんな機能があるザックなのか，写真も動画もいまいち出回っていなくて判断つかなかったから． 結論から言うと買って良かった！ 最近のザックらしくいろんな便利機能が搭載されたザックだった． 他にも迷う人のためにザックの機能を詳しく紹介したい．
外観はこんな感じ．最近流行りのストラップが全部隠れてツルンとした感じでも，防水袋みたいな感じでもない，いわゆる昔ながらのザックっぽいザック． ミリタリーっぽくて格好いいけど，ベルトだらけの重そうな感じでもない．持った感じ最近のザックらしくめちゃくちゃ軽い． 逆に耐久性大丈夫かな？と思ってしまう．横向きの二本のベルトがあるからスノーボード固定できるのは明らか．
Norronaやlofotenなどのロゴを見てみよう．正面の上にはnorronaの「おっさんの顔」マークがあって，下には「lofoten 30」の文字． 背中側にも腰あたりに「おっさんの顔」マークがある．右肩のベルトにもNorronaの文字が．文字やマークが控えめな色で格好いい．
では最も大きな背中側のスペースを見ていこう． ファスナーがショルダーストラップの内側にあるので，ショルダーストラップを外さずに背面を開けることができる背面エントリー仕様は，最近のOspreyのBC用ザックと同じ． ただ形が崩れないように正面側のポケットに毛布を入れているので白い部分がもっこり膨らんでいて，少し狭く見えるかも． じつはこのタイプの背面がガバーッと開く仕組みはあまり使いやすくないと思ってる．ぎゅーぎゅーに荷物を入れるとファスナーが閉まらなくなるし，空けたときに小物がこぼれ出てしまうから． ただ，スノーボードを正面に装着すると上から荷物を出し入れしずらくなるので，無いなら無いでまた困る．
毛布のせいで狭く見えているけど，上と下にはこれだけのスペースがある． ただ，後で紹介する前面側のスペースと仕切りで分け合ってるだけなので荷物の配分を考えないと， 背面側ばかりに荷物を入れすぎると，前面側にはあまり入らなくなると思う．
背中側の収納スペースはこれだけではない．ベロンとめくりあげた背面パッドの裏側にもポケットがある． それも2重になっていてメッシュ状のポケットと，ファスナーで開くともう少し大きなスペースがある． こんなに細かく分割されていて何に使うんだろう．それと上に脱着可能なプラスチックのクリップがついた輪がある．ハイドレーション用のバッグとか鍵ぶら下げとくのかな？
前面のスペースはスコップとかプローブ（ゾンデ）を入れるように設計されているらしい． これは他のBC用のザックでも同じ．こっちもベロンとめくった蓋側の方にファスナーで開くポケットがある． しかもここにもクリップ付きの細いストラップが見える．
ザックの上面にもファスナーがある．これをあけると底の浅いスペースが現れる． 他社のザックの解説をみてるとここにゴーグルを入れるのが定番らしい． 山を登っているときはサングラスにして，滑るときはゴーグルに変える．さっと出し入することを要求されるような小物を入れるんだろうな．
もう一つ大事な点は，このポケットの上部がマジックテーブで開くようになっていて，大収納スペースに繋がっていること． これなんのためだろう？マジックテーブで開く穴は小さくて大収納スペースと一体化する訳ではない． 恐らくここから手を入れて大収納スペースにある小物を取り出すために作られたんじゃないかな？ これは自分の理想に近い！理想のBC用ザックは背面側と上側の両方から大収納スペースにアクセスできること． あとは巾着タイプだと言うこと無いんだけどな．
底の正面からアクセスできるここは，非常に薄くて小さなスペースで，ここにはヘルメットホルダーネットが入ってる． 無くならないように紐で本体と繋がっている．
実は，底にはもう一つスペースがある． これは側面から開く．横からものを出し入れするように作られている． ここは，ファスナーの仕切りを開ければ，大収納スペースとつながるようになっている．
最近のザックらしくストラップにもポケットがつけられている． 一つは右肩のストラップ，もう一つは腰のベルト． 腰のベルトは色々便利そうだけど肩のストラップはどういう使い道があるんだろう． 本体の収納から穴でつながってるのでハイドレーション用のチューブを通すためのところってこと？ 使わないときはチューブを全部ストラップ内にしまっておくのかな？ハイドレーション使わないのでよく分からないです．
まずは腰ベルトについているギアループ．Ospreyの「Stow-On-The-Goポールアタッチメント」機能のようにストックを固定するのは無理そう．カラビナとかちょっと掛ける程度かな？ ストックなら普通に後ろにつけることができる．あと，ザック底にある水抜き穴．水が入らなくても空気抜きの穴として使えそう．
次は，ストラップ．よく見ると，スノボを取り付ける横ストラップの根本がループになっている．これは左右両側にあってアックスを取り付けることができそう． それと，上下2本の横ストラップの下地は，よく見たら空洞になっている．左右が貫通している．何に使うんだろう？なぞが多い．
これは一番最初に紹介すべき機能かもしれないけど，スノボを固定する横ストラップのクリップ部分をよく見るとストラップが滑って緩んでしまわないようにロックする機能がある． 左がロック時で右がオープン時．重いスノボやスキーを固定するためなのでこのロック機能は横ストラップとスキー固定用の上ストラップにしかない． と思ったらショルダーストラップの下部にもあった．いる？
まだあります．胸のストラップのクリップ部分．変な突起があると思ったらなんと笛でした．なるほど．．．非常時にこの笛の存在を思い出すかな？．．必要以上に機能満載のザックに驚きでしたが，謎なのはどうしてこれをもっとアピールする動画とかページとか作らないのだろう．．．英語サイトですらほとんど詳細が紹介されていない．
最後にこの部分をアップ示して終わります． 緑ばかりじゃつまんないので赤いストラップがアクセントになってるというか． カレーライスに添えられた福神漬けと言ったところでしょうか．ミリタリー系のザックとの違いかな？
