「便所の落書き」は楽しかった

　昔、掲示板サイトの2ちゃんねるがよく「便所の落書き」と批判されていて、いまのSNSも同じように言われることがありますが、最近の都会のトイレは総じて綺麗で、落書きというものを見ることがほぼなくなりました。なので、そもそも便所の落書きとはどんなものか、どういう気持ちで書くものなのかを若者に伝える義務があると感じています。

　大学時代、私が主に授業を受けていた建物は1階にトイレがあって、男子トイレはたぶん個室が3つぐらい付いていたと思います。個室の内側の壁にはだいたい何かしら落書きがされているのですが、よく見かけるのは、「○○ちゃんかわいい！」と好きな女の子への愛を叫ぶものや、単なる下ネタや、嫌いな人の携帯番号を書いていたずら電話を誘うものなどでした。電話番号は、大学に限らずトイレの落書ではよく見かけるもので、本当かどうかは分からないけど、クスリの取引やゲイの出会いを匂わせるものもありましたね。

　で、そのトイレの一番左端の個室は、洋式便器に腰を掛けると目の前に見えるあたりに、誰かが「一学2ちゃんねる」と書いた領域がありました。「一学（いちがく）」というのはその建物の名前です。
　「一学2ちゃんねる」は、その個室で用を足す人たちが交換日記のように使っていて、鉛筆やボールペンで

今日もいいウンコが出た
おめでとう、俺もだ
次テストだけど勉強してないわ
がんばれ

というような他愛もない会話が書き込まれていました。自分が何を書いたかは思い出せませんが、だいたい似たようなレベルの、どうでもいいことしか書いていないと思います。

　この「一学2ちゃんねる」には独特な楽しさがありました。読んでいるだけだとあまり何も思わないのですが、自分で書き込むようになると不思議な緊張感を覚えるもので、「次ウンコするときに返信が付いてたらいいな」と少しワクワクするような面もありました。
　匿名のコミュニケーションは「無責任だ」と言われることが多く、それはまぁ確かにそうだと思いますが、お互いの素性が分からないからこそ得られる妙な連帯感というものもあります。少なくとも「一学2ちゃんねる」は、他の個室の壁に比べれば和気あいあいとしていて、荒れることもありませんでした。

　会話の相手が先輩なのか後輩なのかも分からないが、とにかく俺たちはこの同じ便器に座って用を足し、同じ壁を見つめている。そこに「強い絆」というほどのものはないにしても、赤の他人というわけではないし、何より俺たちは互いに話したがっている……。素性を隠しながら同じ目的のために協力する、秘密結社のような感じにも近いかも知れません。

　
　コミュニケーションのチャネルは、太ければ太いほど良いというものでもないのだと思います。家族や同級生のように、頻繁に顔を合わせることで育まれる絆はもちろん大事なのですが、逆に「この人とは、今ここでしか話せないかもしれない」という儚さが、「今ここ」の会話に良い意味での緊張を与える場合もあるんですよね。壁に隔てられて、相手を深く知ることができないからこそ、「話したい」という気持ちが強くなるとも言えます。

　スマホを初めて買った頃、「Qボトル」というアプリを使ったことがあります。ビンに手紙を入れて海に投げると、そのうちどこかの海岸に流れ着いて誰かに読まれるみたいな話がありますが、そのイメージで、誰に届くか分からないメッセージを適当に送信するアプリです。ただしビン投げとは違い、受け取った人がその送り主に返信をすることもできて、少しだけ1対1の会話が成り立つようになっていました。最近は、ランダムに電話をかけて知らない人と話せるアプリなんかもあるようですね。

　そういうものへの需要が「大きい」とは、私は思いません。2ちゃんねる的なものも、たまに使う分にはちょっとワクワクして楽しいというぐらいの話で、毎日入り浸っていたら不愉快な経験のほうが目立つようになってきます。そういう空間に閉じこもってしまうのは、トイレから出られなくなるようなもので、本当に孤立している不幸な人だけでしょう。

　しかし、海に投げた手紙が誰かに届くというエピソードと同じで、「便所の落書き」にも独特なロマンがあるのだということは認めるべきで、それはアカウントによる統合が進んだ今のSNSにはない面白さだったのです。私にとっては、インターネットの「2ちゃんねる」より、本当に便所の壁に書いていた「一学2ちゃんねる」のほうがよりロマンがありましたが、本質は似ていると思います。

　私は高校2年ぐらいのときに「あめぞう掲示板」を使いはじめて、その後、違いがよくわからないまま気づいたら「2ちゃんねる」に移動していました。当時のインターネットには、なんというか、少年の「冒険心」をくすぐるところがありました。知らない人と顔も名前も年齢も分からないまま話し合うっていうのが、ちょっと怪しい、ちょっと危うい、ちょっと悪いことをしているような感覚で、刺激的だったんです。たぶん、その前の時代の「パソコン通信」や、アマチュア無線でたまたまつながった人と雑談をするのも似た感じだったんじゃないでしょうか。

　これはたぶん、インターネットが流行り出した頃、周りよりも早めに使いはじめた人にしか分からない感覚だと思います。使う人が多くなるとトラブルも増えるしコミュニケーションの質も下がってしまって、ドキドキ・ワクワクという感じではなくなったし、そもそもインターネットがあるのが当たり前になると、冒険的ではなくなりますからね。

　私も人に自慢できるほど早かったわけではないし、そもそも「アーリーアダプター」であることが格好いいとも思わなくてむしろ古風な生活スタイルを守っている人に憧れますが、たとえば大学に入った時点で下宿先に常時接続の回線があったのは私ぐらいで、インターネットといえば大学のメディアセンターで使っている人がほとんどだったので（携帯でも接続はできましたが当時は見れるサイトが限られてました）、比較的早いほうだったとは思います。

　で、わずか数年差ぐらいの話ではあるのですが、アーリーアダプター層だけが集まっていた頃の掲示板やチャットルームは、「匿名で怪しいけどなんか凄そうな奴」がいっぱいいて面白かったんです。たとえば、吉野家コピペの元ネタを作ったのは新爆さんという人ですが、吉野家コピペが有名になるよりも前（だったはず）に、お笑い好きが集まる個人サイトの掲示板で新爆さんの書き込みを見て、腹がちぎれるほど笑ったのを覚えています。ちなみに、吉野家コピペが出回るようになってからも15年ぐらいの間、私は元ネタの作者が新爆さんだとは知らなくて、あのコピペとは無関係に新爆さんというキャラクターを記憶してました。

　SNSが登場して以降のインターネットは、「（凄くて）怪しくない人」と「（怪しくて）凄くない人」ばかりになって、私自身はもちろん後者ですが、とにかく当時のような不思議な味わいが得られる機会は減りました。それは残念な面もあるのですが、避けられない道であったとも思うし、それこそ数年限りの「儚い」体験だったからこそいい思い出になっているのかも知れません。

