正規分布の母分散をTensorFlowを使って最尤推定で数値計算してみる

統計

正規分布に従う標本から母集団の分散を最尤推定で計算してみました。最尤推定での計算はTensorFlowで勾配降下法を使いました。

最尤推定で計算した分散は不偏分散ではなく標本の分散と一致することが確認できました。

手順

  1. 正規分布に従う乱数を100個作成
  2. その平均と分散と不偏分散を計算
  3. 母集団が正規分布に従うと仮定したときの母集団の平均と分散をパラメータにとる尤度関数を設定
  4. 平均と分散を勾配降下法で変化させて、尤度関数が最大となる平均と分散を探す
  5. 尤度関数が最大となる平均と分散が2で計算した値と一致するかどうかを見る

Pythonコード

Google Colaboratoryで実行しました。

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import math
import tensorflow as tf

# 標本を作成
n = 100
x = range(n)
y = np.random.randn(n)

# 標本の平均と分散を計算
desc = stats.describe(y)
mean = desc.mean
uvariance = desc.variance
variance = uvariance * (n-1) / n

# 尤度を計算する関数
logpi2 = 0.5 * math.log(2.0 * math.pi)
def prob():
  d = y - l_mean
  d2 = -0.5 * d * d / l_variance
  return tf.math.reduce_sum(d2) / n - 0.5 * tf.math.log(l_variance) - logpi2

# 目的となる損失関数
def loss():
  return -prob()

maxLoopCount = 30

# 適当な初期値
mean_ini = 0.0
variance_ini = 1.0

# 複数のOptimizerを試す
opts = [
  tf.optimizers.Adam(learning_rate=0.2),
  CustomOptimizer(learning_rate=0.1),
]

for i in range(len(opts)):
  # 最初に正解を表示
  print((mean, variance, uvariance))

  # グラフにするための配列
  lossHistory = []

  l_mean = tf.Variable(mean_ini)
  l_variance = tf.Variable(variance_ini)

  for loopCount in range(maxLoopCount):
    l = float(loss())
    # グラフにするために記録
    lossHistory.append(l)
    # 最適化
    opts[i].minimize(loss, var_list = [l_mean, l_variance])
    # 学習の推移を数値で表示
    if loopCount % 10 == 0:
      print((float(l_mean), float(l_variance)))

  # グラフ化
  plt.rcParams['figure.figsize'] = (6.4, 4.8)
  plt.plot(range(maxLoopCount), lossHistory)
  plt.show()

CustomOptimizer はTensorFlowのOptimizerを自作するを参照。

結果

実行結果例

OptimizerがCustomOptimizerの場合

コードの実行結果です。1行目は100個の乱数の平均・分散・不偏分散です。4行目は最尤推定の平均、分散です。最尤推定の分散は乱数の不偏分散ではなく通常の分散と一致していることがわかります。

(0.1080843015362518, 1.0207865528862736, 1.031097528167953) # 平均、分散、不偏分散
(0.010808431543409824, 1.0016233921051025)
(0.10808192193508148, 1.0207873582839966)
(0.1080842986702919, 1.0207866430282593) # 30回ループ後の最尤推定の平均、分散

横軸はループ回数、縦軸は損失関数の値

f:id:suzuki-navi:20200926144719p:plain

OptimizerがAdamの場合

Adamですと収束が不十分でよくわかりません。

(0.1080843015362518, 1.0207865528862736, 1.031097528167953) # 平均、分散、不偏分散
(0.19999417662620544, 1.1999610662460327)
(0.15666580200195312, 1.0132553577423096)
(0.12916100025177002, 1.0091229677200317) # 30回ループ後の最尤推定の平均、分散

横軸はループ回数、縦軸は損失関数の値

f:id:suzuki-navi:20200926144734p:plain

わかること

わかることは2つ

リンク

関連する私の記事

fluentdでS3にJSONLでアップロード

AWS

fluentdでログをS3にJSONL(JSON Lines)でアップロードするところまでを試しました。

前提

  • fluentdはインストール済み
  • credentialsファイルなどでfluentd動作環境からS3アクセスする権限がある

プラグインインストール

fluentdのS3プラグインを別途インストールする必要があります。gemであれば以下のコマンドです。

$ gem install -N fluent-plugin-s3

