『尾瀬―山小屋三代の記』（後藤允，岩波新書）を読む

　先日，研究室の有志メンバーと尾瀬に行ってきた．うちの大学は私立大学にしては珍しく付属演習を所有していて，その一つが群馬県の水上にある．８月の実習で自動撮影カメラを設置したので雪が降る前に回収する必要があった．ついでに足を延ばして尾瀬ヶ原を散策してきたのだ．

　尾瀬と言えば，日本の自然保護運動のメルクマールとなった場所だ．その固有の自然を愛した親子三代がダムや自動車道の建設に反対し続けたという物語は，三代目が保護運動の半ばで遭難死したという悲劇もあわせて，確かにドラマチックである．そして，運動の盛り上がりと成功は時代の転換を象徴するものでもあった．

　私の時代には，小学校6年生の教科書（光村図書）に「守る、みんなの尾瀬を」が掲載されていたこともあり，尾瀬や水芭蕉という言葉の響きに，今でも一種のロマンチシズムを感じてしまう．今回読んだ新書は，その元となった本である．

新書絶版なので古本で手に入れた．

　大人になって改めて読みなおすと，小学校の頃とは全然違った印象を持つ．例えば，平野一家が尾瀬とかかわるようになった始まりは，結構めちゃくだ．一代目長蔵氏は，明治3年の生まれで，青年期に若干怪しげな信仰心と燧ケ岳と尾瀬周辺の自然への愛に目覚め，尾瀬沼近くに勝手に山小屋を建てて移り住んでしまった．どうやら，生まれ育った村の居心地が悪くなったことも背景にはあったらしい．それが尾瀬保護へと至るそもそもの始まりなのだ．

　この本では美化されているが，長蔵氏が「奇人」であったのは確かだろう．寄付を募って燧ケ岳に参道を付け山頂に祠を建てたり，出来そうもない殖産事業（例えば尾瀬沼での養魚事業）を片っ端から試みたり，あるいは，尾瀬に発電計画が持ち上がった時には自ら内務大臣に中止を請願したりと，行動力に長け，頑固一徹，周囲におもねるところのまったくない人柄が目に浮かぶ．遠くから眺めている分にはいいが，身近にいると迷惑な厄介者．おそらくそんな感じだろう．

燧ケ岳を背景に

　しかし，個性豊かな彼の振る舞いとその成功も，やはり時代に根差したものだったことにも確かだ．例えば，信仰の対象も，国家神道の一流派であり，村の素朴な山岳信仰とは必ずしも相いれなかった．立ち上げた組織名が「燧嶽神社附属愛国講社」だったということにも時代を感じる．彼の燧ケ岳の参道の設営や山頂への祠の建設のための寄付集めも，人々の愛国心に訴えることによるものだった．

　一代目のほど強引さはないが，ニ代目、三代目の暮らしと保護運動も，その時々の時代の空気と無縁ではないだろう．二代目長英氏の山小屋経営が成立したのは登山ブームと無縁ではないし，三代目の長靖氏が始めた自動車道の反対運動も経済発展一辺倒の時代からの転換点と重なっていた．運動が本格化した1960年代後半と言えば，70年安保運動の最盛期だ（長靖氏も京大の学生時代は自治会で活発に活動したようだ*1）．また，水俣病がチッソ株式会社による公害と正式に認められたのも1970年だ（環境庁の設立も同じ年）。彼らの生活や尾瀬保護運動は、当時の流れに逆らうもののように見えて、実際にはその時代を独自の在り方で反映した鏡の一部だったのだ．

*1 長英氏は，尾瀬の山小屋に生まれ育ち，毎日往復20キロの道のりを歩いて中学に通ったという．それで京大に行くって相当の秀才な気がする…　

林道でリスが出迎えてくれた

　それでもやはり，強力な個性を持つ人間がいたからこそ，尾瀬という自然が護られたという事実は間違いない．それは「もし彼らがいなかったらどうなっていたか」を想像すれば明らかである．度重なる尾瀬の危機すべてに対応できることはなかっただろう。同じことは，尾瀬の自動車道建設を最終的に認めなかった環境庁長官大石武一氏についても言える．医師でもあり植物学への理解も深かった大石でなければ，強い反対を押し切った決定はなされなかったはずだ．歴史はある方向にしか変わらないが，特異な人間の存在と偶然が組み合わさることで出来上がる世界は確かに存在するのだ．