（参考：このブログのおすすめ記事一覧はコチラ）

天体望遠鏡のトレンドが90年代とはだいぶ変わっていた

　最小二乗法の歴史を調べていたら（先日のエントリ）、かつて統計学を発展させたのは「天体観測」と「測量」だったのだなぁということを改めて実感したのですが、そういえば私は中高生だった1990年代に天体写真を撮っていました。で、「そういえばいまの望遠鏡のトレンドってどうなってるんだろう？」と思ってメーカーのサイトを見てみたら、当時とはかなり考え方が変わっているなと思いました。一言でいえば、「屈折式の地位がめちゃめちゃ上がってる」ということです。

　まず前提として90年代当時の「天体望遠鏡選び」の基本を整理しておくと、

• 屈折式：対物レンズと接眼レンズから構成されるもので、理科の授業とかで触れるのも、漫画とかで描かれるのもこのタイプが多いが、大口径化することが難しい（レンズは鏡に比べて高価）ので、写真撮影ではあまり使っている人がいなかった。
• 反射式（ニュートン式）：対物レンズのかわりに主鏡で光を集めて鏡筒の横から覗くやつ。大口径のものでも比較的安価なので、太陽観測などを除けばどのジャンルでもこれが最もポピュラー。密閉されてないので、主鏡にホコリがついたり結露したりするのを気にする必要があり、光軸調整も必要なので、扱いはけっこう面倒。
• シュミットカセグレン式：屈折と反射のいいとこ取りで、鏡で光を集めるので屈折式よりは大口径化しやすく、鏡筒が密閉されているので扱いがちょっと楽。ただし反射式より高価なので、金持ちの（大人の）贅沢品というイメージなのと、焦点距離が長めでちょっと暗い。
• ドブソニアン式：大口径化に振り切ったやつで、雑誌などで見かけることはあるし面白いが、マニアックすぎてふつうは選択肢にならない。

という感じでした。中高生だった頃（高校時代は競馬のほうが頑張っていたのでほぼ中学時代）の知識なので不十分とは思いますが、天文ガイドのような雑誌は隅々まで読んでいたのと、専門書もけっこう買い揃えていたので（光学などの細かい理論は理解できませんでしたが）、大きく外れてはないと思います。

　この当時の一番のポイントは、とにかく撮影時に長時間露光をする必要があるので、光量をかせぐために「大口径が正義」だったことです。
　当時は、デジタルの冷却CCDカメラというものが少し出始めてはいましたが、フィルム撮影が主流でした。また、パソコンでの画像処理（当時は「電子暗室」とか呼ばれてましたｗ）環境も持っていない人が多かったと思います。色調の調整をするなら、自宅に暗室をつくって現像とプリントで頑張るという感じでした。
　フィルムは、現代のデジカメのセンサーに比べると感度が弱いので、暗い天体を撮る場合は長時間の露光が必須でした。私は月面のクレーターの拡大写真と、木星・土星ばかり撮っていたのですが（スペック的に土星はかなり厳しかったですが）、月のような明るい天体でも強く拡大すると暗くなるので、1/2秒から2秒ぐらいは開けてたと思います。惑星になると、うろ覚えですが、10秒ぐらい露光することもあったような。

　経験がない人にはなかなかイメージできないと思いますが、長時間露光するときは、カメラ本体のシャッターを使わず、文字通り「手動」で露出時間を調整します。どういうことかというと、まず黒く塗った板を手で持って、それを望遠鏡の先っぽにかざして光が入らないように塞ぎます（鏡筒に接触はしないように）。次に、接眼レンズにマウントしたカメラのシャッターを解放します。この時、望遠鏡が微細に振動するので、それが収まるまでしばらく待ちます。そして、もういいかなと思った頃に、かざしている黒い板をどけて露光を開始し、1/2秒とか1秒とか数えて塞ぎます。その後に、カメラ自体のシャッターも閉じる。

　光量確保のために長時間の解放が必要で、何百倍という拡大率で撮影していると少しの振動でもボヤけてしまうので、こういう手順で撮るのが一般的だったと思います。で、フィルム撮影だとその場で写真の確認はできませんから、何十回も撮っておいて、後日現像に出して運よく良いものが撮れているか確認するという流れです。露出時間もいろいろなパターンを試して、何枚目が何秒だったかというのをノートに記録しておき、その後の参考にします。
　ちなみに最近は、惑星や月の撮影で長時間露光はしないみたいですね。デジタルで連続撮影してきれいなものを重ね合わせるという手法が強いようです。

　『天文ガイド』の18歳未満の部に3度掲載してもらいましたが、そのうちの一枚がこれです。少しぼやけていて褒められるようなクオリティではないですが、11.4cmという小口径の望遠鏡で中学生が撮ったものなので、大目に見てください（笑）

　集光力は対物レンズや主鏡の口径の二乗に比例しますので、上述のような環境だと、とにかく「大口径が正義」ということになります。光量があればあるほど、拡大率を引き上げることもできますし、露光時間が短時間で済むので振動（や大気のゆらぎ）の影響を受けにくくなって、綺麗な写真が撮れます。惑星とかを狙うのであれば、だいたい20cmクラスぐらいの望遠鏡を持っていないと雑誌の投稿欄で掲載を争えるような写真は撮れない感じでしたが、私とおなじ機種で、ガイドスコープによる補正などを駆使して木星か何かのきれいな写真を撮っている人もいたので、気合で乗り切れる面もあるのだなと感心したりしていました。

　まぁいずれにしても、大口径原理主義になりがちなので、口径あたりのコストが安いニュートン反射式が圧倒的に主流でした。色収差がないというメリットもあるらしいですが、そんなことよりまずはとにかくコスパでしょう（笑）。屈折式は眼視がしやすいので、写真以外の用途ではよく使われていたと思います。眼視のほうが明るく見えやすいんで。

　今でも、コスパの良さから、大口径の反射式を買ってより本格的な（というか暗くて撮りにくい）天体を……というムーブは当然あると思いますが、デジカメのセンサーがかなり優秀で、かつてのフィルムより感度がいいようです。そうなると、昔ほど大口径にこだわらなくてもよく、扱いが楽な屈折式が写真撮影用途で現実的な選択肢になってくるんでしょうね。
　もともと屈折のほうが、筒内気流や反射時の散乱がなく、コマ収差と言われる滲みもないので、原理的に反射式よりシャープに撮れるという強みがあって、それが活かせるようになったのだと思います。しかも、アポクロマートという技術で色収差の問題がほぼなくなったみたいですし。
　星雲などで長時間露光するにしても、ガイドスコープを使ったオートガイドの精度が上がっているようで、そのことも「小口径でも撮れる」を後押ししてるのかも知れません。色収差の補正技術が向上した結果、短焦点化が可能になり、明るく撮れる屈折式が増えたってのもありそうです。