fluentd設定ファイル

sample.conf という名前で以下のような設定ファイルを作成しました。

<source>
  # ファイルからログを入力
  @type tail
  path input.txt # 入力のファイル名
  pos_file input.pos
  tag test
  <parse>
    @type none
  </parse>
</source>

# fluentdが処理した日時をtimestampという項目名で保存
<filter test.**>
  @type record_transformer
  enable_ruby
  <record>
    timestamp ${Time.now.iso8601}
  </record>
</filter>

<match test.**>
  # 検証しやすいようにログをS3だけでなく標準出力にも出力
  @type copy
  <store>
    @type stdout
  </store>

  <store>
    @type s3

    # 設定詳細はプラグインのREADMEに書いてある
    # https://github.com/fluent/fluent-plugin-s3

    # S3バケット名
    s3_bucket xxxx

    # S3バケットのリージョン
    s3_region ap-northeast-1

    # S3バケット内のパスプレフィックス
    # ディレクトリのつもりならスラッシュを付ける
    path fluentd-sample-log/

    # S3のファイル存在有無をチェックするかどうか
    # デフォルトはtrue
    # falseにすれば余計なS3 APIコールを節約できる
    # そのかわりS3のキーが常にユニークになるように秒をキーに含めるなどの対応が必要
    check_object false

    # S3のキーのルール
    # デフォルトは %{path}%{time_slice}_%{index}.%{file_extension}
    # fluendの動作環境のタイムゾーン設定によらずUTCになる
    # ↓この設定でこういうキーになる: fluentd-sample-log/2020-09-25-084337.gz
    s3_object_key_format %{path}%Y-%m-%d-%{hms_slice}.%{file_extension}

    # AWS credentialsファイルの特定のprofileで実行したいならば
    <shared_credentials>
      profile_name xxxx
    </shared_credentials>

    # 検証しやすいようにログ3件ごとにS3出力
    # ログ3件ごとにファイルが分割される
    <buffer>
      chunk_limit_records 3
    </buffer>

    # JSONLにするための設定
    # これがないと、日時、タグ、JSONのタブ区切りファイルになる
    <format>
      @type json
    </format>
  </store>
</match>

check_object と s3_object_key_format については以下の記事が参考になります。

S3のキーの生成ロジックはソースコードでいうと以下のあたりが参考になります。

https://github.com/fluent/fluent-plugin-s3/blob/v1.4.0/lib/fluent/plugin/out_s3.rb#L276

動かしてみる

入力ファイルを作成し、fluentdを起動します。

$ touch input.txt
$ fluentd -c sample.conf

別のコンソールからテストログを送ります。

$ (echo hello1; echo hello2; echo hello3) >> input.txt

S3には以下のようなファイルができました。

$ aws s3 cp s3://xxxx/fluentd-sample-log/2020-09-25-091258.gz - | gunzip
{"message":"hello1","timestamp":"2020-09-25T18:12:58+09:00"}
{"message":"hello2","timestamp":"2020-09-25T18:12:58+09:00"}
{"message":"hello3","timestamp":"2020-09-25T18:12:58+09:00"}

なお、formatをJSONLにする設定を書かないと、次のような日時、タグ、JSONの順のタブ区切りファイルになりました。（↓testの前後はタブです）

2020-09-25T22:12:52+09:00 test {"message":"hello1","timestamp":"2020-09-25T22:12:52+09:00"}
2020-09-25T22:12:52+09:00 test {"message":"hello2","timestamp":"2020-09-25T22:12:52+09:00"}
2020-09-25T22:12:52+09:00 test {"message":"hello3","timestamp":"2020-09-25T22:12:52+09:00"}

バージョン情報

fluent-plugin-s3 (1.4.0)
fluentd (1.11.1)

リンク

fluentdに関する私の記事

手書き数字(MNIST)をTensorFlowの単純なニューラルネットワークで判別

機械学習

TensorFlowのディープラーニングの練習として、手書き数字(MNIST)をTensorFlowの単純なニューラルネットワークで判別させました。

Google Colaboratoryで実行しました。

TensorFlowのバージョンは2.3.0です。