　「昔の人は立派だった」的な感想は，あまり抱く必要もない．そもそも保護運動のきっかけが奇人の個人的信条に端を発していることは書いた通りだ．しかし，今と比べて，外れ値とも言える個性が確かに存在しえたこと，その外れ値が生み出す奇特な物語があったこともまた確かなことだ．窒息しそうな今の時代、多くの人が尾瀬が紡いだ唯一無二の物語にきっと魅了されるはずだ。私にとってもそれは同じだ。小学生の時には違った尾瀬の豊かさに触れられた今回の散策は，実に有意義なひと時だった．

Megadetectorによる空打ち削除の可能性と限界

最初，このブログに書いていたけど，Noteの方が便利そうなので移行した．こちら→https://note.com/cameratrap/n/na3e99c3e9f95

今年一年を振り返る

　早いものでもう2022年が終わろうとしている。

　今年は事情があって実家のある明石には帰れず、藤沢で年を越すことになった。今日まで研究室に行って普通に仕事をしていたので、あんまり年の瀬の実感がない。

　いい機会なので今年を振り返ってみる。研究面では、正直停滞した1年だった。ここでは書きにくいが、プライベートで大変なことがあり、研究にばかり集中できる状態でもなかった。やっぱり、年と共に、（自分以外の事情も含めて）だんだん研究どころじゃなくなることが増えてくるんだろうな…。今出来ることを少しずつ着実に取り組めるようになりたい。そういうのが1番苦手なのだけど。

　おもに研究面で今年を振り返ると、

・ファースト一本。書いたのは去年だけど。

・共著は、3本。一本は京大のHさんがファーストでJ Applied Ecology。計画の段階から相談してきた研究だったので嬉しい。

・森林総研のIさん、国環研のFさんらと月一回のカメトラ会（カメラトラップ関係の情報交換会）が定着。哺乳類学会の自由集会にも発展させられた。この会は純粋に勉強になっていて楽しい。

・3月の生態学会は、統計関係の自由集会を企画。クリスマスには、それきっかけで繋がりのできた統数研のSさんの菜園にもデビューできた。

・統計本の執筆はそこそこ進んだけど、まだまだ先は長い。

・学会関係の仕事や本の編集、執筆など仕事がどんどん回ってくる。そろそろパンクしそうなので気をつけねば…

・D1のY君が、ようやく学振を当てた。まあ、今年1番のニュースはこれかな。早く論文完成させてほしいけど…

　仕事とは関係ないけど、ストレスからか突発性難聴になって、今も完全には回復していない。今年は、健康に気を遣いながら、マイペースで行きたい。

生存時間解析における「打ち切り」と「切断」

　生存時間解析には、打ち切りと切断という、少々紛らわしい用語が存在する。院生ゼミの時にちょっと関連知識を紹介した。その時の様子をツィートしたら、院生のところに「どう違うのか」と別の大学の院生から質問が来たらしい（すぐに質問するフットワークの軽さは見習うべきだけど。。。）。いい機会なので、理解している範囲でまとめてみる。

　結論から言うと、打ち切りありも切断ありも、観測データが不完全であることを指す点では同じだ。しかし、（少なくとも生存時間解析の文脈では）切断ありが観測対象の選び方の不完全さに関連するのに対し、打ち切りは観測自体の不完全さ（例えば、観測期間の短さ）に由来する。以下に例を挙げながら説明する。

1. 打ち切りについて

　打ち切りは、生態学関連でも出てくるし、RESTでも使うので私にとってが馴染み深い。

　（なんでもいいのだけど、あえてRESTっぽい状況を想定して）動物のある行動の持続時間を測りたいとしよう。行動開始から終了までを常に観察できれば理想的だが、必ずしもそうはいかない。例えば、動物を観察・記録するための装置が「ある時間が経過したら勝手に撮影を終了してしまう」仕様になっているかもしれないし、行動開始時点ではうっかりしていて気づかずに「途中からしか観察出来なかった』ということもあるかもしれない。

　これらのデータをそのまま使って行動持続時間を推定すると、過小評価してしまう。このような制約を伴ったデータのことを打ち切りデータと呼ぶ。しかし、打ち切りデータも「少なくとも〇〇より長い」ということを教えてくれるので、貴重なデータなので捨ててしまうのはもったいない。

2. 打ち切りの種類

　打ち切りには、左打ち切りと右打ち切り、両側打ち切り、そして間隔打ち切りがある。