　いずれにしても、私は天体写真撮影用の望遠鏡としては、「ニュートン反射一強」という認識だったので、屈折式の地位が上がっていることを知って驚きました。市場シェアとしても、屈折式のほうが売れているそうです（参考リンク）。昔は、月でも惑星でも星雲星団でもニュートン反射がバンバン使われてましたが、今は月・惑星では屈折式かシュミットカセグレン式が使われていて、反射式はどちらかというと星雲星団（いわゆる「深宇宙」）撮影用という雰囲気にシフトしてきているようです。時代は変わるものですね〜。

（参考：このブログのおすすめ記事一覧はコチラ）

統計的因果推論の手法を選ぶためのチートシート

　統計的因果推論について、学生とかが「手元にどんなデータがあるか」に応じて使える手法を大雑把に判断するためのチートシートのようなものを作っているのですが、正直私自身もあまり理解してなくて、正確性を気にしていたらキリがないと思い、いったん現行版を掲載しておきます。
　間違いに気づいたら随時修正します。無理やり一枚にまとめようとした感じなので、分類の切り口を含めて、もっといい整理法がありましたら教えてください。

　まず、一番大雑把な表がこちらです。【オリジナル画質で見る場合はコチラ】

　情報をいろいろ補ったものがこちらです。【オリジナル画質で見る場合はコチラ】

　「こういうデータにはこの手法しか使えない」という意味ではないです。工夫すればいろいろなタイプのデータに適用可能な手法も多いので、あくまで「こういうデータを持っているときに、とりあえず検討してみるのはこの手法かな」と見当をつけるための早見表です。

　上の段と下の段の分け方は、理論的にはあまりスッキリしてないとは思います。2SLSや固定効果モデルを、Rubin的な介入効果の分析にも使えるわけなので。
　何というか、ぶっちゃけて言うと私が属している分野（都市計画周辺の工学）の空気感としては、

• 上の段は「最近耳にするようになった統計的因果推論の話題に、意識的についていこうとしないと、よく分からないままになる手法」
• 下の段は「伝統的な統計分析の教育を受けてふつうに生きていれば、どこかで触れることがありそうな手法」

という感じになりまして（笑）、理論的な正確さや明快さより、そういう日常的感覚みたいなもので分けておきたかったというのが大きいです。あまり統計分析の手法について突き詰めて考える人は多くない分野なので、私を含めて認識は雑だと思います。
　ちなみに私自身、実際に回したことがあるのはSCM・DID・PSM・2SLS・GMMぐらいで、論文で使ったことあるのはたぶんSCM・GMMだけだなので、とても正確に理解してるとは言えません。自分が今後勉強するための、頭の整理という感じでもあります。

　画像だけというのもイマイチなので、以下にテキストでも書いておきます。

Rubin型

　処置＝介入の概念が明確で、「反実仮想との差」を評価することを強く意識した手法です。最近流行りの「統計的因果推論」というと、これらをイメージすることが多いと思います。

処置群と対照群	時系列データ	割当ルール	反実仮想の構成方法	手法
別個体	あり	―	並行トレンドを仮定し対照群との差を評価	差分の差分法（DID）
別個体	あり	―	対照群の加重平均や共変量との関係から合成反実仮想を構成	合成コントロール法（SCM）
別個体	なし	明示的ルールなし	1次元化した指標で個体同士をペア比較	傾向スコアマッチング（PSM）
別個体	なし	明示的ルールなし	交絡因子ごとに層別し比較後に統合	層別解析
別個体	なし	量的変数による閾値ルール	閾値近傍の局所比較	回帰不連続デザイン（RDD）
同一個体	あり（前後のみ）	―	切片や傾きの変化を評価	中断（分割）時系列解析（ITSA）
同一個体	あり（前後のみ）	―	ベイズ構造時系列で柔軟に反実仮想を構成	Causal Impact

構造推定型

　伝統的な経済モデルとかで使われているもので、理論モデルを前提として、内生性を統制して構造パラメータの推定を行うものです。反実仮想を明示的には構成せず、介入の概念が明確でないケースも扱いますが、反実仮想との差を評価していると解釈することができるケースもあると思います。

データ構造	ラグ項の利用	推定アプローチ	手法
クロスセクションデータ（複数個体・1時点）	なし	外生的な操作変数を用いる	2段階最小二乗法（2SLS）
1個体の時系列データ	あり	ラグ項を操作変数として利用	ラグ項を用いた2段階最小二乗法
パネルデータ（複数個体の時系列）	あり	ラグの外生性を仮定しGMM推定	動学的パネルデータ分析（GMM）
パネルデータ（複数個体の時系列）	なし	個体固定効果により時間不変の未観測要因を除去	固定効果モデル

各手法のイメージ

表に載せている各手法のイメージに関するメモも、テキストで書いておきます。
代表的な文献の引用は、AIにチェックさせたら「それはワーキングペーパー版であって、査読論文は年号がXXXXになる」という指摘が3つぐらいありましたが、面倒なのでそのままにしています。

手法	特徴・本質的ポイント	代表例
差分の差分法（DID）	トレンドの平行性をチェックする手続きが本質的に重要（介入前に複数時点がないとチェックできない）	ニュージャージー州の最低賃金引き上げはファストフード店の雇用に影響を与えたか（Card & Krueger, 1993）
合成コントロール法（SCM）	万能ではないが反実仮想の構成の手法としては最もシステマティック	テロの発生はバスク地方の1人あたりGDPに影響を与えたか（Abadie & Gardeazabal, 2003） 阪神淡路大震災は兵庫県・神戸市の1人あたりGRP等に影響を与えたか（duPont & Noy, 2015）
傾向スコアマッチング（PSM）	比較対象の適切なマッチングによって交絡を統制する（確率的なモデリングを行う）	失業者、受刑者、若者などへの就業支援のための試用プログラムの効果の再分析（Dehejia & Wahba, 1999）
層別解析	比較対象の適切なマッチングによって交絡を統制する（記述統計的なグループ分けを行う）	喫煙が肺がんの原因であることを主張（Cornfield et al., 1959）
回帰不連続デザイン（RDD）	閾値近傍の局所ランダム化の仮定により、PSMや層別解析より強力に因果を識別	アメリカの高校卒業資格試験取得（試験の点数が強制変数）が市民的行動（投票参加や新聞購読）に与える影響（Dee, 2004）
中断（分割）時系列解析（ITSA）	伝統的な回帰モデルを用いて、トレンドや季節性を考慮し、介入前のデータの外挿を反実仮想とする	シチリア島での禁煙政策が急性冠動脈の発生に与えた効果（Bernal et al., 2017）
Causal Impact	ベイズ構造時系列モデルを用いることで、トレンド、季節性、共変量との関係から反実仮想を柔軟かつ不確実性を含めて構成	オンライン広告がサイトのアクセス数に与える影響（Brodersen et al., 2015）
2段階最小二乗法（2SLS）	外生性を主張できるよう適切な操作変数を探す（操作変数の妥当性は理論依存）	教育年数が収入に与える影響を推定するために、州ごとの義務教育年数と個人の生まれ月を操作変数とした（Angrist & Krueger, 1991）
ラグ項を用いた2段階最小二乗法	操作変数にラグ項を用いることも可能だがマイナー（ラグの活用はむしろ動学パネルの手法というイメージ）	一応やっている研究者もいる程度
動学的パネルデータ分析（GMM）	ラグの外生性（強い仮定である点に注意）を利用し操作変数としてGMM推定量を求める	企業の賃金水準、生産高、資本金が雇用数に与える影響（Arellano & Bond, 1990）
固定効果モデル	時間不変・個体不変の未観測交絡因子は考慮できるが同時性は統制できない	社会的移動性が何に影響されるかを分析（アメリカの741地域×17年のパネルデータ）（Chetty et al., 2014）