!pip list | grep tensorflow
tensorflow                    2.3.0          
tensorflow-addons             0.8.3          
tensorflow-datasets           2.1.0          
tensorflow-estimator          2.3.0          
tensorflow-gcs-config         2.3.0          
tensorflow-hub                0.9.0          
tensorflow-metadata           0.24.0         
tensorflow-privacy            0.2.2          
tensorflow-probability        0.11.0

以下はGoogle Colaboratoryで実行したコードです。

まずはインポート。

import numpy as np
from matplotlib import pyplot as plt
import tensorflow as tf
from tensorflow.keras.datasets import mnist

出力画像の大きさを指定。デフォルトでは小さいため。

plt.rcParams['figure.figsize'] = (16.0, 7.0)

MNISTのデータをロード。

(x_train, y_train), (x_test, y_test) = mnist.load_data()

Downloading data from https://storage.googleapis.com/tensorflow/tf-keras-datasets/mnist.npz
11493376/11490434 [==============================] - 0s 0us/step

教師データの画像の一部を表示してみる。

for i in range(32):
  plt.subplot(5, 16, i + 1)
  plt.imshow(x_train[i], cmap='gray')

f:id:suzuki-navi:20200923204325p:plain

単純なニューラルネットワークで学習。

# 入力と出力サイズ
in_size = 28 * 28
out_size = 10

# 入力データを 28*28 の2次元配列から 784 の1次元配列に変換
x_train_reshape = x_train.reshape(-1, in_size).astype('float32') / 255
x_test_reshape = x_test.reshape(-1, in_size).astype('float32') / 255

# one-hot encoding
y_train_onehot = tf.keras.backend.one_hot(y_train, out_size)
y_test_onehot = tf.keras.backend.one_hot(y_test, out_size)

# モデル構造を定義
def createModel():
  hidden_size = 64
  model = tf.keras.models.Sequential()
  model.add(tf.keras.layers.Dense(hidden_size, activation='relu', input_shape=(in_size,)))
  model.add(tf.keras.layers.Dense(out_size, activation='softmax'))
  return model

model = createModel()

# モデルを構築
model.compile(
    loss = "categorical_crossentropy",
    optimizer = "adam",
    metrics=["accuracy"])

# 学習を実行
result = model.fit(x_train_reshape, y_train_onehot,
    batch_size=50,
    epochs=10,
    verbose=1,
    validation_data=(x_test_reshape, y_test_onehot))

# 学習の様子をグラフへ描画
def plotLearning():
  # ロスの推移をプロット
  plt.plot(result.history['loss'])
  plt.plot(result.history['val_loss'])
  plt.title('Loss')
  plt.legend(['train', 'test'], loc='upper left')
  plt.show()

  # 正解率の推移をプロット
  plt.plot(result.history['accuracy'])
  plt.plot(result.history['val_accuracy'])
  plt.title('Accuracy')
  plt.legend(['train', 'test'], loc='upper left')
  plt.show()

plotLearning()

Epoch 1/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.3328 - accuracy: 0.9073 - val_loss: 0.1855 - val_accuracy: 0.9454
Epoch 2/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.1622 - accuracy: 0.9529 - val_loss: 0.1375 - val_accuracy: 0.9593
Epoch 3/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.1192 - accuracy: 0.9649 - val_loss: 0.1159 - val_accuracy: 0.9668
Epoch 4/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0936 - accuracy: 0.9723 - val_loss: 0.0971 - val_accuracy: 0.9718
Epoch 5/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0761 - accuracy: 0.9773 - val_loss: 0.0950 - val_accuracy: 0.9711
Epoch 6/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0647 - accuracy: 0.9806 - val_loss: 0.0878 - val_accuracy: 0.9729
Epoch 7/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0552 - accuracy: 0.9834 - val_loss: 0.0947 - val_accuracy: 0.9714
Epoch 8/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0484 - accuracy: 0.9858 - val_loss: 0.0822 - val_accuracy: 0.9759
Epoch 9/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0423 - accuracy: 0.9874 - val_loss: 0.0825 - val_accuracy: 0.9747
Epoch 10/10
1200/1200 [==============================] - 2s 2ms/step - loss: 0.0362 - accuracy: 0.9891 - val_loss: 0.0837 - val_accuracy: 0.9738

f:id:suzuki-navi:20200923204346p:plain

f:id:suzuki-navi:20200923204400p:plain