　左打ち切りというのは、先程あげた例の2番目が該当する。すなわち、開始時点が不明な状況だ。右打ち切りとは、先ほどの1番目の例のように、終了時点が不明な状況を指す。そして、両側打ち切りとは、この両方が生じている場合を指す。

　最後の間隔打ち切りは、前提とする観測条件が少しだけ違う。例えば、動物がある行動を続けているかどうかを1分間隔で記録したとする。この場合、持続時間の最小値と最大値は分かるが、ピッタリこれだけの長さというデータを得ることは原理的にできない。このようなデータを間隔打ち切りという。野外データのほとんどは間隔打ち切りな気もするが、意外にそれとして扱われることは少ない気もする。

　確率分布を仮定した解析を行う場合、打ち切りデータの尤度は割と簡単に定式化できる。それについては、こちらでも説明しているので参考にしてほしい。リンク先でも述べている通り、理屈は簡単だが、実際の解析でRの関数を使おうとすると、打ち切り（観測なし）を1として指定するものや、観測あり（打ち切りなし）を1と表現するものがあったりして、結構面倒くさかったりする。こういうので躓くのは避けたいところだ。

3. 切断とは？

　では、切断とはどのような状況を指すのだろうか。こちらは、生態学ではあんまり出てこないが、医学の世界では頻出なようだ。

　ちょうどゼミでも使った「厨二病の発症年齢を推定する」という状況を考えてみよう。中学生を入学時点から追跡調査し、いつ発症したかを記録したとする。まあ、厨二病だと発症のタイミングが不明確だが、それが仮に分かるとしよう。また、後で図にしやすくするために、13歳の誕生日で中学に入学することにしておく。

　発症年齢を推定する上で、この調査方法には大きな制約があることは明らかだ。調査対象を中学生の厨二病患者に絞っているので、①小学生の時に発症した人と、②高校生になってから発症する人が対象から抜け落ちているからだ。

　図示すると、以下のようになる。①に該当するのデータを左切断データと呼び、②に該当するのが右切断データだ。切断という表現は、真の分布のが13歳と16歳の時点で「切れている」ことを表現したものだ。

　最初に「切断ありが観測対象の選び方の不完全さに関連するのに対し、打ち切りは観測自体の不完全さ（例えば、観測期間の短さ）に由来する」と言ったが、その意味が分かってもらえるはずだ。

4. より複雑な状況

　実際には、もっと複雑な状況が起こりうる。先ほどの厨二病の発症と同時に、治癒についても記録していたとする。そして、発症年齢、治癒年齢、治癒に要する期間を推定したいとする。この場合、様々な状況を考慮しないといけなくなる．

　全部で上の図の6通りが存在する．何を推定する場合に，どれが「打ち切り」や「切断」に相当するのかを考えてみると頭を整理できるはずだ．医学の世界では，多分，「発症のタイミングによって治癒までの期間がどうなるか」みたいな解析をする必要があるんだろうと思う．しかも，途中で治療法を変えたりしながら…これは結構大変な状況だけど，生態学では，幸いそこまで複雑な状況はまれだろう（多分）．

検出エリア内の滞在時間から移動速度を推定する

はじめに

　ちょっとしんどいことがあって，仕事にもあまり手がつかない感じなのだが，研究モチベーションの維持のためにも最近取り組んでいたことをまとめてみた．もう少し落ち着いたら後で振り返ろう．

　昔から「RESTで使う正三角形内での滞在時間から移動速度を出せないのか？」と聞かれることが多い．「やればできるだろうなあ」と思いながら，類似の問題（目的と言うより手法が類似）を，元筑波大のAさんが取り組んでいたはずなので手を出さないでいた（あの研究はどうなったんだろう？）．考えてみると，まあ，我々が使っている正三角形バージョンに限定してしまえば，とくに重なるところもないはず．ちょっとやってみることにした．

問いの確認

　取り組みたい問題は，以下の通り．まず，自動撮影カメラの前に設定した正三角形の検出エリアがある．調査の後に，動画を見ながら，動物がその検出エリア内にとどまった時間（あるいは，通過するのに要する時間．以下，滞在時間）を測定したとする（REST法を使えば，滞在時間と撮影頻度から動物の密度が推定できる）．この滞在時間には，動物の移動速度の情報が含まれていることは自明だろう．もし動物が，1. ランダムな場所から侵入しランダムな場所から脱出する，2. 検出エリア内では真っすぐ動くならば，滞在時間から動物の移動速度の推定も可能な気がする．では，具体的にはどうしたらいいのだろうか？