注記

　表の下に言い訳のようにダラダラ書いてるメモもコピペしておきます。

• この表は、どんな性質のデータを持っているかを前提として、左から右にみていく想定で作ったもの。たとえばCausal Impactは、対照系列を用いて識別力を高めることもできるので、「同一個体の処置前後のデータしかない場合」以外にも使える。
• 一般的に統計的因果推論と呼ばれるのはRubin型のものだが、ここでは未観測の交絡因子（原因と結果の双方に影響する因子）や同時性等からくる内生性を統制する戦略、つまり原因変数の外生性を確保するための明示的な戦略を持った方法で構造パラメータを推定する手続きも因果推論とみなした（行動モデルや経済モデルなどで、介入の概念が明確でないケースでも因果の議論が可能という利点がある）。
• 単にOLS回帰（に基づくSEM、ケインズ型マクロモデル、一般均衡モデル等）で効果を求めるだけだと、観測可能な共変量についてしか統制できないので、理論依存が強く、因果関係を主張できる度合いは、識別のための仮定を吟味するこれらの手法より劣る。（計算の手続き自体はDIDもOLSの一種ではある。問題は、因果の識別のための条件を意識しているかどうか。）
• Rubin型のうち、PSMと層別解析、ITSAとCausal Impactは、観測された交絡しか統制できない（Causal Impactは共変量の工夫で一部統制可能）が、反実仮想を構成する明示的な枠組みを持ち、その妥当性を確認する手続きがあることが、単なるOLS回帰にない利点だと言える。
• 因果の識別（モデルの識別とは別概念なので注意）能力は結局程度の問題で、どれも完璧ではなく、一概にどの手法が優れているとは言えない。
• 構造推定型は、反実仮想との対照を明示的には行わず、「理論を信じてパラメータを推定」という感じだが、2SLSで適切な操作変数が存在するケースでは、強力に因果を識別できるといえる。操作変数法は、潜在アウトカムとの比較として定式化すれば、介入効果（局所平均処置効果）を見ていることになる。
• 構造推定と言っても、ケインズ型の同時方程式モデルでは内生性の統制が中心課題だが、CGEやDSGEにおいては原因の外生性は理論的に仮定する場合が多く、内生性の統制はあまり問題にしないらしい。
• 2SLSやGMMは操作変数（GMMのラグ操作変数はかなり強い仮定だが）を使って原因変数の外生性を確保するのに対し、固定効果モデルにはそれがないのだが、時間不変・個体不変の未観測の交絡因子を統制するのも立派な因果識別戦略だと思う。

（参考：このブログのおすすめ記事一覧はコチラ）

なぜ誤差を「二乗」するのか？

　何年か前に、統計学の勉強会をしている学生の会話を聞いていたら、パラメータの推定に「二乗誤差」を用いる理由を「誤差の符号を正にするため」というふうに先輩が後輩に説明していました。理由を考えようとする姿勢は素晴らしいと思いつつ、「そう単純でもないんだよな」ということで、分散の加法性がどうたらこうたらとコメントした気がしますが詳細は忘れました。

　観測誤差は、真の値からプラスの方向にもマイナスの方向にも振れることが多く、仮に正規分布のように対称な分布になっているとすると、単純な足し合わせではゼロに近づいてしまって、「誤差の大きさ」（観測の精度）を議論できなくなります。「じゃあ、絶対値を取ればいいじゃないか」という発想はあり得て、実際にロバスト回帰という手法では絶対偏差（誤差の絶対値の総和）が使われるのですが*1、一般的には「二乗誤差」で評価することがほとんどです。

　では、それは何故なのかという話になると、Wikipediaにも色々なことが書かれてますが、2つの方向での説明があり得ます。1つは、最小二乗法が使われるようになった歴史的経緯をたどり、それを考案した人たちにとってどういう理由があったのかを考えるということ。もう1つは、現代の統計学に接続する意味で、理論的な「扱いやすさ」を確認するということです。これらは、「使われ出した理由」と「使われ続けている理由」に対応するとも言えるかも知れません。

　理論的な扱いやすさのほうを先に挙げておくと、二乗誤差というのはL2ノルム、つまり内積を定義できる空間を考えることにあたり、いろいろと便利な性質を持ちます。
　たとえば、内積が定義できる空間を考えればデータを直交分解できるってのがあります。直交分解というのは、あるベクトルを、「ある空間への射影」と「それに直交する成分」に分離するということですが、これは「回帰直線と残差への分解」もそうですし、分散分析でデータの分散（情報量）を「群内平方和」「群間平方和」「残差平方和」にキッチリ分けることができるのも、内積が定義できる空間を考えているからこそです。主成分分析で各主成分を直交するように求めていくのも、同じ手続きだと言えますね。
　もう一つの理論的な理由は、二乗誤差の最小化が、正規分布の最尤推定に一致するということが挙げられます。これはガウスが1809年に主張したことですが（私はちゃんとは読んでないですが）、最小二乗法の登場よりは後の話なので、後づけの正当化の一種だと言えます。
　あと、絶対値だと原点で角になるので微分できませんが、二乗誤差は微分できるというのも大きいと思います。二乗誤差を偏微分して正規方程式の連立を解いて最適化という流れは、後述のルジャンドルも使っていますね。

