Drastic changes before the 2011 Tohoku earthquake, revealed by exploratory data analysis
https://arxiv.org/ftp/arxiv/papers/2302/2302.02326.pdf
図がここではあんまきれいに出ないので、こっちのPDFをみてください。
地震を予測することは、特にリスクの高い国々にとってきわめて重要であり、これまで多くの努力がなされてきたが、まだ実現には至っていない。それにもかかわらず、古い理論が検証ないままにただ信じられているほど、地震学における統計的アプローチは乏しい。そこで、日本における地震の時間とマグニチュードの記録を探索的データ解析（EDA）により分析した。EDAはデータの特徴に基づいたパラメトリックな統計手法であり，データ駆動型の調査に適している．各データセットの分布様式を決定し、重要なパラメータを見出した。これにより、データ中の異常値を特定し、評価することができるようになった。2011年の震災以前には、大きな群発地震がありえない頻度で発生していた。１年以上前から、すべての地震の発生頻度と規模が大きくなっていた；この2つの変化により、より大きな地震が発生しやすくなり、M9の地震でも2年に1回ほどの確率で発生すると予想された。このように、EDAは単純な計測から有用な情報を引き出すことに成功した。同様な手法で、もっと多くのデータの異常を検出し、評価することができれば、より正確な地震予知につながるだろう。
Introduction
日本はいくつかのプレートが衝突して沈み込む場所にあり地震の回数が多く、また強い地震がおおい(Maruyama 1994; Stern 2002; Bird 2003; Fujii & Kazuki 2008; Japan Meteorological Agency (JMA) 2023a)。その被害は甚大であり、大地震の予知が急務となっている。例えば、2011年の東北地方太平洋沖地震では、死者・行方不明者が2万2千人以上となった(Cabinet Office Japan 2023)。とある小学校では校庭に避難した子供たちと引率の先生が津波にのまれ誰も帰ってこなかった。ある原子力発電所は外部電源を喪失して爆発し、まだその撤去の方法さえわかっていない。放射能汚染で立ち入りできない地域が広範囲に広がっている。これらの悲劇は、ほんの数時間前にわかっていれば避けられたかもしれない。しかし、様々な要因を正確に測定しようとする多くの努力にもかかわらず、地震予知は実現されていない(Mogi 1998; Hayakawa 2012; Stark 2022)。この失敗の理由の一つとして、統計的アプローチの少なさが考えられる。統計学はランダムな事象の確率を正確に予測することができ、この特性を生かさなければ予知は困難である。
このような状況は、探索的データ解析（EDA）を用いることで解決できるかもしれない。これは、データ分布などのデータの特徴に基づいて、データから有効な情報を抽出するパラメトリックな統計手法である（Tukey 1977; Croarkin & Tobias 2012)。分布についての予見をもたずにデータからそれを調べる点で、古典的な統計学とは異なる。データがなぜその性質をもつのか、隠されたデータの構造を明らかにしつつ、データがもつ意味を最大限に抽出することを目的とする。データドリブンな分析を可能にする利点をもつ。具体例として私の経験をひとつ紹介させてほしい。さまざまな遺伝子の発現を測定するマイクロアレイという技術がある。得られるのは相対値であるために標準化を必要とする。このデータが対数正規分布するという性質を利用して標準化を行い(Konishi 2004)、異なるチップでもデータを比較可能にした(Konishi et al. 2008)。分布は普遍的なものだからひとつの単位系として使用できるからだ。またこの性質が、細胞がmRNAを調節する機構に由来することを明らかにし、細胞がどうそれを調節するのかを熱力学をつかって説明できるようにした(Konishi 2005)。
ここでは、日本の地震のデータセットにEDAを適用した。発見されたパラメータは大地震の前に変化し、より大きな地震が発生しやすいことを示していた。このような不安定な性質は、いま誤って信じられている地震の性質とはかなり異なっている。EDAを使えば、このような変化を簡単に検出し、評価することができるようになるはずだ。
Materials and methods
（ふつうの論文にくらべて教科書的な書き方になっている、読者に統計学の知識を期待できるかどうか不安だったからだ。エディターから削るように勧められたがそのまま残してある）。
気象庁が公開しているpdfファイルから、有感地震がおきた時刻とマグニチュードを抽出して分析した。これは彼らが測定している様々なデータのごく一部であるため、たとえば地域によってわけることなしにまとめて取り扱うことにする。
見つかった分布のひとつはThe Exponential Distributionである。これはそのdensity が x≥0 においてf(x) = λe−λx であらわされる、λを唯一のパラメータにする分布だ。左に大きくかたよっていて、大きな値ほど小さいdensityをもつ。最頻値はゼロで、平均とSDはどちらも 1/λ である。
ランダムにおきる事象の時間間隔を表す分布として知られている 地震がもしランダムにおきているのならば時間間隔はこれに従う しかしなんらかの理由でそのランダムネスの条件が変わればλが変化したり分布がまざったりする。
もうひとつは良く使われる正規分布だ。これはそのdensityがf(x) = (2πσ2)-1/2exp(-(x-μ)2/2σ2)であらわされる分布で、μ と σのふたつのパラメータをもつ。Independent and identically distributed random numbersの和がこの分布になることが知られている。xの対数値がこの分布に従う場合は対数正規分布という。
ロバスト推定
データに外れ値が含まれる場合、平均値やSDなどのパラメータの推定はその外れ値に影響される。これを回避する計算方法がいくつかあり、ロバスト法と呼ばれる。例えば、平均を求める方法として中央値やトリム平均があり、SDを求める方法として中央値絶対偏差(MAD)がある。また、R には line や lm などの関数があり、ロバストに直線近似直線を求めることができる。
Quantile-Quantile (QQ) plot
分布の確認にはヒストグラムを用いることが多いが，より厳密に確認するために，QQプロットを用いた．データ数nとパラメータを指定することで、分布の理想値を求めることができる。本来、QQプロットはデータセットの各分位のデータと理想値を比較するが、ここではソートされたデータ全体と理想値を比較した。データセットが分布モデルに従えば、直線関係が形成される。指数分布の場合、λ=1の理想値が用いられた。QQプロットの傾きは、データの1/λを表す。正規分布の場合、σ=1、μ=0の理想値が用いられた。QQプロットの傾きはσ、切片はμを示す。パラメータはRのline関数を用いたロバスト計算で同定した。
P値、累積分布における確率、分位数
P値はある値よりも極端な値が得られる確率のことである。統計的検定でよく使われるが、ここでは値を評価するために使う。より小さい値の場合は、P値は累積分布の確率であり、逆の場合は累積分布を1から引いた値となる（図S1）。P値が小さいほど、異常が強い（偶然に得られる可能性が低い）ことを意味する。zスコア＝（データ-μ）/σ＝-2であれば、累積確率はp値と同じで、Rの記法で書くならばpnorm（-2）＝0.023である。同様に、zスコアが2であれば、pnorm(2)=0.977、したがって、p値は1-0.977=0.023となる。ちなみに、Rで累積確率から分位値をもとめるには、qnorm(0.023)=-2とする。
平均値の分布
中心極限定理によれば、任意の乱数の平均は、σ = SD/n1/2 で正規分布すると予想され、ここで n は平均のデータ点数である (Stuart & Ord 1994; Bohm & Zech 2010)。この性質を利用すると、移動平均が年に一度だけ到達する水準を推定することができる。例えば、図3Cの本震前の有感地震回数は１年で1,246回であった。 これはRの表記法では、分位数はqnorm(1-1/1246)=3.15となる。5日間の平均測定回数は14.3回であった。大きさの中央値は3.38で、MADは0.276であった。ここから、その水準は中央値+ MAD/sqrt(14.3)*3.15 = 3.57と算出される。これは図3Cの緑色の線で示されている。
逆に、移動平均値がどの程度に極端なものであるかを、p値を推定することで評価することができる。図2Cの2011年1月12日に移動平均が4.0になったが、この時の平均のzスコアは(4.0-中央値)/(MAD/sqrt(14.3))=8.5であった。これはp値を推定するには大きすぎていて、そのR表記は1 - pnorm(8.5) = 0である。
計算は，EDA に適した統計計算とグラフィックスのためのフリーソフトウェア環境である R ソフトウェアパッケージ (R Core Team 2023) を用いて行った．インストールと計算のためのスクリプト一式はsupplementで提供されている。
Results
予想どおり時間間隔はThe Exponential Distributionに従った。1998年の10月はあまり大きな地震がおきていないちょっと珍しい月であった。ここで測定するとλは1/9.3であった。しかし、λは時々増加する（図1B）。過去20年間に3回極端に増加したが、いずれも余震が連続する大地震に対応したものである（図1C）。これらの地震は通常の状況のランダム性から外れているため、ランダム性の変化が生じたのである。このパネルでわかるように、これらは冪級数（Utsu et al 1995）に従って急激に減少している。また、これは本震が発生する前からデータ上に現れる。群発地震があるとQQplotでは原点で直線関係に歪みが生じ（パネルB）、これは時間間隔が短いほど顕著になる（図1D、1E）。あるいは、両軸の対数をとることで、原点付近の異常が検出されやすくなる（図1F）。対数軸のプロットでは、線形回帰直線とλの関係が変化する（図1の凡例を参照）。また、対数をとることで、着目するデータの範囲が若干変化するため、パラメータ推定も若干変化した。
なお、この変化は群発地震を除外しても生じていた、ほとんどのデータはλを大きくした直線上にあった（パネルB,D,E,F）。これは、2011年の震災前に、全ての地震の発生頻度が上昇していたことを示している。
Fig.1 図1 地震の時間間隔(時間)。(A)指数関数的な分布に厳密に従っている。λの逆数は時間間隔の平均値とSDである。パラメータλは安定していない。(B)2011年東北地方太平洋沖地震の1ヶ月前。λだけでなく、QQプロットの形状も原点付近で変化している。(C) λは時々大きく変化する。2000年の三宅島噴火、2011年の東日本大震災、2016年の熊本地震の3つのピークが関係している。(D, E) QQプロットの周期を短くすると変化が明確になる。3月11日の本震直前の分布。 (F) 原点付近の異常に着目する一つの方法として、両軸の対数（底＝10）を取る方法がある。このデータでは、高いλをもつ分布が混じると、傾きが急になる。直線近似の傾きは原理的に1である。1/λは10^切片となる。
タイムコースを図2に示す。青い線は測定5日前までの移動平均である。緑色の線は平均値より2σ低いレベルを示している。平均値が緑より低下したのは3回である。パネルAの”1” （青）では、2011年震災の震源地よりも内陸部で群発地震が発生した。ここでは、1ヶ月以上前から小さな地震が多数発生しており、それがオープンデータに記録されるほど強くなったものだ(JMA 2023b)。同時に、新潟でも群発地震が記録されており、ここでは3月12日にM5.9の地震が発生した。”2” は別の場所で発生したM5の地震とその余震、"3"は2011年の震災の直前に発生した群発地震である。パネルBは、2000年の三宅島の噴火の前の年のデータである。5つの群発地震が検出された。”1” と “3” は4/29に最高M4の地震が発生し、その後2001年4/25にM5.6の本震が発生した大分、”2” は噴火の2日前から始まった有珠山の火山性微動、"4 "は三宅島の噴火の前に起きた群発地震、”5” は噴火に伴う余震である。ちなみに、Aの平均間隔時間は3.7時間、Bは8.7時間で、2011年の震災前に頻度が高まっていることがわかる。
Fig.2 図2 地震の発生時間間隔の変化。(A)2011年東北地方太平洋沖地震発生前。薄赤色の点は東北の太平洋岸地域。青い線は5日移動平均。緑色の水平線は、本震前に推定された時間間隔の平均値より2σ短いレベルである。1〜3(青)で示した時間帯では、移動平均がこのレベルより低くなっている。本震後の小地震の記録は公表されていない（JMA 2023b）。 (B)2000年の三宅島噴火前のレベル。北海道の有珠山噴火もこの期間に発生している（左の灰色の縦線）。
マグニチュードの分布は正規分布に従った（Fig. 3A）、外れ値があったので、パラメータはロバストに求めている。また、2つの異なる傾きが観測される場合があった（図3B）。これは、日本にlocationとscaleが小さい下半分と、その逆の2つの状態が同時に存在する場合に起こる。実際、東北地方の2011年のデータでは、locationとscaleが大きいデータだけであった（パネルB、入れ子）。
マグニチュードの分布は、時に大きく変化する。地震が少なかった1998年１０月はscaleは0.59程度であった。ところが震災前には二倍以上に大きくなっている。これは、より大きな地震が起きやすくなっていることを示す。表1は、各マグニチュードのp値、つまり、地震発生ごとにそのマグニチュードより大きくなる確率と、その間隔から推定した年間の期待値を示したものである。状態1は1998年10月、状態2は震災前であり、状態1ではM9はありえないが、状態2では2年に1回発生すると予想される。実際、この年のマグニチュードのQQプロットを見ると、M9はさほどの異常値ではないことがわかる（パネルB）。しかしパネルAでは不可能な値である。
移動平均の大きさを評価することができる。例えば、図3Cの本震前の有感地震の回数は1,246回である。この数から、平均が年に1回だけ超える可能性のある分位数を推定することができる（Materials and Methods）。これがパネルCの3.57の緑線である。パネルDの値は3.16と低い。
驚くべきことに、2010-11年には移動平均がこれを何度も超えていたのである（図3C）。このような現象は、通常では考えられない。比較のために、2000年の三宅島の噴火の前の年（図3D）；噴火の前は、いくつかの異なる場所で発生した地震が重なった以外は、平均はほとんど緑の線を超えていないことがわかる。パネルCでは、東北地方の沿岸部（薄赤色）でM4-6クラスの地震が多く発生し、平均値を上げている。この地域では、この1年間に地震が活発化していて、月報でも何度かコメントされていた（JMA 2023b）。ちなみに、それらの移動平均がどの程度稀なものであるかは、p値という形で推定することができます（Materials and Methods）。2011年1月12日に移動平均が4.0になったが、この値のzスコアは8.5で、p値は0であった。つまり、このような大きな異常値は、通常ありえないということである。
特に重要なのは、本震の数日前に発生した地震（パネルC）で、本震の震源域の北側の海底で発生した地震です (Ide et al. 2011; Editorial Committee 2013; JMA 2023b)。減衰のない平均 M5.2 の地震が 43 回記録されている。このzスコアはなんと14である。これは群発地震が独立した事象ではなく、互いに関連していることに起因するアーチファクトであるはずだ。このことは、多くのアスペリティが互いに支え合うことで、近くにある他の多くの破壊を誘発することを強く示している。周辺に蓄積された弾性エネルギーが不安定化していたはずだ。
Fig.3 図3 マグニチュードの分布と時間経過。(A)マグニチュードの分布。正規分布の理想的な分布と比較している。パラメータが変化する傾向がある。(B)東日本大震災前。この場合、2つの状態が混在しているように見える。入れ子は、2011年1月1日から3月11日までの東北地方太平洋沖地震。(C)2011年東北地方太平洋沖地震前の時間経過。緑は青の移動平均が1年に1回だけ到達すると予想されるレベル。(D)2000年三宅島噴火前。
Table 1
表1 観測された分布から予測されるP値と年間地震数の期待値。
ちなみにこのマグニチュードが正規分布するという知見は、ながく信じられてきたGutenbert-Richter law、つまり地震のエネルギーの対数値と、それがおきる頻度の対数値が比例するという理論である、とは一致しない。Oct2015-Mar2016は熊本地震の前の期間で、この間はあまり大きな地震がおきていない。ここをNormalQQplotでみるときれいな直線が得られるFig.4A。モデルの比較を公平にするために、この法則の表現を更新した（図S2）。図4は2つのQQプロットで、パネルAは正規分布、パネルBはGutenberg-Richter則である。正規分布モデルの方がより広い範囲のデータをカバーしていることは明らかである。
Fig. 4 図4 Gutenberg-Richter則との比較。(A)大きさの正規のQQプロット、データの全範囲で直線関係が得られている。(B）別の表現でのGutenberg-Richterの法則。指数分布のQQプロット。パネルAとの比較のために表示を変えている。詳細は図S2の注釈を参照。
残念ながら、いつも必ずしも明確な地震の兆候があるわけではない。例えば、2016年の熊本地震の前には、マグニチュードの変化はわずかであり（図1C、4A）、マグニチュードの移動平均には異常は見られず、直前に群発地震が発生しただけだった。特に火山性地震では、2000年の三宅島や有珠山（図3D）、熊本地震に伴うであろう2016年の阿蘇山、2021年のトカラ列島などがそうだった。
Discussion
現象に適した分布モデルを知ることは、ランダム性が持続する限り何が起こるかを推定でき、ランダム性の変化をいち早く検出できることを意味する。EDAは、この現象に対応するための理想的なアプローチだった。Tukeyの時代は計算が非常に難しく、彼の著書を見ると、いかにして計算を簡単にするかということに多くの努力が払われていることがわかる（Tukey 1977）。しかし、現在では、Rを使って誰でも簡単に実行することができます(Croarkin & Tobias 2012; R Core Team 2023)。もう一つの利点は、扱う確率事象のp値がわかるので、何かおかしなことが起きたときに、そのp値と期待値を推定できることである。これを利用すれば、それぞれの事象を客観的に評価することができる。避難勧告を出す場合、それがどの程度の経済的損失をもたらすかを考えなければならない。そのためには、今起きている現象がどの程度異常なのかを把握する必要がある。
ここでは、オープンデータを使用したが、これは観測データ全体のごく一部に過ぎない。もしもっと多くのデータがあれば、例えば10kmごとのメッシュを解析して、おそらく地震の発生しにくい空白域を明らかにすることができただろう。そのような空白域は、スティックスリップ現象（Scholz 1998; Matsukawa 2009）のピン止め領域かもしれないし、よく滑る領域かもしれない。前者であれば、弾性エネルギーの貯蔵庫となるはずで、次の大地震の震源地となる可能性がある。また、他の地域からの情報は基本的にノイズであるため、ノイズレベルを大きく下げることになる。そのため、異常の発見がより敏感になる。
群発地震は、間隔のランダム性が変化することで識別できる。余震は通常、急速に弱まる（図1C）（Utsu et al.1995）。おそらく、本震で放出されたエネルギーの一部が一時的に蓄積され、それが放散するときに発生するのだろう。しかし、孤立した群発地震はこのような減衰をしない。これは、別々にエネルギーを蓄えながらも、互いに支え合うアスペリティが連続して破裂するためである。この場合、広い範囲で発生する可能性があり、実際、2011年震災の”3”の場合（図3C）、その範囲は約50kmに及んだ (Ide et al. 2011; Editorial Committee 2013; JMA 2023b)。これは、プレートの面積が移動し、条件が変化したことを示しています。これは，より大きなエネルギーの貯蔵庫を不安定にすることになる．ここで確認された事例では，多くの孤立型群発地震のあと1年以内に本震がおきている。特に、群発地震のマグニチュードが大きい場合は、大きな本震が近づいていることを示している。震災では、本震とそれに続く長さ100km以上の大破壊が連鎖的に発生し、広範囲の海底変動を引き起こし、津波を引き起こした (Ide et al. 2011; JMA 2023b)。
マグニチュードが概略で正規分布するということは,地震のエネルギーがどう蓄積されるかということに関連している。地震が放出するエネルギーは、大陸を運んでいるプレート移動のエネルギーのほんの一部にすぎない。プレートがするであろうStick-slip phenomenonで(Scholz 1998; Matsukawa 2009)、あるピンしているasperityがどの程度に移動エネルギーを弾性エネルギーとして蓄えるかは、そのasperity (Lay & Kanamori 1980)のいくつかの性質に依存するだろう、たとえば形状、大きさ、かたさ、温度、速度。またひとつのasperityが崩壊するまでの時間はexponential distributionのようにして分布するだろう、これはプレートの硬さとも関係するかもしれない。これらファクターのいくつかが掛け算的に影響しあってエネルギーを決定するなら、それは対数正規分布するはずだ。おそらく私たちがみているのはこの現象だ、マグニチュードは1000^.5を底とする対数値だからだ(JMA 2020)。この正規性はエネルギー蓄積のモデルのヒントになると思われる。
時間間隔と大きさの分布はともに変化し続けている。これは、プレートが不均質な物体であり、接触面で硬度が変化すると、パラメータが変化するのだろう。実際、2011年東北地方太平洋沖地震の前年には、マグニチュードのσと間隔のλが非常に大きくなっていた（図3B）。このような条件下では、大きなマグニチュードの地震が発生しやすくなる（表1）。これは長期間続く可能性がある（図1B、3B）。
古典的なGutenberg-Richter則（Gutenberg & Richter 1944）は今でも信頼されているが（Fujii & Kazuki 2008; Stark 2022）、有感地震のデータとは一致しなかった（図4）。この法則は、M3.5 未満の測定値がまだ正確でなかった頃に考案されたものである。データとの整合性の確認は原始的であった（図 S2）。また、この式が物理的に何を意味するのかも不明であり、現象の理解を深めるものではなかった。これらは、歴史的な背景を考えると、やむを得ない欠点である。欠測した小地震が、この法則が予測するほどに多くある場合にだけ、この法則はまだ適切だと言えるだろう。しかしこれは現時点では検証不可能な仮定であるため、法則を維持するためにこれを導入することはできない。知覚可能なデータ（図3A、図3B、図4A）では正規性はしっかりとしており、マグニチュードの分布を扱うのにうまく機能している。
80年間もの間、この法則の検証を怠ったことに対して、歴史的背景は言い訳にならない。まるで中世のようなこのエピソードは、この分野がいかに統計学を欠いてきたかの証拠だろう。特に、線形回帰の傾きであるパラメータb（図S2A）がこの法則では常に一定であるという認識は、観測された事実 (Hirose et al. 2002)（図3AB、表1）と矛盾しており、修正されなければならない。ここで挙げた例から明らかなように、実際には大地震の可能性は常に変化しており、これらは統計的アプローチによって容易に推定できる。
ここで、2011年の震災の前に、3つの異常が検出された。まず、地震発生間隔が短くなったこと（図1B）。これは群発地震を伴うことが多かったが、群発地震が短くなった原因の全てではない。むしろ、間隔が短くなった状況が群発地震の引き金になった可能性もある。第二に、これらのパラメーターの大きさが大きくなったことである（図3B）。1つ目と2つ目の変化により、大地震の予想が拡大した（表1）。3つ目は、マグニチュードの移動平均がありえないほど大きくなっていたことだ（図3C）。特に、本震の直前に大きな群発地震が発生していた。残念ながら、この3つが常に全ての地震で前兆として現れるわけではない。これは、これらがプレートの沈み込みに伴う地震の特徴であるためと思われる。おそらく、異なるメカニズムで起こる地震を予知するためには、特性の異なる他の測定値を用いるべきだろう。しかし、少なくとも、孤立した群発地震に対しては、注意を喚起する必要がある。同時に重なる異常が観測された場合は、さらに警戒レベルを上げる必要がある。ただし、過去のデータを分析するという本研究のアプローチは、予測という点ではやや不公平かもしれない。ここで明らかになった考え方の信憑性は、今後、研究者が統計的アプローチによって何を予測するかで判断すべきだろう。
References
1 Maruyama, S. Plume tectonics. Journal - Geological Survey of Japan 100, 24-49, doi:10.5575/geosoc.100.24 (1994).
2 Stern, R. J. Subduction zones. Reviews of Geophysics 40, 3-1-3-38, doi:https://doi.org/10.1029/2001RG000108 (2002).
3 Bird, P. An updated digital model of plate boundaries. Geochemistry, Geophysics, Geosystems 4, doi:https://doi.org/10.1029/2001GC000252 (2003).
4 Fujii, T. & Kazuki, K. Encyclopedia of Earthquakes, Tsunamis and Volcanoes (Japanese). (Maruzen, 2008).
5 Japan Meteorological Agency. How Earthquakes Occur
6 Japan Cabinet Office. Special Feature: The Great East Japan Earthquake (Japanese),
7 Mogi, K. Thinking about Earthquake Prediction (Japanese). (Iwanami, 1998). ISBN: 978-4004305958
8 Hayakawa, M. The frontier of earchquake prediction studies (Kinokuniya, 2012). ISBN 9784931507166
9 Stark, P. B. Pay no attention to the model behind the curtain. Pure and Applied Geophysics 179, 4121-4145, doi:10.1007/s00024-022-03137-2 (2022).
10 Tukey, J. W. Exploratory data analysis. (Reading, Mass. : Addison-Wesley Pub. Co., 1977). ISBN: 0-201-07616-0.
11 Croarkin, C. & Tobias, P. NIST/SEMATECH e-Handbook of Statistical Methods,
12 Konishi, T. Three-parameter lognormal distribution ubiquitously found in cDNA microarray data and its application to parametric data treatment. BMC Bioinformatics 5, 5, doi:10.1186/1471-2105-5-5 (2004).
13 Konishi, T. et al. Coincidence between transcriptome analyses on different microarray platforms using a parametric framework. PLoS ONE 3, e3555, doi:https://doi.org/10.1371/journal.pone.0003555 (2008).
14 Konishi, T. A thermodynamic model of transcriptome formation. Nucleic acids research 33, 6587-6592, doi:10.1093/nar/gki967 (2005).
15 Utsu, T., Ogata, Y. & Matsu'ura, R. S. The Centenary of the Omori Formula for a decay law of aftershock activity. Journal of Physics of the Earth 43, 1-33 (1995).
16 Stuart, A. & Ord, J. K. Distribution theory. (Hodder Arnold, 1994). ISBN: 0340614307
17 Bohm, G. & Zech, G. Introduction to statistics and data analysis for physicists. (2010).
18 Japan Meteorological Agency. Summary of seismic activity for each month,
19 Ide, S., Baltay, A. & Beroza, G. C. Shallow dynamic overshoot and energetic deep rupture in the 2011 Mw 9.0 Tohoku-Oki Earthquake. Science 332, 1426-1429, doi:doi:10.1126/science.1207020 (2011).
20 Editorial Committee for the joint investigation report on the great east Japan earthquake. Fundamental Aspects 1., in report on the great east Japan earthquake disaster, The Japan Society of Earthquake Engineering (2013), ISBN 978-4-907517-00-7 (2013).
21 Gutenberg, B. & Richter, C. F. Frequency of earthquakes in California. Bulletin of the Seismological Society of America 34, 185-188 (1944).
22 Scholz, C. H. Earthquakes and friction laws. Nature 391, 37-42, doi:10.1038/34097 (1998).
23 Matsukawa, H. Friction of macroscopic systems − from paper and rock to eatrthquake − (Japanese). Hyomen Kagaku 30, 548-553, doi: https://doi.org/10.1380/jsssj.30.548 (2009).
24 Lay, T. & Kanamori, H. Earthquake doublets in the Solomon Islands. Physics of the Earth and Planetary Interiors 21, 283-304, doi:10.1016/0031-9201(80)90134-X (1980).
25 Japan Meteorological Agency. Technical reference material on the outline and processing method of Earthquake Early Warning (Japanese),
26 Hirose, F., Nakamura, A. & Hasegawa, A. Changes in b-values associated with the destruction of asperities (Japanese). Jishin, 2nd series 55, 249-260, doi:10.4294/zisin1948.55.3_249 (2002).
27 R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2023).
Supplementary Information
Fig. S1 and S2の他、 “data&script.zip”は次の内容である
これらは figshareからダウンロードできる。 doi: 10.6084/m9.figshare.22010279
Rscript.txt
Rコードとインストールの解説
data2202.txt
データ、タブ区切りテキスト、Feb 2022.
data2203.txt
データ、カンマ区切りテキスト, Mar 2022. 漢字データを含む。
.RData
Rの空データ。
ダブルクリックでディレクトリを指定しつつRを立ち上げられる。
Fig. S1 図S1 正規分布の累積分布関数、μ＝0、σ＝1。分位数が2のとき、これより大きい確率は1-pnorm(2)=0.023で、点線の横線より上の部分である。
Fig. S2 図S2 Gutenberg-Richterの法則。(A)頻度の対数値とマグニチュードの関係。直線はM3.5以上のデータを直線近似したもの（a=4.0, b=0.65）。これがGutenberg-Richter則の本来の表現であるが、この式の表現は理想的とは言い難い。データを大きさで分類する必要があるため、まとめたデータが小さくなり、正確な検証ができないのである。この法則をn = 10^(a - bM)と書くと、nは各クラスの頻度、aとbは定数で、a = 0のときλ = b/log(e)で指数分布するはずで、実際に指数関数的に分布する乱数はこの法則を満たす；（B）データシミュレーション、指数分布する乱数を本来のGutenberg-Richter法則にあてはめた（rateは図4Bから推測したものである）。大雑把なものであるが、必要十分条件を満たすだろう。以上から、この法則は「マグニチュードは原点にずれがある指数分布に従う」と表現できる。これを図4Bに用いた。マグニチュードのゼロは測定感度の限界であるので、原理的にはaはゼロであるはずだ。このずれaは、得られなかったデータの量に起因していると考えられる。ここではデータポイント数は201であるが、理想的な数値をn=3.1E4とし、比較のために最大の201を用いるとaはキャンセルされる。したがって、(3.1E4 - 201)が大きさの小さい欠測データの数ではないか。