テストデータを判別させてみる。（赤字が判別結果。たまに誤判定があるのがわかる）

pred = model.predict(x_test_reshape[0:64])
for i in range(64):
  plt.subplot(5, 16, i + 1)
  plt.imshow(x_test[i], cmap="gray")
  plt.text(0, -2, str(np.argmax(pred[i])), size="xx-large", color="#FF0000")

f:id:suzuki-navi:20200923204421p:plain

ニューラルネットワークの構造を画像で表示。

tf.keras.utils.plot_model(model, show_shapes=True)

f:id:suzuki-navi:20200923204440p:plain

TensorFlowのOptimizerを自作する

機械学習

TensorFlowのOptimizerを自分で独自に実装することができるようなので、実装してみたソースコードをメモとして残しておきます。以下の私の記事で使ったアルゴリズムを実装してみました。これで自分のアルゴリズムをTensorFlowで動かせます。

以下の記事ではこのOptimizerを実際に使っています。

これらの記事のPythonのソースコード中に CustomOptimizer というのがありますが、全部以下のソースコードで定義しています。

import numpy as np
import math
import tensorflow as tf

# オリジナルのOptimizer
class CustomOptimizer(tf.optimizers.Optimizer):
  def __init__(self, learning_rate=0.1, name='CustomOptimizer', **kwargs):
    super(CustomOptimizer, self).__init__(name, **kwargs)
    self.learning_rate = kwargs.get('learning_rate', learning_rate)

  def get_config(self):
    config = super(CustomOptimizer, self).get_config()
    config.update({
        'learning_rate': self.learning_rate
    })
    return config
  
  def _create_slots(self, var_list):
    for v in var_list:
      self.add_slot(v, 'grad', tf.zeros_like(v))
      self.add_slot(v, 'diff', tf.zeros_like(v))

  def _resource_apply_dense(self, grad, var):
    k = math.log(3.0)
    prev_grad = self.get_slot(var, 'grad').numpy()
    prev_diff = self.get_slot(var, 'diff').numpy()
    grad = grad.numpy()
    diff = np.copy(grad * self.learning_rate)
    replace_mask = (prev_grad != 0.0)
    diff[replace_mask] = prev_diff[replace_mask] * (0.5 + np.tanh(k * (1.5 * grad[replace_mask] / prev_grad[replace_mask] - 0.5)))
    self.get_slot(var, 'grad').assign(grad)
    self.get_slot(var, 'diff').assign(diff)
    return var.assign_sub(diff)

Google Colaboratoryで動かしていますが、TensorFlowのバージョンは2.3.0です。

Raspberry Piにキーボード・マウス・モニタなしでOSセットアップ

ちょっとやりたいことを思いついてRaspberry Piをセットアップしました。Raspberry Piを触るのは、たぶんですが2016年2月ごろ以来です。

キーボード・マウス・外部モニタのないRaspberry PiにOSをインストールしてSSH接続できるまでを設定します。

OSイメージをSDカードに書き込み

SDカードにOSのイメージを書き込みます。いまはRaspberry Pi Imagerというツールで書き込むのがいいようです。WindowsやMacで動くアプリですので、これを使いました。

「圧倒的に速い」──ラズパイにOSをインストールする新ツール「Raspberry Pi Imager」- ITmedia NEWS

このツールを起動すると書き込むOSを選択できますので、Raspberry Pi OS (32-bit) Liteを選択しました。Liteは2GBの小さいSDカードにも入ります。Liteでないほうは4GB必要です。手元にあったSDカードが2GBしかなかったので、Liteにしました。

このツールを使って書き込みが終わるとPCからアンマウントされるようです。Windows PCで作業したのですが、SDカードを挿し直すとWindowsのExplorerからSDカードの中が見れました。

この中に ssh という拡張子なしの3文字のファイル名で空ファイルを作成します。理由は以下のページに書いてありますが、Raspberry PiでSSHサーバを自動起動させるためです。

SSH (Secure Shell) - Raspberry Pi Documentation

WindowsのExplorerからは次のように見えます。

f:id:suzuki-navi:20200920220656p:plain

自宅LAN内のIP割り当てを確認