　では歴史的にはどうだったのかというと、まずニュートンは1700年前後の時期に、いろいろな天体観測を行う中で、同じ量を何度も計測した場合は、「観測値の平均」を取って代表値としていたそうです。計測機器の性質から、観測精度に差があったので、その差を重みとして考慮した重み付き平均も使っていたらしい。平均値というのは、じつは結果的に「二乗誤差を最小にする点」になりますので、ニュートン自身が「二乗誤差を最小にしよう」というような発想を持っていなかったとしても、結果的に最小二乗法が実践されていたことになります。
　ちなみにWikipediaの記事では、ニュートンが「正規方程式」を初めて書いたと書かれていますが、arxivの論文をたどると、ニュートンはOLSの2本の正規方程式（回帰直線がデータの重心を通り残差の和がゼロになるという条件と、説明変数と誤差の相関がゼロになるという条件）のうち、1本目には到達していたという話です。これがどれだけ「惜しかった」のかは、私にはにわかに判断が付きませんが、2本目の方程式のほうが本質的でかつ到達するのが難しい気はします。

　その後の時代になりますが、ラプラスは1780年代から1790年代にかけて、絶対偏差を最小化する方法を使っていたそうです。現代でいうロバスト回帰の手法ですね。ラプラスは、誤差の確率密度関数（ラプラス分布という指数分布の一種を考えていたらしい）を使って正当化しようとしたらしいですが、理論的にうまくはいかなかったようです。
　この時代に重要だったのは、絶対偏差と二乗誤差のどちらが優れているかというよりも、「たくさんの観測値から代表値を求めると誤差を打ち消すことができる」という発想が定着したことですね。1つの観測の精度を高めることも大事だけど、誤差が対称に分布しているのだとすれば、たくさん測ることで誤差を大幅に減らすことができる。私も、大学で測量の野外実習を教えていたことがあるのですが、精度の低い古風な手法で測っても、平均を取れば思った以上にキッチリ特定の値に近づいていくので感動しました。

　ルジャンドルというのは、1805年に最小二乗法をはじめて理論化したとされる人で、Stiglerという人の論文に詳しく書いてあるように、ガウスとの間で「どっちが先に考案したか」で争いがあったそうです。一応、先に公表したのは間違いなくルジャンドルということになるらしいのですが、ルジャンドルがなぜ最小二乗法を推奨したかというと、二乗誤差を最小化する問題を偏微分すると線形の正規方程式が得られて、この解を求めるというアプローチであれば、未知のパラメータが複数あっても最適化ができるからです。要するに、歴史的には、「パラメータが複数あるケースでも解析的に解ける」という実用上の理由が大きかったのだということです。

　ガウスは自分では「1795年から最小二乗法を使っていた」と主張したらしいのですが、公に発表したのは間違いなくルジャンドルということになります。ただ、ガウスが1809年に自身の手法を公表した際には、正規分布との関係性が議論されており、理論的には間違いなくガウスのほうが洗練されていたようです。そして、上述したように最尤法との一致やら直交分解やらの便利な性質から、現代に至るまで「最小二乗法」が使われているというわけですね👍️
　ガウスの論考は、もとはラテン語らしいですがここに英訳があった（長いけど該当箇所は後ろのほうの172節〜179節）ので、理解できる気はしないけどあとで眺めておこうと思います。178節で正規分布が導出されていますね（で、が観測誤差です）。

（参考：このブログのおすすめ記事一覧はコチラ）

*1:二乗誤差よりも外れ値の影響を受けにくい。これは原理的には、平均値というものが二乗誤差を最小にする代表点、中央値というものが絶対偏差を最小にする代表点であり、中央値のほうが外れ値の影響を受けにくいことに対応している。

95%信頼区間は何が「95%」なのか？

　10年以上前に、「信頼区間の意味と、Rのpredict()関数の使い方の注意点」というエントリを書いてましたが、余計な記述が多かったのであらためて論点を整理しておきます。

　統計学の教科書をまじめに読んでいる人にとっては今さらな話なのですが、「95%信頼区間」というのは、「この範囲に真の値（母数）がある確率が95%」ということを表しているのではありません。ところが、学生のような初学者だけでなく、理系の大学教員でも「この範囲に母数が含まれる確率が95%」というふうに間違った説明をしてしまっていることがあるので、意外と厄介な概念です。

　たとえば単回帰分析を行って、約2.5という回帰係数が得られ、その95%信頼区間が約1.8〜3.2だったとします。

　ここで、「真のパラメータが1.8〜3.2である確率が95%」と言ってしまうと間違いだということになります。頻度主義の統計学は、「真の値が所与の場合に、データがどう変動しそうか」を考えるものなので、「真の値がこの範囲で動く」という考え方は基本的には持ち出しません。ベイズ統計学の場合は、「データが所与の場合に、真の値はどのあたりにあると考えるのが確からしいか」と考えるので、真の値が動きますが、ベイズ統計学の立場を取っているわけではない文脈でというか、おそらくベイズ統計学を意識的に勉強したわけでもないのに、そういうイメージで語ってしまっている人がたまによく居ます。

　頻度主義における「信頼区間」を正しくイメージするなら、

• 真のパラメータ（母数）を設定してやると、それに従っサンプルデータを生成する、乱数発生器を考える。
• サンプルを生成するたびに、一定の手続きに従って、そのサンプル（データセット）から母数を推定し、信頼区間も計算する。
• この作業を何回も繰り返すと、信頼区間そのものは、毎回異なるものが得られる。「サンプリングと信頼区間計算」を100回繰り返すと、100個の信頼区間が得られる。
• この100個の信頼区間のうち、95個ぐらいは、最初に設定した真のパラメータを含んでいる。

　ということになります。一番のポイントは、手元に得られたデータから計算された1回1回の信頼区間の値そのものには、あまり意味はないということです。「今回得られた区間に95%収まるように何かが動く」と考えるから、間違えるわけです。そうではなく、サンプリングを繰り返すたびに「信頼区間自体が動く」んですよね。

　別の言い方をすると、「95%」というのはなんとなく強そうな数値ですが、「今回得られた信頼区間」が95%の強さを持っているわけではないんです。そうではなく、「信頼区間を計算する手続き」（ルール）が、何度も実行されたときに、95%の確率で真の値を捕捉することができるという強さを持っているということです。95%というのは、いわば「信頼区間算出ルールの勝率」です。この、「算出ルールの勝率」という概念を持っておくと、うまく頭が整理されるのではないかと思います。

　まぁ、「真のパラメータが1.8〜3.2である確率」を意味するベイズ信用区間も、標準的な設定（共役事前分布を使うとか）で計算すると結果的には一致するし、そっちのほうが直観的に理解しやすいのですが、頻度主義の統計手法を使っている限りは、明確に異なる概念であるということは理解しておく必要があります。