推定の手順

滞在時間の分解

　最初に考えなければならないのは，測定される滞在時間と速度の関係である．これは，難しい話でもなんでもなくて，\(滞在時間t = 距離y / 移動速度x\)になる．

　さて，距離yや移動速度xは，動物が通過するたびに確率的に変わる値，すなわち，確率変数である．今，移動速度xが従う確率分布をf(x)，距離yが従う確率分布をg(y)としよう．すると，滞在時間tは，2つの確率変数の商が生成する確率分布d(t)に従う確率変数と見なせる．一般に，確率変数yとxの「商」は，以下のような確率分布に従う（導出は省略．途中，ガンマ関数とか出てきてややこしいが，以下のような形に落ち着くはず¹）．

\[ d(y/x) = \int y f(x) g(y) dy \]

すなわち，f(x)とg(y)が決まれば，yを積分消去することで，滞在時間が従う確率分布をxの関数として定義できることになる．そこまでいけば，滞在時間データから速度が従う確率分布のパラメータを決めることも可能だろう（今回は，最尤法で推定する）．

移動速度の確率密度関数f(x)

　では，移動速度は，どのような確率分布に従うと想定できるだろうか？これについては，イギリスのRowcliffe大先生がパナマの熱帯雨林で詳細に調べた研究から，対数正規分布に従うことが多いということになっている．おそらく，動物は，時々立ち止まるので，裾野長い分布になるということだろう．対数正規分布は，その逆数も対数正規分布に従う（v ~ dlognorm(mealog, sdlog)の場合，1/v ~ dlognorm(-meanlong, sd)になる）という便利な性質があるので，これはうれしい（今回は使わない）．

移動速度の確率密度関数g(x)

　一方，距離の方も（既存の確率分布にはもちろん従わないが）「ある距離になる確率密度」を与える関数を定義することは不可能ではなさそう．最初，これが結構面倒くさそうに感じていたが，超優秀な東大の院生M君と話していたら，割とエレガントな導出方がありそうだったので，やってみることにする．

　正三角形を考える．動物が底辺と左斜辺の任意の点を通過したとする（下図の赤線がその軌跡）．

この時，動物の軌跡bの長さは，正弦定理から以下のようになる．

\[b^2 = a^2 + b^2 - 2 \times a \times c \times cos(pi/3)\]

このbの長さが，どのような分布になるのかをシミュレーションで確かめてみる．辺の長さは，重要ではないので，以下では１辺の長さ1の正三角形の場合を考えてみる．

n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1) # 辺の任意の点
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y) # bの長さ
}

h <- hist(z, freq = F, breaks = 100)
lines(density(z), col = "blue")

　結構面白い形になる．0.9のちょっと前くらいまでは，直線的に増加していき，そのあと，滑らかに減衰していく．

　このようになる理由は，次のように考えるとわかりやすい．長さの決まった棒があり，その棒の両端を正三角形の辺の上に重ねるという遊びを考える．棒の長さが短いと以下の図のように底辺にぎりぎり沿う線から垂線まで動かせる．棒は，左図の黒丸Dから黒丸Bまで動かせると考えるとわかりやすい．棒の置きやすさは，DBの長さに比例するので，棒が長いほど，うまい具合に重ねやすくなる．

　棒がさらに長くなると（三角形の高さ = sqrt(3)/2より長くなると），垂直の棒が置けなくなってしまう（右図の垂線）．この場合，辺ABとの接点である黒丸が動ける範囲は，AからFの範囲になる．

　以上のように考えると，相対的な棒の置きやすさ（=確率密度に比例する値）は，棒の長さをyとすると，

• 0 < b < sqrt(3)/2の範囲では，左図のDBの長さ．すなわち，\(g(y) = 2y/\sqrt{3}\times 定数\)

• sqrt(3)/2 <= b < 1の範囲では，右図のAF

ということになる．

　正弦定理を使えば，AFの長さも容易に求まる．最終的にAFは以下のような形になる（綺麗な式になったらうれしいなと思ったけど，そうでもなかった）．

\[ g(y) = (1 - \sqrt{4y^2 - 3})\times 定数 \]

なお，定数と書いているのは，式で算出されるのは，相対的な生じやすさで，確率密度に比例する値だから．まあ，推定上は無視していい値だが，以下では一応値を出している．

念のためRで確認