このあと起動するRaspberry PiはIPアドレスがわからないとSSH接続できません。起動する前に自宅LANに接続しているIP一覧を見ておきます。Raspberry Pi起動後にも同じものを見て、増えているIPがRaspberry Piということです。

arp -a というコマンドで見れるようです。このコマンドはMacでもWindowsコマンドプロンプトでも動きました。

$ arp -a

f:id:suzuki-navi:20200920220627p:plain

Raspberry Pi起動

SDカードをRaspberry Piに入れて、LANケーブルと電源を接続します。

しばらくしたから上記の arp -a をもう一度実行して、増えたIPを確認します。

IPがわかったらssh接続します。

$ ssh pi@192.168.XXX.XXX

ユーザ名は pi、パスワードは raspberry です。以下のページに書いてあります。

Installing operating system images - Raspberry Pi Documentation

これでRaspberry PiにSSH接続できます。

pi@raspberrypi:~ $ ls -al
total 28
drwxr-xr-x 3 pi   pi   4096 Sep 20 11:03 .
drwxr-xr-x 3 root root 4096 Aug 20 11:31 ..
-rw------- 1 pi   pi      7 Sep 20 11:03 .bash_history
-rw-r--r-- 1 pi   pi    220 Aug 20 11:31 .bash_logout
-rw-r--r-- 1 pi   pi   3523 Aug 20 11:31 .bashrc
drwx------ 3 pi   pi   4096 Sep 20 11:00 .gnupg
-rw-r--r-- 1 pi   pi    807 Aug 20 11:31 .profile

pi@raspberrypi:~ $ cat /etc/os-release 
PRETTY_NAME="Raspbian GNU/Linux 10 (buster)"
NAME="Raspbian GNU/Linux"
VERSION_ID="10"
VERSION="10 (buster)"
VERSION_CODENAME=buster
ID=raspbian
ID_LIKE=debian
HOME_URL="http://www.raspbian.org/"
SUPPORT_URL="http://www.raspbian.org/RaspbianForums"
BUG_REPORT_URL="http://www.raspbian.org/RaspbianBugs"

PostgreSQLの各テーブルのレコード数とデータ容量の概数を調べるには

PostgreSQLのレコード数を見るには SELECT COUNT(*) ... すればわかりますが、大きいテーブルですと時間がかかります。またディスク上のデータサイズはこれではわかりません。概数でいいので、レコード数とデータサイズを見るにはどうしたらよいか？

システムカタログにテーブルごとのレコード数が保存されています。データサイズはブロック数として保存されています。ブロック数にブロックのサイズを乗算すればデータサイズの概数がわかります。

まずは、PostgreSQLのブロックサイズを見ます。この値はPostgreSQL構築時のパラメータで、テーブルによらず同じ値です。

SHOW block_size;

psqlコマンドですと、次のように表示されます。

=> SHOW block_size;
 block_size 
------------
 8192
(1 row)

1ブロックが8192バイトです。

次に、pg_class という名前のシステムカタログに、テーブルごとのデータサイズなどの統計が入っていますので、それを見ます。

pg_class の reltuples というカラムには、レコード数があります。

pg_class の relpages というカラムには、ディスクを占めるブロック数があります。これにブロックサイズを乗算すればバイト単位のディスク上のデータサイズになります。

いずれも実際の正確な値ではなくPostgreSQLのプランナが使用する推測値です。VACUUMコマンド、ANALYZEコマンドなどで値が更新されるようです。

pg_class の relname はテーブル名です。

次のSQLで全テーブルのサイズを確認することができます。

SELECT relpages * 8192 AS table_size, reltuples, relname FROM pg_class ORDER BY table_size DESC;

1テーブルだけ見たい場合は relname で絞ればよいです。

SELECT relpages * 8192 AS table_size, reltuples, relname FROM pg_class WHERE relname = 'sample1';

 table_size | reltuples | relname 
------------+-----------+---------
       8192 |         0 | sample1
(1 row)

テーブルサイズが大きくて2GB以上の場合はバイト単位ですと32ビットの整数に収まらずに次のようなエラーになります。

ERROR:  integer out of range

1000000で割ってMB単位にするか、もしくはデータ型を64ビット整数であるBIGINTにするとよいです。

SELECT relpages * 0.008192 AS table_size_mb, reltuples, relname FROM pg_class ORDER BY table_size_mb DESC;