　じつは上述の信頼区間は、が真のパラメータだとして、という帰無仮説からという帰無仮説までは、今回のデータによって棄却されないという意味にもなります。また、真のパラメータが1.8や3.2のとき、そのモデルからサンプリングを繰り返してパラメータを推定する操作を繰り返すとした場合に、真のパラメータからサンプリングされたデータから計算されるパラメータ推定値のばらつきの両側95%の範囲内に、今回の手元のデータから推定されたパラメータが入っているということでもあります（通常の線形回帰の場合）。
　以前のエントリではその説明を先にしていたのですが、この言い方だと、真の値が1.8から3.2までの間で動くかのような印象を与えるので、やめといたほうがいいですね。この説明においても、あくまで「真のパラメータが1.8に固定されている場合」から「真のパラメータが3.2に固定されている場合」までを考えているということであって、真のパラメータを変動するものとして考えているわけではないのですが、ややこしいです（笑）

　上の図を描いたコードは下記のとおり。回帰係数と信頼区間が特定の値になるように無理やり選んでいます。

library(ggplot2)

set.seed(1)

# 設定
target_beta <- 2.5
target_lwr  <- 1.8
target_upr  <- 3.2
target_half <- (target_upr - target_lwr) / 2  # 0.7

tol_beta <- 0.002   # 傾きの許容誤差
tol_half <- 0.002   # 半幅の許容誤差
max_tries <- 20000  # 試行回数

# x と ノイズ規模
n <- 50
x <- seq(1, 25, length.out = n)

beta0 <- 0
beta1_true <- 2.5
half_width <- target_half  # 0.7

Sxx <- sum((x - mean(x))^2)
tcrit <- qt(0.975, df = n - 2)

# sigma を逆算
sigma <- half_width * sqrt(Sxx) / tcrit


# 探索ループ（中心と半幅で評価）
best <- list(score = Inf)

for (i in 1:max_tries) {
  y <- beta0 + beta1_true * x + rnorm(n, 0, sigma)
  df <- data.frame(x = x, y = y)
  
  fit <- lm(y ~ x, data = df)
  beta_hat <- unname(coef(fit)["x"])
  
  ci <- confint(fit)["x", ]
  lwr <- unname(ci[1])
  upr <- unname(ci[2])
  half <- (upr - lwr) / 2
  
  # 距離（中心と半幅だけで評価）
  score <- abs(beta_hat - target_beta) + abs(half - target_half)
  
  # ベスト更新
  if (score < best$score) {
    best <- list(
      score = score,
      df = df,
      fit = fit,
      beta_hat = beta_hat,
      lwr = lwr,
      upr = upr,
      half = half,
      iter = i
    )
  }
  
  # 十分近ければ終了
  if (abs(beta_hat - target_beta) <= tol_beta &&
      abs(half - target_half) <= tol_half) {
    message("Found at iter = ", i)
    break
  }
}

# 結果表示
best$iter
best$beta_hat
c(best$lwr, best$upr)
best$half


# 予測値＋信頼区間（塗りつぶし）
newx <- data.frame(x = seq(min(best$df$x), max(best$df$x), length.out = 200))
pred <- predict(best$fit, newdata = newx, interval = "confidence", level = 0.95)
pred_df <- cbind(newx, as.data.frame(pred))  # fit/lwr/upr


# プロット
ggplot() +
  geom_point(data = best$df, aes(x = x, y = y)) +
  geom_ribbon(
    data = pred_df,
    aes(x = x, ymin = lwr, ymax = upr),
    alpha = 0.25,
    inherit.aes = FALSE
  ) +
  geom_line(
    data = pred_df,
    aes(x = x, y = fit),
    linewidth = 1,
    inherit.aes = FALSE
  ) +
  labs(
    title = paste0(
      "Slope = ", round(best$beta_hat, 4),
      " / 95% CI = [", round(best$lwr, 4), ", ", round(best$upr, 4), "]"
    ),
    x = "x",
    y = "y"
  )

（参考：このブログのおすすめ記事一覧はコチラ）

「検定を繰り返すと第一種の過誤が増える」は間違い

　やや煽り気味のタイトルを付けてみましたが、そんな大層なことをいいたいわけではありません。
　以前、検定の繰り返しと多重比較についてというエントリを書いていたですが、こんなにくどくど書かなくてもいいよなと思ったので、改めて話を簡潔に整理しておきます。

　私が言いたいのは、「繰り返してはダメ」という言い方をしてしまうと、本質が見えなくなって何が問題なのかよく分からなくなるということです。たとえばこの記事に、「この検定を1回だけ行えば、通常は表が多く出たか少なく出たかを正しく判断できます。しかし、同じ検定を何度も繰り返すと、たまたま表が多く出たケースが出てくる可能性が高くなります。」と書かれているのですが、こういう、「繰り返すから悪い」式の説明はよくないと思います。

　心理学をやる人は必ず、頻度主義の統計学で「検定の繰り返し」や「多重比較」について習っていると思うのですが、いちばん分かりやすいのは分散分析の下位検定として、テューキーの多重比較検定やボンフェローニ補正をさせられるケースですね。
　たとえば、小学校にA・B・Cという3つのクラスがあって、それぞれのクラスのテストの成績が有意に異なっているかを検定するとします。そこで1要因3水準の分散分析によって「有意差あり」という結論が出ると、当然、A⇔B、B⇔C、A⇔Cのどこに差があるのかを知りたくなります。そのときに、A⇔B、B⇔C、A⇔Cのそれぞれについて個別にt検定をやると、これは検定の繰り返しにあたり、第1種の過誤（本当は受け容れるべき帰無仮説を棄却してしまう、つまり本当は有意ではないのに有意だと結論付けてしまう）をおかす危険性が増えてしまうので、やめなさいというわけです。しかしこの説明は、誤解を生みやすいです。

　まず、分散分析の下位検定の話の前に、もっと単純な検定の繰り返しの話を整理しておきましょう。たとえば、先ほどの記事で論じられている、コイントスの実験をイメージすると分かりやすいと思います。
　ある実験を行って有意水準5%で「統計的に有意です」と結論づけるのは、「本当は差がないとすれば20回に1回ぐらいの確率でしか見られないほど極端な差が観察されています」と判断することに等しいわけですが、これは逆に言えば、今回観察されたのは「本当は差がなくても20回に1回ぐらいは観察され得る程度の差である」ということです。つまり、実験を20回やれば、本当は差がなくても「有意差あり」と主張できてしまい、それはよくないということです。

　しかし、単に「検定を繰り返すのが悪い」と言われてしまうと、「じゃぁ、異なるテーマでいろんな実験をやって論文を何本も書くと、第1種の過誤をおかす確率が上がってしまうから悪いの？」というような無意味な議論が生じ得ます。もちろん、そんなことはありません。