　Rで確認する．確率密度にするためには，面積を１にする必要があるので，最初に積分して面積を出して調整する．

# 規格化するための定数の算出（総面積を1にするための準備）
# yが0 ~ sqrt(3)/2の範囲
scale.s <- function(x){
  return(2*sqrt(3)/3*x)
}
integrate(scale.s, 0, sqrt(3)/2) # 念のため確認 sqrt(3)/4に一致

## 0.4330127 with absolute error < 4.8e-15

ss <- sqrt(3)/4

# yがsqrt(3)/2 ~ 1の範囲
scale.l <- function(x){
  return(1 - sqrt(4*x^2 - 3))
}
sl <- integrate(scale.l, sqrt(3)/2, 1)$value # これも頑張れば手計算可能

scale.value <- 1/(ss + sl) # これを掛ければ面積1に


# 距離分布の作成
p.dis.s <- function(x){ # y < sqrt(3)/2
  return(2*sqrt(3)/3*x*scale.value)
}
p.dis.l <- function(x){ # sqrt(3)/2 < y < 1
  #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10)))*scale.value) 
}

# 一応，シミュレーション結果と一致することを確認
n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y)
}

h <- hist(z, freq = F, breaks = 100)
x0 <- seq(0, sqrt(3)/2, by = 0.001)
y0 <- p.dis.s(x0)
lines(x0, y0, col = "red", lwd = 2)

x1 <- seq(sqrt(3)/2, 1, by = 0.001)
y1 <- p.dis.l(x1)
lines(x1, y1, col = "red", lwd = 2)

赤の線が式によって得られたモノ．シミュレーションとぴったりと一致したので，大丈夫そう．

尤度関数の定義

　では，相対的な生じやすさが分かったので，最初の式に戻って，滞在時間の確率密度関数を定義することにする．積分が出てくるが，多分陽には解けないので，Rで近似計算する．関数の中でfor関数を回すのは結構気持ち悪いけど，この辺りは慣れの問題．あと，最尤推定の際に必要なので，対数尤度にマイナスをつけた形で定義している．

# マイナス対数尤度関数の作成
minus.loglf <- function(stay){
  lik <- numeric(0)
  return(function(para){
    for(i in 1:length(stay)){
      int.fs <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.s(d)
      }
      int.fl <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.l(d)
      }
      bon.d <- sqrt(3)/2
      lik[i] <- integrate(int.fs, lower = 0, upper = bon.d)$value + integrate(int.fl, lower = bon.d, upper = 1)$value
    }
    loglik <- -sum(log(lik))
    return(loglik)})
}

最尤推定

　optim関数を使って，最尤推定してみる．

# データ準備
nsim <- 300 # 多めに300個の滞在時間を作成
meanlog <- 1 # これが正解
sdlog <- 0.2 # これが正解
v <- rlnorm(nsim, meanlog = meanlog, sdlog = sdlog) # 速度が対数正規分布に従う

# 移動距離の作成
dist <- rep(0, nsim)
for(i in 1:nsim){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  dist[i] <- sqrt(x^2 + y^2 - x*y) # 距離の乱数
}
hist(dist)

stay <- dist/v # 速度と距離から滞在時間を作成
hist(stay,  breaks = 10)

# 距離分布の作成
p.dis.s <- function(x){ # d < sqrt(3)/2
  return(2*sqrt(3)/3*x)
}
p.dis.l <- function(x){ # sqrt(3)/2 < d < 1
   #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10))))
}
res <- optim(par = c(0.8, 0.15), fn = minus.loglf(stay)) # parは初期値
res

## $par
## [1] 1.032135 0.201001
## 
## $value
## [1] 899.061
## 
## $counts
## function gradient 
##       55       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

hist(v, freq = F, ylim = c(0,5), breaks = 30)
curve(dlnorm(x, meanlog, sdlog), add=TRUE, col = 'blue') # 正解
curve(dlnorm(x, res$par[1], res$par[2]), add = TRUE, col = 'green') # 推定線

　青が真の確率分布，青が推定されたパラメータによるもの．無事に一致した．

　ということで，無事，滞在時間から移動速度のパラメータを決めることが出来た．実際の動物は，直線的に移動するとは限らないので，実用的な価値は分からないが，とりあえずできてうれしい．

---

1. 「はず」というのは，色んなサイト見ながら最後は自分でやったから．割と間違った（少なくとも，推定が上手くいかない）式が書かれていることもあるので要注意．↩︎