SELECT relpages :: BIGINT * 8192 AS table_size, reltuples, relname FROM pg_class ORDER BY table_size DESC;

参考 51.11. pg_class - PostgreSQL 12.0文書

2020/10/26追記

データベースごとのデータ容量を見る方法も書きました。

金星・水星の最大光度のタイミングの違いを解く

この記事は、数ヶ月前にニューラルネットワークをScalaでゼロから書いてたときに勾配降下法の応用としてやっていた内容です。

金星・水星の光度変化を数式で表現し、その最大値のタイミングを調べました。最大値は勾配降下法を使わなくても解けそうですが、せっかく勾配降下法を作ったので試したものです。

金星・水星の最大光度とは

水星と金星は地球よりも内側を周回している惑星です。内惑星といいます。内惑星は地球から見ると月と同様に満ち欠けをします。

ただし、内惑星の満ち欠けが月と大きく違うのは、地球との距離が大きく変わるため見た目の大きさが変わることです。内惑星が満月状のときが地球から最も遠く、見た目も最も小さくなります。新月状のときが地球に最も近く、見た目も最も大きくなります。

満月状のときを「外合」、新月状のときを「内合」といいます。

f:id:suzuki-navi:20200906221643p:plain

月が最も明るくなるのは満月のときですが、内惑星はいつ最も明るくなるのでしょうか。遠いけど満月状のとき（外合）か、細い姿だけど地球から近いとき（内合前後）か。

実際は、2つの内惑星で光度変化の挙動が大きく異なります。水星は遠いにもかかわらず満月状のとき（外合）が最も明るくなり、金星は細いけど地球に近い三日月状のとき（内合の前後）が最も明るくなります。

水星は満月状のとき（外合）が最も明るくなります。外合は太陽と同じ方向で、太陽が眩しいため、最も明るいにも関わらず水星を人が観測することはできません。

金星は三日月状のとき（内合の前後）が最も明るくなります。金星が最大光度のころは、宵の明星または明けの明星として、夕方または日の出前にとても明るく輝いて見えます。最大光度の金星は、太陽、月についで明るく見える天体です。

この記事では、内惑星の光度を適当な式で表現し、その式から勾配降下法で光度が最大になるタイミングを計算しました。

式で表現したならばその微分が0になるところを直接計算してもいいのですが、勾配降下法の例題として考えたので、この記事ではその計算していません。

計算式を考える

位置

まず、位置関係についての計算式です。

太陽を座標の原点として、視点となる地球の位置を \displaystyle x=1, y=0 に固定して、平面上を内惑星が円軌道で周回するとします。

\displaystyle r_{vs} は内惑星ごとの定数です。 \displaystyle \theta_s は内惑星の位置を表すパラメータです。他の変数の値はすべてこの2つから計算できます。 \displaystyle \theta_s=0 が内合、 \displaystyle \theta_s=\pi が外合です。

  • \displaystyle \theta_s : 太陽から見て地球と内惑星がなす角
  • \displaystyle \theta_e : 地球から見て太陽と内惑星がなす角
  • \displaystyle \theta_v : 内惑星から見て太陽と地球がなす角
  • \displaystyle x_v, y_v : 内惑星の位置の座標
  • \displaystyle r_{vs} : 内惑星と太陽の距離（日心距離）（一定値）
  • \displaystyle r_{ve} : 内惑星と地球の距離
 \displaystyle
\begin{align}
x_v &= r_{vs}\cos \theta_s \\
y_v &= r_{vs}\sin \theta_s =  r_{ve} \sin \theta_e \\
{r_{ve}}^2 &= (1-x_v)^2+y_v^2 \\
y_v &= r_{ve} \sin \theta_e \\
\pi &= \theta_s + \theta_e + \theta_v \\
\end{align}

実際には地球も内惑星も楕円軌道を描き、もっとずっと複雑です。単純化のため円軌道を仮定しています。

明るさ

次に、明るさに関する計算式です。明るさを計算する一般的なモデルはわかりません。この記事では次の計算式に従うものと仮定しました。実際には惑星の表面や大気の構成によっても変わりそうです。

  • \displaystyle a : 距離を一定にし、内惑星の見た目の形( \displaystyle \theta_v )から計算される明るさで、満月状のときを基準にした相対値
  • \displaystyle b : 地球からの距離と \displaystyle a から計算される明るさ（距離1で満月状の場合に明るさ1とした相対値）
 \displaystyle