　上の例の場合、本質的には何が悪いのかというと、有意差が得られなかった19回分のデータを無視して、有意差が得られた1回分のデータだけをみて結論付けようとしているのが悪いんです。「19回は有意差がなかったんですけど1回は有意差がありました」と報告しているなら、べつに問題はないと思います。
　普通は、20回実験を繰り返してデータが得られているのであれば、その20回分のデータを統合した分析をしなければならないということになります。統合にひと工夫必要な場合、「メタ分析」と呼ばれる一連のノウハウを勉強する必要があり、ボンフェローニ補正もそのノウハウの一種だと言えます。
　さらによく考えると、20回分のデータセットからランダムに1セットを選んでそれを分析するなら公平だと言えるのですが（そんな変なことは普通しませんが、別々の研究者がお互いを知らずに同じ実験をしている場合はこれに相当するでしょう）、有意差が得られたやつだけを選んで報告するとなると、自分の研究意図にとって都合のいいサンプルを恣意的に選んでいることになるので、そこに大きな問題があるわけです。同じ原理で、有意差が得られた論文だけ査読を通るようなシステムだと、全体として第一種の過誤が増えてしまうという「出版バイアス」の問題が生じるわけですね。

　分散分析の下位検定としての、多重比較の話に戻りましょう。たとえば、当初はクラスAの成績とクラスBの成績だけが得られていて、両者のあいだでt検定を行って「5%水準で有意差あり」と結論付けたあとに、クラスCでも同じテストが行われてA⇔CやB⇔Cの比較が可能になったとすると、当初のA⇔Bのt検定における「有意差あり」という結論は、間違っていたことになるのでしょうか？

　じつはそんなことはなくて、研究上の関心があくまで「A=B」という帰無仮説を棄却することにあるなら、依然としてA⇔Bのt検定をやっておけば問題はないんです。問題が生じるのは、ABCをあわせた分散分析に接続しようとする場合です。
　ABCをあわせた分散分析というのは、帰無仮説が「A=B=C」になっています。で、「A=B=C」という帰無仮説を棄却する目的でA⇔Bの検定を使うなら、ふつうのt検定ではいけません。というのも、A⇔B、B⇔C、A⇔Cのうち「少なくとも1つの組み合わせ」に有意差があれば「A=B=C」という仮説は棄却できてしまうので、棄却のハードルが下がってしまうことになるからです。

　つまり、問題は「検定を繰り返すこと」にあるわけではないし、「部分的な比較を行うこと」にあるわけでもない。検定対象となる帰無仮説が「A=B=C」であるときに、「A=B」という別の帰無仮説の検定結果をそのまま使ってしまうと、整合性がなくなるよというだけの話です。

（参考：このブログのおすすめ記事一覧はコチラ）

不偏分散の分母がn-1であることの直観的な理解

　同じタイトルで数ヶ月前にエントリを書いていたのですが、よくよく考えると全然意味のある説明になってなかったので、まとめ直します。
　不偏分散が

であることの、証明というか導出は、統計学の教科書を見れば書いてあるのですが、「なんで1引くの？」というのを直観的に理解したいと思ったことがある人は多いと思います。「データのばらつきは、独立したランダム（確率的）な変動が生じる情報空間の次元（広さの一種）を意味する自由度に応じて評価すべきもので、パラメータを1個（ここでは母平均）推定するとデータの変動に決定論的な仮定（拘束）が持ち込まれて自由度が1下がるから」というような説明で、直観的に腑に落ちるという人はいいのですが、そういう人はたぶん少数派でしょう。
　それで整理の仕方を考えていて、ひとまず以下のようにまとめましたが、これで腑に落ちる人もあまりいないかも知れません（笑）

　母平均と母分散を設定してやると、それにしたがって正規乱数を生成するような、乱数発生器を考えます。もちろん、観測値はすべてiid（独立同分布）とします。
　分散について考える時、ふつうはこういう数直線上の分布をイメージすると思います。これは母平均2、母分散1に設定しているのですが、標本平均は少しズレていて2より少し大きくなっていますね。

　ここで少しデータの見方を変えて、乱数発生器から生成されたサンプル（データセット）を、「サンプルサイズNを次元数とする空間」上にマッピングすることにします。イメージしづらいかも知れませんが、データポイント1つが1本の座標軸に対応しているような、N次元の空間を設けてやると、N個のデータからなる「データセット」そのものを、その空間上の1点として表現できるわけです。
　クラスの生徒が40人いるなら、生徒1人1人が座標軸となって、それが40本直交している、40次元の超空間を考えます。軸には0から100までの目盛りが売ってあり、たとえば算数のテストを実施したら、そのテストの「40人分の成績」を、この40次元空間上の1点としてプロットできます。国語の成績、社会の成績なども、同様にプロットできますし、算数のテストを何回も実施したなら、それも全部プロットできる。
　40次元の空間を図示することはできないので、3人の生徒しかいないという前提で、3次元の空間に50回のテストの結果をプロットするとこんな感じになります。

　この場合、観測値が「独立にばらつく方向」が、サンプルサイズ分（つまり3本）だけ存在すると言え、この捉え方の下では、「分散」とはすなわち1方向あたりの平均的なばらつき具合ということになります。乱数実験であれば母平均が分かるので、1つ1つの座標軸にその母平均の印を付けておくと、プロットされた1つのデータセットが「それぞれの軸に関して母平均からどれぐらいズレているか」を測ったものが母分散に一致し、それを軸の数だけ足し合わせると、総二乗誤差になります。
　そして、3個の観測値がすべて等しくなる点を結んだ線（2次元平面の場合でいうy=xの線）をマクロ経済学にならって45度線と呼ぶことにすれば（対角線でもいいが）、45度線からの距離が、そのサンプルの二乗誤差を表すことになりますね。

　ばらつきの総和としての総二乗誤差（偏差平方和）は、ばらつきの方向の数に比例します。つまり「分散×サンプルサイズ」が、総二乗誤差の期待値となります。しかしこれはあくまで、「乱数発生器に設定した母平均からの二乗誤差」を考える場合の話ですね。
　一方、サンプリングのたびに得られる「そのサンプルの標本平均からの二乗誤差」を考えるとすると、その期待値は上述の値よりも必ず減ります。なぜなら、じつは標本平均というのは理論上、「そこから測った二乗誤差を最小にする点」だからです。ちなみに中央値は、二乗しない「偏差」（の絶対値の和）を最小にする点です。この、「標本平均から測ると二乗誤差が最小になる」という話は、上述の特殊な空間で考えるよりも、元の数直線で考えたほうがイメージしやすいと思います。

　では、具体的にどれだけ「減る」のかというと、ちょうど「方向1本分」だけ減ります。したがって、総二乗誤差を方向の数で割って「1方向あたりの二乗誤差」（すなわち分散）に直すという処理をする際に、割る数をnではなくn-1にしておかないと、期待値が本来の「1方向あたりの二乗誤差」（分散）に一致しなくなるわけです。