\begin{align}
a &= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
b &= \frac{a}{{r_{ve}}^2} \\
\end{align}

\displaystyle a の式はちょっと複雑な形をしていますが、次の積分を計算して式を作りました。この積分も私が適当にモデルを考えたもので、合っているのかどうかはよくわかりません。

 \displaystyle
\begin{align}
a &= \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta \sin(\theta_v + \theta) d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta ( \sin \theta_v \cos \theta + \cos \theta_v \sin \theta ) d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \sin \theta_v \sin \theta \cos \theta + \cos \theta_v (\sin \theta)^2 d \theta \\
&=  \frac{2}{\pi} \int_0^{\pi-\theta_v} \frac{1}{2} \sin \theta_v \sin 2\theta + \frac{1}{2} \cos \theta_v (1 - \cos 2\theta) d \theta \\
&=  \frac{2}{\pi} \left[ -\frac{1}{4} \sin \theta_v \cos 2\theta + \frac{\theta}{2} \cos \theta_v - \frac{1}{4} \cos \theta_v \sin 2\theta \right]_{\theta=0}^{\theta=\pi-\theta_v} \\
&=  \frac{2}{\pi} \left[ \frac{\theta}{2} \cos \theta_v - \frac{1}{4} \sin (\theta_v + 2\theta) \right]_{\theta=0}^{\theta=\pi-\theta_v} \\
&= \frac{2}{\pi} \left( \frac{1}{2} (\pi - \theta_v) \cos \theta_v - \frac{1}{4} \sin (2\pi - \theta_v) + \frac{1}{4} \sin \theta_v \right) \\
&= \frac{2}{\pi} \left( \frac{1}{2} (\pi - \theta_v) \cos \theta_v + \frac{1}{2} \sin \theta_v \right) \\
&= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
\end{align}

ちなみに \displaystyle a は \displaystyle \theta_v の関数ですが、グラフにすると次のような形をしています。1つ目は横軸 \displaystyle \theta_v の単位をラジアンにしたもの、2つ目は横軸 \displaystyle \theta_v の単位を度にしたものです。

f:id:suzuki-navi:20200906221715p:plain

f:id:suzuki-navi:20200906221727p:plain

ここまでの式をもとに、明るさ \displaystyle b は \displaystyle \theta_s から以下の順番に計算できます。 \displaystyle r_{vs} は定数とします。

 \displaystyle