　減る量が「ちょうど1方向分」だけであるのはなぜなのか。標本平均というのは、「そこから測った誤差がn-1個分だけ分かれば、残り1個のデータは自動的に決まってしまうような点」であり、これはすなわち、「標本平均から測ると、ばらつきの方向が1本減る」ということに等しい。このことを指して「残差の自由度が1下がっている」と言っていることになります。
　実際には、観測そのものが1個減るわけではないので、変動の方向も減ってるわけではないのですが、全体としてちょうど1本分減るように、すべての観測値が平均に引き寄せられます。……いや、じつはよく考えるとこの言い方にも語弊があって、標本平均というのはそこからデータが生成される「母」（母平均）ではなく、「本当の母がどこにいるかはわからないので、サンプルが得られるたびに仮にそこに母がいることにしよう」と決めている点に過ぎないので、正確には、平均からの偏差がトータルで方向1本分減るように、平均のほうが観測値に近づいていくというほうが、イメージとしては正しいです。

　「n-1で割る理由」を直観的に理解したいのなら、ここまでの説明で十分ですが、もう少し続けます。
　上述の話は、元のN次元空間のまま考えて、標本平均から測った場合の総二乗誤差の期待値が、乱数発生器に設定した母平均から測った総二乗誤差の期待値に比べて、方向1本分だけ減っているということです。一方、何回もサンプリングを繰り返してそのたびに残差（偏差）のベクトルを記録していくとすると、もちろん残差ベクトル自体はN次元のベクトルになるのですが、平均が特定の値になるように拘束を与えると（もちろん実際の観測ではそんな拘束はできませんが）、そのN次元空間の中に収まるn-1次元の超平面上に分布するようになります。
　つまり、基底を取り直せばn-1次元のデータとして再解釈できるということで、これが、「残差の自由度が1下がる」ということの直観的なイメージに対応します。
　下の図は、3次元空間上に分布するはずのデータに、平均が特定の値になるような拘束を与えた場合のグラフです。

　実際の計測では、サンプリングのたびに標本平均自体も変わるので、同じ平面上に分布するようにはならないのですが、いわば1回1回のサンプリングで得られる計測値のセットが、「3次元で動くものを2次元に押し込める」のと同じ程度に、自由を失うということですね。（ここ、意味がわかりにくいと思いますが。）

　「自由度」って、統計の勉強をするときに「なんかそんな概念あるらしい」程度で済ませている人も多いと思うのですが（うちの学生はほとんどそうです）、「自由度あたりのばらつきを評価する」という形で登場する、けっこう重要な概念だということを理解しておく必要があります。「独立に変動する情報の個数」（自由度）が情報のばらつきを捉える上で重要だということは、考えれば分かるはずですけど、考えなくても分析の実務はできてしまうので放置されがちな気がします。

　作図のためのRのコードは以下のとおりです。

set.seed(1)

n <- 100  # 点の数
x <- rnorm(n, mean=2, sd=1)

col_pt <- rgb(0.1, 0.4, 0.9, alpha=0.40)  # 青で半透明

plot(x, rep(0, n),
     pch=16,
     col=col_pt,
     cex=1.0,  # 点の大きさ
     ylim=c(-0.1, 0.1),
     yaxt="n", ylab="",
     xlab="Value",
     main="Samples on a Number Line")

abline(h=0, col="gray")

# 平均
abline(v=mean(x), col="red", lwd=2)

library(plotly)
set.seed(1)

T <- 50          # テスト回数（点の数）
mu <- 60          # 平均点
sigma <- 10       # 点数のばらつき（標準偏差）

# 3人分の得点（各行が1回のテスト、各列が生徒）
X <- matrix(rnorm(T * 3, mean = mu, sd = sigma), ncol = 3)
colnames(X) <- c("Student1", "Student2", "Student3")

# 0～100に丸める
X <- pmin(pmax(X, 0), 100)

plot_ly() |>
  add_markers(
    x = X[,1], y = X[,2], z = X[,3],
    marker = list(
      size = 2,
      color = "rgba(30, 30, 200, 0.35)"  # 青の透明
    )
  ) |>
  layout(
    title = "Each test = one point in 3D (3 students)",
    scene = list(
      xaxis = list(title = "Student 1 score", range = c(0, 100)),
      yaxis = list(title = "Student 2 score", range = c(0, 100)),
      zaxis = list(title = "Student 3 score", range = c(0, 100)),
      aspectmode = "cube"
    )
  )

library(plotly)
set.seed(1)

T <- 50
mu <- 60
sigma <- 10

# データ生成
X <- matrix(rnorm(T * 3, mean = mu, sd = sigma), ncol = 3)
colnames(X) <- c("Student1", "Student2", "Student3")
X <- pmin(pmax(X, 0), 100)

# 対角線用データ（0〜100）
diag_vals <- seq(0, 100, length.out = 100)

plot_ly() |>
  add_markers(
    x = X[,1], y = X[,2], z = X[,3],
    marker = list(
      size = 2,
      color = "rgba(30, 30, 200, 0.35)"
    )
  ) |>
  add_trace(
    x = diag_vals,
    y = diag_vals,
    z = diag_vals,
    type = "scatter3d",
    mode = "lines",
    line = list(
      color = "black",
      width = 5
    ),
    name = "45-degree line (x=y=z)"
  ) |>
  layout(
    title = "Each test = one point in 3D (3 students)",
    scene = list(
      xaxis = list(title = "Student 1 score", range = c(0, 100)),
      yaxis = list(title = "Student 2 score", range = c(0, 100)),
      zaxis = list(title = "Student 3 score", range = c(0, 100)),
      aspectmode = "cube"
    )
  )

library(plotly)
set.seed(0)

N <- 1500
mu <- 2
sigma <- 1

X <- matrix(rnorm(N*3, mean=mu, sd=sigma), ncol=3)
colnames(X) <- c("X1","X2","X3")

R <- X - rowMeans(X)

grid <- seq(min(R), max(R), length.out=25)
z_mat <- outer(grid, grid, function(x,y) -(x+y))

# ふつうにプロット
p1 <- plot_ly() |>
  add_markers(x=X[,1], y=X[,2], z=X[,3],
              marker=list(size=1.5, opacity=0.35)) |>
  layout(title="Raw samples (R^3)",
         scene=list(
           xaxis=list(title="X1"),
           yaxis=list(title="X2"),
           zaxis=list(title="X3")
         ))

# 中心化した場合
p2 <- plot_ly() |>
  add_markers(x=R[,1], y=R[,2], z=R[,3],
              marker=list(size=1.5, opacity=0.5)) |>
  add_surface(x=grid, y=grid, z=z_mat,
              opacity=0.15, showscale=FALSE) |>
  layout(
    title="Centered residuals (R1+R2+R3=0)",
    scene=list(
      xaxis=list(title="R1"),
      yaxis=list(title="R2"),
      zaxis=list(title="R3"),
      aspectmode="cube"
    )
  )

p1  # 表示する
p2  # 表示する

（参考：このブログのおすすめ記事一覧はコチラ）