\begin{align}
x_v &= r_{vs}\cos \theta_s \\
y_v &= r_{vs}\sin \theta_s \\
{r_{ve}}^2 &= (1-x_v)^2+y_v^2 \\
r_{ve} &= \sqrt{r_{ve}^2} \\
\theta_e &= \sin^{-1} \frac{y_v}{r_{ve}} \\
\theta_v &= \pi - (\theta_s + \theta_e) \\
a &= (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\
b &= \frac{a}{{r_{ve}}^2} \\
\end{align}

光度の変化をグラフで見る

\displaystyle r_{vs} に実際の水星と金星の値を入れてみて、内惑星の位置 \displaystyle \theta_s と明るさ \displaystyle b の関係をグラフにしました。1つ目は水星、2つ目は金星です。横軸は \displaystyle \theta_s の単位を度にしたものです。

f:id:suzuki-navi:20200906221749p:plain

f:id:suzuki-navi:20200906221803p:plain

水星と金星とでグラフの形が大きく違います。

勾配降下法を使う

最も明るくなる箇所を勾配降下法で探します。

勾配降下法は、与えられた関数の勾配（微分）を数値的に計算して、関数が最大または最小となる箇所を探す方法です。

微分の値を計算するには、式全体を微分してから数値を当てはめればいいのですが、複雑な式を微分するのは面倒なので、ニューラルネットワークのバックプロパゲーション（誤差逆伝播法）と同様に、またTensorFlowの自動微分と同じように、計算ノードごとに微分を計算することにします。

いくつかの微分式をここに書いておきます。

 \displaystyle
\begin{align}
\frac{d}{dx} \sin x &= \cos x \\
\frac{d}{dx} \cos x &= -\sin x \\
\frac{d}{dx} \sin^{-1} x &= \frac{1}{\sqrt{1-x^2}} \\
\frac{d}{dx} \frac{1}{x} &= -\frac{1}{x^2} \\
\frac{d}{dx} \sqrt x &= \frac{1}{2\sqrt{x}} \\
\end{align}

これで、定数 \displaystyle r_{vs} に対して \displaystyle b が最大となる \displaystyle \theta_s を勾配降下法で計算できます。

内惑星の日心距離と最大光度のタイミングの関係を見る

内惑星の日心距離（太陽と内惑星の距離 \displaystyle r_{vs} ）から、内惑星の光度（ \displaystyle b ）が最大となる \displaystyle \theta_s を勾配降下法により計算します。

横軸に \displaystyle r_{vs} 、縦軸に最大光度の \displaystyle \theta_s を置いたものが以下のグラフの青い線です。ついでにそのときの \displaystyle \theta_e 、 \displaystyle \theta_v の線もひいています。

f:id:suzuki-navi:20200906221822p:plain

青: 太陽から見て地球と内惑星がなす角 \displaystyle \theta_s 橙: 地球から見て太陽と内惑星がなす角 \displaystyle \theta_e 緑: 内惑星から見て太陽と地球がなす角 \displaystyle \theta_v

日心距離0.45付近に大きな崖があります。

日心距離が0.45付近より小さい内惑星では、 \displaystyle \theta_s=180° つまり外合のときに最大光度となることがわかります。日心距離がそれより大きくなると、最大光度となる \displaystyle \theta_s が急激に内合のほうに近くなることがわかります。

日心距離が閾値より小さい内惑星は、地球から見ると太陽のすぐ近くを小さく周回しており、地球との距離の変化が小さくなります。そのため、月と同じように満月状のときが最大光度になります。

逆に日心距離が閾値より大きい内惑星は、地球から見ると太陽の周りを大きく周回しており、地球に近いときは極端に近くなります。地球との距離の変化が激しくなるため、距離のほうが光度の変化に大きく寄与することになり、内合に近いタイミングで最大光度になります。

水星は日心距離0.3871で0.45付近のしきい値より太陽に近いため、外合のときが最大光度となります。金星は日心距離が0.7233ですので、内合に近いときが最大光度となります。

惑星の表面の違いによる光度変化の違いを考察する

今回明るさを計算するモデルとして以下の式を使いました。

 \displaystyle
a = (1 - \frac{\theta_v}{\pi}) \cos \theta_v + \frac{1}{\pi} \sin \theta_v \\

実は最初はもっとシンプルな以下の式を使っていました。

 \displaystyle
a = \frac{1 + \cos \theta_v}{2}

\displaystyle a はグラフにすると以下のような形をしています。

f:id:suzuki-navi:20200906221843p:plain

f:id:suzuki-navi:20200906221854p:plain

もっと複雑な最初に紹介した \displaystyle a のグラフを再掲します。

f:id:suzuki-navi:20200906221749p:plain

f:id:suzuki-navi:20200906221803p:plain

見た目ほとんど同じですが、シンプルな \displaystyle a のほうが、真ん中あたりの値が少し大きいです。

シンプルな \displaystyle a を使って計算したグラフを並べておきます。

水星の光度変化。

f:id:suzuki-navi:20200906221947p:plain

金星の光度変化。

f:id:suzuki-navi:20200906222002p:plain

日心距離と最大光度の位置の関係。

f:id:suzuki-navi:20200906222015p:plain

シンプルな \displaystyle a のほうが閾値が小さいです。

金星はどっちであっても結果に変わりありませんが、水星の光度変化は \displaystyle a をシンプルにすることで大きく変わってしまいました。事実とはだいぶ違ってしまっています。

\displaystyle a の計算式のモデルの説明は省略しますが、シンプルな \displaystyle a は濃い大気に覆われている天体を想定したモデル、複雑な \displaystyle a は大気が少なく地表がクレーターなどごつごつしている天体を想定したモデルです。あくまで素人が勝手に考えたモデルなので、違っている可能性高いのですが、水星の光度変化が水星の表面を推測する根拠になりうることには驚きました。

リンク

過去の勾配降下法に関する私の記事