平均と分散の物理モデル

以下は、前エントリー「なぜ分散は2乗の和なのか」d:id:rikunora:20190409 から出たオマケの与太話。

なぜ分散は2乗なのか、と問われたとき、「分散を最も小さくする点が平均値だから」というのが1つの答でした。
この答を物理的なモデルに当てはめると、こんな絵が描けます。

f:id:rikunora:20190409113610p:plain

・各データと、とある1点をバネで結んだとき、1点は平均点で止まる。

このモデルで分散は、「バネに溜まったエネルギー」に対応付けられます。
分散が足し算できること(加法性)は、「エネルギー保存則」として理解できます。

ところが分散は、こんな物理モデルに当てはめることもできます。

f:id:rikunora:20190413105347p:plain

・各データに重みを付けてバランスをとったとき、重心は平均点と一致する。

このモデルの場合、分散は「慣性モーメント」に対応付けられます。
分散が足し算できること(加法性)は、「角運動量保存則」として理解できます。

以下、確率・統計入門(小針あき宏)という本からの引用。

f:id:rikunora:20190413105443p:plain

これだけを見ていると何でもないようだが、これは‘平均の持つ意味’についての、一つの思想を物語っている。
・・・平均値をその分布の代表とする理由は、そのまわりの分散が最小ということによって特徴づけられているわけである。
ρ(x)を質量密度と考えると、分散は慣性能率に対応していることは、説明を待たないであろう。
そして自由な剛体の運動は、加えられた力をF,質量をMとするとき、
   F=Mx''
によって記述される、重心の曲線運動と、
慣性能率をI,重心からの力の加えられたところまでの距離をrとすると
   Fr=Iθ''
によって記述される、重心回りの回転運動によって完全に決定されるが、
どうして《重心のまわり》に回転がおこるか、というと、そのまわりの慣性能率が最小だから、に他ならない。
そして重心~平均、慣性能率~分散、の対応を見るとき、平均の持つ意味について、新たな理解が生ずるだろう。

f:id:rikunora:20190413105501p:plain

エネルギーとは、「物理的に、時間変化に対して不変な保存量」のことです。
角運動量とは、「物理的に、角度の変化に対して不変な保存量」のことです。
物理学には“ネーターの定理”と名付けられた定理があります。
標語的に言うと、「対称性あるところに保存則がある」といった主張です。
ネーターの定理によると、
 ・時間並進の対称性から「エネルギー保存則」が、
 ・回転の対称性から「角運動量保存則」が、
 ・空間並進の対称性から「運動量保存則」が、
それぞれ導かれます。

ならば、分散を「運動量保存則」に結びつけるような物理モデルは無いものだろうか・・・
少し考えてみたのですが、うまいモデルは思いつきませんでした。
それもそのはず、エネルギーと角運動量が2次式なのに対し、運動量保存則は1次式なので、土台あてはまらないのです。
※ 運動量保存則は、”分散”ではなく”平均”に結びつけられます。
※ 平均は、確かに平衡移動について対称、一斉にゲタを履かせても不変ですが、あまりにもアタリマエなので意識されません。

話はこれだけなのですが、以上は何も分散だけに限った話ではなく、2次式で表されるような量全般に言えるように思えます。
もともと2次式に当てはまるような量であれば何であれ、
エネルギーや角運動量に関連付いたもっともらしい物理モデルとして表される、ということです。
一般に、理解の仕方や解釈は、個々人の感性やセンス、才能といったものに委ねられると言われがちです。
しかし、この分散の物理モデルを見て思うことは、
「思考の筋道とは全くの自由奔放ではなくて、自ずと規定されている」ということです。
不自然さを感じさせない発想は、文字通り“物の理”に叶っているのです。
そうした“物の理”に叶った発想だけが、分かりやすく、人に伝えやすく、結果的に生き残るのではないか。
そんな風に思えるのです。

確率・統計入門

確率・統計入門

なぜ分散は2乗の和なのか

Q.なぜ分散は、単純な差(偏差の絶対値)ではなく、差の2乗を計算するのか?
A.分散を最も小さくする点が平均値だから。(単純な差を最も小さくする点は中央値となる。)

“分散”というキーワードは統計学の基礎中の基礎であり、どんな教科書にも“平均”の次くらいに載っていることがらです。
しかしながら、いきなり登場する“分散”の意味が分からず、統計学の入り口で挫折する人は少なくありません。

偏差の2乗の平均、つまり、各値と平均との差の2乗の平均を分散といい、
分散の平方根の正の方を標準偏差という。
統計で、ちらばりを表すものとして、標準偏差や分散が多く用いられる。
      -- 高校の教科書(啓林館)より.

教科書にはこのように書かれているのですが、これで分かった気になるでしょうか。
 ・なぜ、差の2乗を計算するのか?
 ・差そのものであってはいけないのか?
 ・なぜ、分散と標準偏差の2種類があるのか?
最後の疑問については「分散が2乗した結果なので、元の単位に戻すために平方根をとる必要がある」のだと説明されます。
  (標準偏差) = √(分散)
だとすると、そもそも値を2種類作らなければならなかった理由は2乗を計算するからであって、
最初から差そのものを採用していれば1種類で済んだはずです。

各データと平均値との差のことを「偏差」と言います。
分散とは、偏差を2乗した値の平均のことです。
一方、偏差の値そのもの(絶対値)の平均は「平均偏差」と呼ばれています。
分散(と標準偏差)に比べて、平均偏差はあまり有名ではありませんし、全ての教科書に載っているわけでもありません。

f:id:rikunora:20190409113422p:plain

もし何の予備知識も無しに「ばらつきを数字で示せ」と言われたなら、
データの差そのものを持ってくる平均偏差の方がよほど自然に思えます。
なのになぜ、今日の統計では、ややこしい分散(と標準偏差)の方が主流なのでしょうか。
歴史を振り返ってみると、偏差の2乗が一般化する以前に、差そのものをばらつきの指標と捉えていた時期がありました。
例えばラプラスは、差そのものを誤差損失の指標として扱っていたという経緯があります。
はっきりと2乗を指標とするようになったのは、ガウスの登場以降のことなのです。

Laplaceはこのことを同じような方法で考察した。しかし彼は、つねに正として扱われる誤差自身を損失の量として選んだ。
    -- ガウス誤差論より.



■ 最小2乗法のこころ

ではなぜ、ガウスは2乗を考えたのか。
その心を知るには、ガウスお家芸である“最小2乗法”をひもとくのが早道です。
最小2乗法とは、誤差を含んだデータに対して、もっともよく当てはまる関係を見出す方法です。
例えば、とあるグループの身長と体重のデータが、こんな風だったとしましょう。

f:id:rikunora:20190409113458p:plain

ここで「身長が高いほど体重が重い」という関係を1本の直線で示すには、どのように線を引くのがもっとも合理的でしょうか。
最小2乗法では、次のように考えます。

f:id:rikunora:20190409113520p:plain

一本の棒と、それぞれのデータの間をバネで結んで引っ張ると、棒はほどよいところに落ち着きます。
この落ち着いた先の状態を「もっとも当てはまりが良い」と見なすのが、最小2乗法の考え方です。
なぜ2乗なのかというと、バネは引っ張れば引っ張るほど、長さに比例した力がかかるからです。
バネの長さを2倍に引き伸ばすには、
 ・長さが2倍
 ・力も2倍
合わせて2×2=4倍のエネルギーが必要です。

f:id:rikunora:20190409113541p:plain

これが2乗の由来です。
最小2乗法とは、つまり「バネに溜まったエネルギーが最小となる位置を探す方法」のことだったのです。
* ∫ニョロっと伸びるはエネルギー、バネは答を知っている
>> http://miku.motion.ne.jp/multi/01_MinSquare.html

いま、2次元のグラフで見た方法を、さらに単純化して1次元としてみましょう。
複数個のデータと、自由に動かせる、とある1点をバネで結んだとしたら、とある1点はどこに落ち着くか。
落ち着く先は、ちょうどデータの平均値となります。

f:id:rikunora:20190409113610p:plain

※ 2次元の場合、当てはめた直線(回帰直線)は必ず平均値を通ります。
平均値とは、「バネに溜まったエネルギーが最小となる点」のことだったのです。
そして「バネに溜まったエネルギー」とは、各データと平均との差の2乗、つまり分散に比例しています。
かくして手元のデータから、もっとも当てはまりの良い点を探そうと考えたとき、分散の最小化というアイデアに導かれるわけです。


■ 分散の加法性

ここで1つクイズ。
 『厚さが平均1cmの月刊誌を1年分積み上げたら、高さは平均で何cmになるか?』
答は当然、12cm。
では次に、
 『その雑誌のばらつきは±1mmだった。1年分積み上げたとき、ばらつきは±何mmになるか?』
ここでの「ばらつきが±1mm」の意味ですが、話を単純化して、雑誌の厚さには 9mm と 11mm の2種類しか無くて、
その2種類が毎月ランダムに発行されるものとしましょう。

答の選択肢:1年分のばらつきは・・・
 (a) ±12mm、 ±1mm × 12ヶ月だから。
 (b) ±0mm、  +と-が打ち消し合って0になるので。
 (c) ±12mm でも ±0mm でもない、中途半端な数。

 

 

 

 

 

    ・・・ちょっと考えてみよう・・・
 

 

 

 

 

正解は (c)。
試しに ±1 の乱数を12個、何回も足し合わせた平均の絶対値をとると 2.7mm くらいの値となります。
(2項係数を用いて詳しく計算すると 2.70703125... といった値になります。)
つまり、ばらつきに関しては単純な足し算ができないということです。

何とか上手い具合に“ばらつきの足し算”ができないものだろうか・・・
そこで編み出されたのが「2乗すれば足し算できる」というアイデアです。
±1mm の乱数が12個あったとして、それぞれをただ足し合わせるのではなく、2乗して足し合わせてみましょう。
  (±1)^2 + (±1)^2 + (±1)^2 + ・・・ ±1)^2 = 12
これなら個々のばらつきがプラスであってもマイナスであっても、合計は確実に+12となるでしょう。
このように、2乗した値によってばらつきを数えたのが、分散という指標だったのです。
 ・雑誌1冊あたりの分散   =  1
 ・12冊積み重ねた全体の分散 = 12
   ※12冊のサンプルがたくさんあったときの分散。つまり毎年の雑誌の山をたくさん並べたときのばらつき。

分散は足し算することができる。
この性質を分散の加法性と言うことがあります。
むしろ、足し算ができるような指標を探したところ、うまくいったのが分散であった、というのが実情でしょう。

では、雑誌1年分の分散が12だったとして、実際に雑誌の山の高さは±何mm程度なのでしょうか。
そこで分散の平方根をとって、元の次元に戻したのが標準偏差です。
今の場合、
  (標準偏差) = √12 = 3.46mm
といった値となります。
標準偏差 3.46mm は、少し上に書いた平均偏差 2.71mm とは異なります。
標準偏差は偏差2乗和の平均、平均偏差は偏差絶対値の平均です ← しつこいようですが、ここんとこよく間違うので。


■ 分散のエネルギー解釈

 ・分散とは「バネに溜まったエネルギー」のこと。
 ・分散とは「足し算できるように」編み出されたばらつきの指標。

この2つが「なぜ分散が2乗なのか」についての主な答になると思います。
さらに穿った見方をするならば、この2つは別々の事柄ではなく、共通の内容を含んでいると解釈することもできるでしょう。

そもそもエネルギーとは何だったかといえば、
 「物理的に、時間変化に対して不変な保存量」
のことでした。
何らかのデータ操作について、例えば上に述べた“バネと棒”のような物理的なモデルを考えたなら、
そのモデル上で足し算できる保存量は、自ずとエネルギーに対応付けられるはずです。
こと分散については、「変位の2乗で溜まるエネルギーのような保存量」というイメージで把握できると思います。
この線で標準偏差に解釈を加えるなら、
 「各データにエネルギーを等しく分配したとき、1個あたりのデータが変位する大きさ」
ということになります。


■ 最小絶対値法ではダメなのか

以上、分散(と標準偏差)が、いかにばらつきの指標として合理的であるかを見てきました。
ならば、ばらつきの指標は分散だけが唯一無二であり、それ以外にはあり得ないのでしょうか。
実はそうでも無いのです。
分散は、一般的には最も優れた指標なのですが、場合によっては分散以外の指標が有用なこともあります。
先に、分散 ~ 最小2乗法に基づいてデータの当てはまりを考えたのですが、
次に、偏差の値そのもの ~ “最小絶対値法”に基づくデータの当てはまりを考えてみましょう。

1次元の場合、もしデータが2つだけだったなら、
もっとも当てはまりの良い点は、2つのデータの間のどこでも良いことになります。

f:id:rikunora:20190409113648p:plain

このように最小絶対値法には、結果が1点に落ち着かないといった不定性があります。
これが最小絶対値法 ~ 平均偏差が流行らない1つの理由です。
データの数をもう少し増やしてみましょう。
もしデータが5つだったなら、最小絶対値法によるもっとも当てはまりの良い点は、中央の3番目の点になります。

f:id:rikunora:20190409113707p:plain

平均偏差を基準に考えたなら、
 ・データが奇数個の場合、順位が中央の点となる。
 ・データが偶数個の場合、順位が中央に最も近い2点の間のどこでもよい。
つまりこれは「中央値」ということです。
※ 中央値の定義としては、データが偶数個の場合、中央に最も近い2点のちょうど中間(平均)となります。

・分散   :偏差の2乗の平均値  : 最小2乗法  : 最小となる点は平均値
・平均偏差 :偏差の絶対値の平均値 : 最小絶対値法 : 最小となる点は中央値

このように対比させてみれば、決して“分散が優れていて、平均偏差が劣っている”わけでは無いのだと分かるでしょう。

さらに興味深いのは、2次元の場合です。
一般によく用いられている最小2乗法ではなく、最小絶対値法による直線の当てはめも、それはそれでアリなのです。

* 最小絶対値法による回帰分析
>> http://www.orsj.or.jp/~archive/pdf/e_mag/Vol.40_02_261.pdf
* 最小二乗と最小絶対値:ロバスト回帰で外れ値分析
>>https://biolab.sakura.ne.jp/robust-regression.html
* Excelで操る!誤差の絶対値和を最小にする
>> http://godfoot.world.coocan.jp/Abs_Min.htm

最初のリンク先に、ラプラスの考えた最小絶対値に基づく方法が載っていました。

この回帰分析規準(1)を最良のものと考え,そのアルゴリズムを最初に提案したのはLaplace[PiereSimmonMarquisdeLaplace,1749-1827]である。

リンク先の説明を参考に、ラプラスの考えた直線当てはめの方法を試してみました。

・原点を通る直線を考える。
・各データを傾きの順番に並べる Y1/X1 >= Y2/X2 >= Y3/X3 ・・・ Yn/Xn
・X1 + X2 + X3 ・・・ + Xr-1 < Xr + Xr+1 + Xr+2 ・・・ + Xn かつ
 X1 + X2 + X3 ・・・ + Xr > Xr+1 + Xr+2 ・・・ + Xn となる r を探す。
 つまり、傾きの大きい方から X を合計していって、合計値が残りの合計をちょうど越える番号 r を探す。
・最も当てはまりのよい傾きとして、Yr/Xr を採用する。

f:id:rikunora:20190409113731p:plain

これは当てはめ結果の1例です。
青線が最小2乗法、赤線がラプラスの考えた最小絶対値法です。
この例のように極端な外れ値が混じっていた場合、最小2乗法は外れ値に引っ張られるため当てはまりが悪くなります。
中央値が平均値に比べて外れ値に強いように、最小絶対値法には外れ値に対して頑健な性質があります。

回帰分析のあてあまり評価の指標として、RMSE, または MAE という値がよく用いられます。
これらは何なのかというと、分散と平均偏差、あるいは平均値と中央値、という並びに位置づけられるものです。

* 平均平方二乗誤差 = RMSE(Root Mean Square Error)
・RMSE の最小化 = 二乗誤差の最小化 = 誤差に正規分布を仮定した場合の最尤推定
* 平均絶対誤差 = MAE(Mean Absolute Error)
・MAE の最小化 = 絶対誤差の最小化 = 誤差にラプラス分布を仮定した場合の最尤推定
 -- 精度評価指標と回帰モデルの評価
 >> https://funatsu-lab.github.io/open-course-ware/basic-theory/accuracy-index/

あるいは最近流行りの機械学習に登場する L1正則化(Lasso), L2正則化(Ridge) といったキーワードも、やはりこれらの概念の延長上に位置づけられます。


ガウス誤差論より

なぜ2乗なのか、その源流を遡ると、ガウスの誤差論にたどり着きます。

f:id:rikunora:20190409113752p:plain

確かに2乗は優れた概念なのだが、ならば2乗を支える根拠・原理は何なのか。2乗以外の概念はあり得ないのか。
私は長らくこの疑問を抱いてきたのですが、元祖ガウスの言葉にようやく答を見たように思います。
今さらですが、(推測)統計学とはデータから帰納的に導かれるものであり、
他の数学のように原理原則から演繹的に導かれるものではありません。
その意味で、2乗が唯一絶対の方法では無かったのです。
ただ、いったん2乗を離れて他のあらゆる方法と比べてみると、改めて2乗こそが最も単純で適していると認めざるを得ません。
これが分散を2乗で数える理由であり、今日に至るまで統計学の根底に流れ続ける基盤だったのです。

誤差論

誤差論


Python による Laplaceの絶対値回帰

# -*- coding: utf-8 -*-
#
# Laplaceの絶対値回帰

import numpy as np
import matplotlib.pyplot as plt

# 直線回帰を行う
def fit_line(x, y):
    mx = np.mean(x)
    my = np.mean(y)
    myx = np.mean(y * x)
    mxx = np.mean(x * x)
    w0 = (myx - my * mx ) / (mxx - mx ** 2)
    w1 = my - w0 * mx
    return ( w0, w1 )

# 適当な直線+乱数データを用意
N_SAMPLE = 10
ALPHA = 0.0  # 切片は0に限定する
BETA = 0.3
EPS = 0.1

X = np.random.rand(N_SAMPLE)
# Y = ALPHA + BETA * X + np.random.normal(0, EPS, N_SAMPLE)     # 正規乱数
# Y = ALPHA + BETA * X + EPS * (np.random.rand(N_SAMPLE) - 0.5)   # 一様乱数
Y = ALPHA + BETA * X + EPS * np.random.standard_cauchy(size=N_SAMPLE) # コーシー分布
    # この場合、Laplace回帰はけっこういい感じだ
# Y = ALPHA + BETA * X + np.random.laplace(loc=0.0, scale=EPS, size=N_SAMPLE) # ラプラス分布
    # Laplace回帰が最尤値になるという。

####
# まずは普通の線形回帰をやってみよう
lin, seg = fit_line( X, Y )
print( "線形回帰の傾き,切片 = {:.4f}, {:04f}".format( lin, seg ) )

####
# 以下、Laplaceの考えた方法で回帰をやってみる
YperX =Y / X   # 個々の点の傾きリスト

Slp_index = {}  # index->傾きリストを作る
for idx in range(N_SAMPLE):
    Slp_index[idx] = YperX[idx]

# 傾きの値で降順にソートする
x_list = []
idx_list = []
for key, val in sorted(Slp_index.items(), key=lambda x: -x[1]):
    # print( "{{ {} : {} }}, X={}".format( key, val, X[key] ) )
    idx_list.append(key)
    x_list.append( X[key] )

# Xの前半の合計値が後半の合計値を越える境界を探す
r_index = 0
for r in range( N_SAMPLE+1 ):
    pre  = x_list[ 0 : r ]              # 前半
    post = x_list[ r : len(x_list)+1 ]  # 後半
    if sum(pre) >= sum(post):
        # print( "Pre:", pre, "len=", len(pre) )
        # print( "Post:", post, "len=", len(post) )
        r_index = r -1  # インデックスなので1引く。
        break

x_esti = idx_list[ r_index ]
# print( "x_ans={}, X={}, YperX={}".format( x_ans, X[x_esti], YperX[x_esti] ) )

lap = YperX[x_esti]
print( "Laplace回帰の傾き = {:.4f}".format( lap ) )

# グラフ出力
xlist = np.arange(0, 1, 0.01) # 0.01単位で細かくプロットし、線として出力させる
ylist_lsq = [ (lin * x + seg) for x in xlist ] # 線形回帰
ylist_lap = [ (lap * x) for x in xlist ] # Laplaceによる回帰

plt.plot(X, Y, 'o') # 元データ
plt.plot(xlist, ylist_lsq, color="blue", label="lsm")       # 線形回帰
plt.plot(xlist, ylist_lap, color="red", label="laplace")    # Laplaceによる回帰

plt.legend(loc='lower right')

plt.show()

サイコロの和と多次元立方体の関係

複数個の一様分布の和は、厳密にはどのような分布に従うのか。

2個のサイコロを同時に振ったとき、最も出やすいのは合計が7のとき。
なぜなら、1+6, 2+5, 3+4, 4+3, 5+2, 6+1 と、6通りの出方があるからです。
それに比べて合計が2となるのは、1+1 の1通り。合計が12となるのは 6+6 の1通りしかありません。
グラフにすれば、このようなテント型になります。

f:id:rikunora:20190211235435p:plain

同じようにして、3個のサイコロの和をグラフにまとめると、こんな形になります。

f:id:rikunora:20190211235517p:plain

* 計算サイト・サイコロの和の確率
>> http://www.calc-site.com/probabilities/dice_total

さらにサイコロの数を10個、20個と増やしていくと、グラフは滑らかな釣り鐘型曲線へと近づきます。
この釣り鐘型のカーブが「正規分布」なのだ、という話は統計学の教科書に譲るとして、
ここで問題にしたいのは正規分布に至るまでの途中の段階です。
たとえば上のサイコロが3個 ~ 一様分布を3個足し合わせた結果は、厳密にはどのようなカーブなのでしょうか。
あるいは一様分布を7個足し合わせた結果は、正確にはどのような分布なのでしょうか。

先に結果を示します。
以下は、Mathematicaというソフトで一様分布を1~7個まで足し合わせた結果を描いたものです。
(※ただしグラフの高さは正規化されていません。)

f:id:rikunora:20190211235547p:plain

f:id:rikunora:20190211235620p:plain

f:id:rikunora:20190211235636p:plain

f:id:rikunora:20190211235651p:plain

f:id:rikunora:20190211235706p:plain

f:id:rikunora:20190211235723p:plain

一定の美しい規則に従って、釣り鐘型曲線が削り出される様子が見て取れたでしょうか。
N個の一様分布の和は、N本の(N-1)次曲線の組み合わせから成り立ちます。
たとえばN=3(サイコロ3個)の場合、次の3本の2次曲線で、3つの区間をカバーしています。

 y = x^2   -- {x = 0~√3/3 の区間}
 y = x^2 - 3 ( (x-1(1/√3))^2 )   -- {x = √3/3~2√3/3 の区間}
 y = x^2 - 3 ( (x-1(1/√3))^2 ) + 3 ( (x-2(1/√3))^2 )   -- {x = 2√3/3~√3 の区間}


では、どうやってこれらの曲線を求めたのか。
実はこれらの曲線は、N次元の立方体を、対角線に直行して切った切断面の面積のグラフなのです。
まず2個のサイコロで考えてみましょう。

f:id:rikunora:20190211235811p:plain

1個目のサイコロの目を横軸に、2個目のサイコロを縦軸に、全ての出目を図示すると6×6の正方形となります。
この正方形の上で、2つの目の合計が一定となる線は、右下がり45度の直線となります。
たとえば、上の図で青で引いた線は合計が4となる直線。
赤で引いた線は合計が7となる直線です。
そして、これらの合計が一定となる直線の長さは、その合計値が出る場合の数=確率に比例しています。
そこで、改めて正方形の黄色い対角線を横軸にして見れば(つまり図を45度傾けて見れば)、
正方形がそのままテント型の分布を形作っていることに気付きます。

次に次元を1つ上げて、3個のサイコロの和を考えてみましょう。

f:id:rikunora:20190211235829p:plain

3個のサイコロを3次元の縦、横、高さに配置すれば、6×6×6の立方体となります。
この立方体の上で、3つの目の合計が一定となる面は、立方体の対角線に直行する面となります。
立方体の対角線を軸にして、切断面の面積をグラフに描けば、3個のサイコロ(一様分布)の和の分布となります。
立方体の切断面は、図の手前から順に3つのパートに分かれます。
1番手前の青色のパートは、1つの頂点からスタートし、三角錐がどんどん大きくなる過程。
2番目の中間にある黄色のパートは、断面の正三角形から、もう一方の反対向きの正三角形に向けて、6角形の断面が徐々に変形する過程。
3番目の赤色のパートは、1番目と反対に三角錐がどんどん小さくなる過程です。
1番目と3番目のパートが2次曲線となることは、断面積が辺の長さの2乗となることから想像が付くでしょう。
ちょっと悩むのは2番目の、中間にある6角形断面のパートです。

f:id:rikunora:20190211235849p:plain f:id:rikunora:20190211235907p:plain

この2枚のスケッチは、6角形断面のパートを取り出して描いたものです。
特に2枚目のスケッチをよく見てください。
断面の6角形とは、拡大する1つの正三角形の3つの頂点から、3つの正三角形を差し引いた残りだと見なせます。
拡大する1つの正三角形と、差し引く3つの正三角形は、実は全く同じ図形を空間的にずらしたものです。
というのも、差し引く3つの正三角形は、元の立方体の隣に位置する、同じ形の立方体の断面だからです。
空間を埋め尽くしている立方体の群れを切断することを想像すれば、
とある1つの立方体の断面積は、同じ形の正三角形を足したり引いたりすることで算出できるのです。
より具体的に、6角形断面の中間パートは次のように定式化できます。

 (1番目のパートの正三角形) - 3×(途中から始まった2番目のパートの正三角形)

立方体の対角線の長さは √3 です。
そして、対角線に対する3つのパートの幅はそれぞれ等しく √3/3 ずつです。
なので (途中から始まった…)とは、数字の上では √3/3 だけずらすことを意味します。
以上に基づいて式を完成させたものが、当初に挙げた「一様分布3個」のグラフだったのです。

※3個のサイコロについては、以下がわかりやすいです。
※『ちょっとした「つり鐘型曲線」を描いてみる』
https://teenaka.at.webry.info/201704/article_12.html
※『ちょっとした「つり鐘型曲線」と正規分布の関係』
https://teenaka.at.webry.info/201704/article_13.html

4個のサイコロ(一様分布)の足し合わせは、4次元立方体を対角線に直行するように切断した断面として表わされます。

f:id:rikunora:20190211235935p:plain

4次元立方体を直接描くことはできませんが、断面である3次元図形を想像することは可能です。
4次元立方体の断面は、以下の4つのパートから成り立ちます。

 第1パート:正4面体がだんだん大きくなる過程。
 第2パート:正4面体の各頂点が削れ、正8面体となる過程。
 第3パート:正8面体の4つの面が成長して、正4面体となる過程。
 第4パート:正4面体がだんだん小さくなる過程。

これらの過程は、3次元の場合と同じように、最初の正4面体の足し引きによって定式化できます。

 第1パート:(正4面体)
 第2パート:(正4面体) - 4×(1つずらした正4面体)
 第3パート:(正4面体) - 4×(1つずらした正4面体) + 6×(2つずらした正4面体)
 第3パート:(正4面体) - 4×(1つずらした正4面体) + 6×(2つずらした正4面体) - 4×(3つずらした正4面体)

式中に登場する正四面体の数、4、6、4 は、二項係数となっています。
二項係数とは、n個の中からk個の要素を取り出す組み合わせの数のことです。
なぜ二項係数となるのか。
4次元立方体の各頂点に、4桁の0または1で番号を付けることを考えます。
4つの頂点は4つの次元に対して、それぞれ手前と奥のどちらか(あるいは上下、右左などの)一方に位置付けられます。
最も手前にある頂点 0000 から出発して、1111に向かって辺で結ばれた頂点を線で結んでたどってゆくと、以下のような図形が描かれます。

f:id:rikunora:20190211235958p:plain

(このような図を“ハッセ図”と言います。)
この図から、同一切断面に並ぶ頂点の数は「4つの桁の中に何個かの1が入る組み合わせの数」になっていることが分かります。
そしてこの頂点の数は、そのまま足し引きすべき図形(正四面体)の数となっているわけです。

以上で見てきた関係は、そのまま次元数が上がっても変わりません。

・最初のパートはN次元の角錐から始まる。グラフの上で、最初の区画は(N-1)次の曲線となる。
・最初のパートの図形(N次元の角錐)を、二項係数に従って足し引することによって、2番目以降のパートをつなげる。
 (足し引きは、プラスマイナスを交互に繰り返せばよい。)

こうして、N次元立方体の断面から、N個の一様分布の和が描き出されます。
※ 最初のパートとなるN次元の角錐の体積は、前の記事に記載しました >> d:id:rikunora:20190210

サイコロの数をうんと増やした極限が正規分布になる、ということは、正規分布とは“無限次元立方体の切断面”の形だということになります。
単なるサイコロの足し合わせを、多次元立方体の切断面として捉えるならば、また違った風景が見えてくる。

多次元正三角錐の体積

2次元の正3角形、3次元の正4面体を延長した、4次元空間にある正三角錐の体積はいくつになるか。
さらに延長して、一般にN次元の正三角錐の体積はいくつになるか。

f:id:rikunora:20190210194659p:plain

4次元空間にある正三角錐のような図形は、正五胞体というのだそうだ >> wikipedia:正五胞体
ウィキペディアには『超体積: (√5/96)a^4』と書かれているが、これはどうやって計算したのだろうか。

いきなり4次元では想像が付かないので、目に見える2次元、3次元から類推しよう。
まず2次元の正方形、3次元の立方体を考える。
とある1つの頂点から伸びる辺を全て含むように平面で切り取った角錐のことを“コーナーブロック”と呼ぶことにしよう。

f:id:rikunora:20190210194737p:plain

“コーナーブロック”とは、上の図で青く塗った部分のことだ。
2次元の場合は直角2等辺三角形、3次元の場合は切り口が正三角形となるような三角錐である。

f:id:rikunora:20190210194802p:plain

4次元の場合、コーナーブロックの切り口は正4面体になる。
そして5次元であれば、コーナーブロックの切り口は4次元空間にある正三角錐(正五胞体)となる。
つまり、コーナーブロックの切り口が、知りたかったN次元の正三角錐の体積となっているわけだ。
そこで、[コーナーブロックの体積] -> [コーナーブロックの切り口] -> [N次元の正三角錐] という方針で、
N次元の正三角錐の体積を計算してみよう。

コーナーブロックが、もとの立方体の何分の1を占めるかを考えてみる。
2次元の正方形の場合には、見ての通り 1/2 だ。
3次元の立方体の場合には、(1/2)×(1/3) = 1/3! となる。
この延長で、4次元のコーナーブロックを考えると (1/2)×(1/3)×(1/4) = 1/4! になる。

なぜ (1/N次元) が掛け合わされるのか。
たとえば3次元空間の角錐は、1つの立方体を3つの方向に、同じ形に切り分けた一部分なので、元の立方体の 1/3 だと分かる。

f:id:rikunora:20190210194826p:plain

3次元のコーナーブロックは、2次元のコーナーブロックの底面に対して (1/3) を掛けるので、
結局のところ元の立方体の (1/2) × (1/3) となる。
4次元空間は絵には描けないが、きっと4つの方向に同じ形に切り分けた部分なのだから、
3次元のコーナーブロックの底面体に対して 1/4 を掛けたものになる。
以下、N次元空間ならコーナーブロックの体積は、元の立方体の 1/N! になる。

次に、コーナーブロックの、立方体の対角線に対する高さを考える。
2次元の(1辺の長さが1の)正方形の対角線の長さは√2、
3次元の(1辺の長さが1の)立方体の対角線の長さは√3、
N次元の(1辺の長さが1の)超立方体の対角線の長さは√N 。
コーナーブロックの高さは、その対角線の長さの 1/N となる。
2次元であれば、(高さ) = √2/2、
3次元であれば、(高さ) = √3/3 となる。

なぜ、コーナーブロックの高さが、ちょうど対角線の 1/N となっているのか。
その理由は、N次元立方体が等間隔に空間を埋め尽くすからである。
3次元で考えてみよう。

f:id:rikunora:20190210194846p:plain

図の立方体[A]に着目して、1つの頂点から反対側の頂点まで、→a, →b, →c と、3本の辺を伝って移動することを考える。
この →a, →b, →c 3本の辺の長さを、立方体[A] の対角線に投影したとき、3本の長さは等しくなる。
たとえば →b を取り上げてみたとき、立方体[A] の隣にある立方体[B] を想像すれば、
→b とは、ちょうど「立方体[B] にとって、立方体[A] の →a 」と同じ位置にある。
なので、→a と →b の、対角線に投影した長さは等しい。
→c についても同じである。
→c とは、ちょうど「立方体[C] にとって、立方体[A] の →a 」と同じ位置にある。
なので結局のところ、元の立方体[A] の対角線は、→a, →b, →c の投影によって3等分される。
この、「立方体が空間を埋め尽くす」という事情は、次元が幾つ上がっても変わらない(だろう)。
なので、N次元のコーナーブロックの高さは、いつでも対角線の 1/N となるわけだ。

以上で、コーナーブロックの体積と高さが分かった。
この体積を高さで微分することによって、切り口の面積(体積?!)を知ることができる。
まず馴染み深い(?!)3次元空間でやってみる。

f:id:rikunora:20190210194902p:plain

グラフの横軸xは対角線方向の長さ(高さ)、縦軸は切り口の面積S、Sの積分がコーナーブロックの体積Vに相当する。
3次元の場合、切り口の面積は高さの2乗に比例するので、S = a x^2 と置く。(今の時点で a は未知の定数。)
あとは高さxと、体積Vに、上の値を入れて a を計算すればよい。
  ∫[0~(√3/3)]{ a x^2 }dx = 1/3!
  a/3 (√3/3)^3 = 1/6
  a = 3√3 / 2
切り口の面積Sは、
  a x^2 = (3√3 / 2) (√3/3)^2 = √3/2

これが切り口の正三角形の面積に一致するはずだ。
一辺が1の正三角形の面積は √3/4 なのだが、今考えているコーナブロックの切り口にある正三角形の1辺の長さは√2である。
なので、
  S = (√3/4) (√2)^2 = √3/2
確かに一致しているので、どうやらこの方法で良さそうだ。

上と同じことを4次元でやってみよう。
この場合、切り口S には正4面体の体積が現れる。
  ∫[0~(√4/4)]{ a x^3 }dx = 1/4!
  a/4 (1/2)^4 = 1/24
  a = 8/3
  S = a x^3 = (8/3) (√4/4)^3 = 1/3

1辺が1の正4面体の体積は √2/12。
今求めた切り口S は、一辺が√2の正四面体なので、(√2/12) (√2)^3 = 1/3
確かに一致している。

さらに、5次元でやってみよう。
これで、知りたかった4次元正三角錐の体積が得られる。
  ∫[0~(√5/5)]{ a x^4 }dx = 1/5!
  a/5 (√5/5)^5 = 1/120
  a = 25√5/24
  S = a x^4 = (25√5/24) (√5/5)^4 = √5/24

これは一辺が√2の4次元正三角錐なので、一辺が1の体積は、
  S = (√5/24) / (√2)^4 = √5/96
確かに、Wikipediaに書かれている値と一致した。

以上を一般化して、(N-1)次元の正三角錐の体積が知りたかったなら、
  ∫[0~(√N/N)]{ a x^(N-1) }dx = 1/N!
  a = (√N)^N / (N-1)!
  S = a (√N/N)^(N-1) = N^N / {(N-1)! N^(N-1) √N}   <-- 一辺が√2 の(N-1)次元の正三角錐
一辺の長さを1にするには、このS を (√2)^(N-1) で割ればよい。
  V = S / (√2)^(N-1) = N^N / {(N-1)! (N √2)^(N-1) √N}

正規分布以外の区間推定

t分布を用いた平均値の区間推定には、「正規母集団を仮定する」というただし書きが付いています。
ということは、母集団が正規分布でなかった場合、区間推定は使えないのでしょうか。
母集団の分布は未知であることが多いので、もし正規母集団という条件が絶対だったなら、
区間推定は実際の現場ではほとんど使えない理屈ということになってしまいます。

結論から言えば、たとえ母集団が正規分布で無くても、
母集団から抜き出した標本の平均値が正規分布に従うなら、区間推定の方法を使うことができます。
中心極限定理によれば、ほとんどの分布は(コーシー分布やべき分布といった“平均の無い”分布を除けば)
標本の平均が正規分布に近づくので、区間推定は正規母集団以外にも広く適用可能なのです。
ただし、正規分布から外れるほど大きなサンプルサイズ(データの数)が必要です。
目安は以下の通りです。
 ・一様分布なら30個以上
 ・ベータ分布(a=4,b=2)くらい歪んだ分布なら100個以上
 ・指数分布のように一方的に偏った分布なら300個以上

これらの目安はどうやって得られたのかというと、ひたすらシミュレーションを繰り返した結果です。
以下に掲げる図は、正規分布以外の特定の分布から、何度も繰り返しサンプルを抜き出して区間推定を行った結果です。
図に示す曲線が元になった母集団の分布。
赤でぼんやり描かれているのは、標本を抜き出す操作を100回繰り返したときの、標本平均のヒストグラムです。
100回の繰り返しを10セット行い、ヒストグラムは10枚重ねて描いてあります。

正規分布


まずこれが正規分布、サンプルサイズ30個の場合。
95%信頼区間であれば、100回のうち95回は真の平均値が信頼区間の中に含まれます。逆に言えば、5回は外れます。
データ30個の場合、外れ回数は5.2回。標本平均の標準偏差=0.18。

●一様分布



一様分布、データ30個の場合と、100個の場合。
データ数が増えるほど平均値の推定幅が小さくなる様子が見てとれます。
データ30個の場合、外れ回数は4.5回。標本平均の標準偏差=0.10。
データ100個の場合、外れ回数は5.2回。標本平均の標準偏差=0.06。
データ数30個であっても区間推定は十分有用でした。

●ベータ分布(4,2)

ベータ分布(4,2)、真の平均値は0.666...、歪んだ分布の例として取り上げました。
データ30個の場合、外れ回数は4.2回。標本平均の標準偏差=0.031。
データ100個の場合、外れ回数は5.5回。標本平均の標準偏差=0.018。
データ300個の場合、外れ回数は4.9回。標本平均の標準偏差=0.010。
30個だと結果が不安定(ばらつきが大きい)でしたが、100個あればかなり安定します。

●指数分布

指数分布、真の平均値は1。
データ30個の場合、外れ回数は8.7回。標本平均の標準偏差=0.18。
データ100個の場合、外れ回数は6.6回。標本平均の標準偏差=0.10。
データ300個の場合、外れ回数は5.3回。標本平均の標準偏差=0.06。
元の分布がここまで偏っていたなら、安定した結果には300個のデータが欲しいところです。

●コーシー分布

コーシー分布には“平均値がありません”。なので、区間推定は使えません。
上はデータ30個の場合ですが、確かに見ての通り、区間推定は適用できません。
* コーシー分布とarctan微分 >> [id:rikunora:20171210]

検索したところ、ここに詳しく書かれていました。
* 大人の知らない統計学 6  非正規分布の信頼区間
>> https://blogs.yahoo.co.jp/pironotakarabako/63655916.html

以下にシミュレーションのソースを示します。
(分布や数値の付け替えはコメントで行っています。わかりずらい。)
シミュレーションといっても、要は区間推定をひたすら繰り返しているだけです。

# -*- coding: utf-8 -*-
"""
区間推定の頑健性、正規母集団以外では、どれほどサンプルがあれば妥当なのか。
"""
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

class IntervalCheck:

    N_TRIAL_SET =  10   # 試行セット回数
    N_TRIAL     = 100   # 試行回数
    N_SAMPLE    = 30    # 100   # 300   # サンプルサイズ=データ数
    
    up_miss   = 0   # 上側に失敗した回数
    down_miss = 0   # 下側に失敗した回数
    
    x_min = -2.0    # -1.5    # 0.0     # グラフの範囲
    x_max = +2.0    # +1.5    # 1.0
    y_scale = 20    # 10    # 分布曲線をy方向に拡大して描く
    
    # N_TRIAL_SET だけ試行セットを繰り返す
    def run(self):
        misses = np.zeros( self.N_TRIAL_SET )
        
        for i in range( self.N_TRIAL_SET ):
            misses[i] = me.each_run()
        
        print( u"信頼区間からの外れ, mean={:.03f}, stdev={:.03f}".format( np.mean(misses), np.std(misses) ) )
        
        plt.xlim( self.x_min, self.x_max )
        self.draw_dist()    # 分布曲線を描き足す
        plt.show()
        
        return
    
    # 分布曲線を描く
    def draw_dist(self):
        n = np.linspace( self.x_min, self.x_max, 1000 )
        p = []
        for i in range(len(n)):
            x = stats.norm.pdf( x=n[i], loc=0, scale=1 )    # 正規分布
            # x = stats.uniform.pdf( n[i], loc=-1, scale=2 )  # 一様分布
            # x = stats.beta.pdf( n[i], 4, 2 )    # ベータ分布
            # x = stats.expon.pdf( n[i], loc=0, scale=1 )  # 指数分布
            # x = stats.cauchy.pdf( n[i] )    # コーシー分布
            p.append( self.y_scale * x )
        
        plt.plot( n, p )
        return
    
    # 1回の試行セット
    def each_run(self):
        # print( 'sample_size={}'.format(self.N_SAMPLE) )
        self.up_miss = self.down_miss = 0
        
        means = np.zeros( self.N_TRIAL )
        
        for i in range(self.N_TRIAL):
            means[i] = self.trial()
        
        miss_ratio = (self.up_miss + self.down_miss) / self.N_TRIAL
        print( 'miss_ratio={:.02f}, up_miss={}, down_miss={}, '.format( miss_ratio, self.up_miss, self.down_miss ), end="")
        
        plt.hist( means, color='r', alpha=0.1, histtype='stepfilled' )  # 半透明で重ね書き
        print( 'mean_means={:.05f}, std_means={:.05f}'.format( np.mean(means), np.std(means) ) )
        
        return miss_ratio
    
    # 1回の試行
    def trial(self):
        sample = np.random.normal( 0, 1, self.N_SAMPLE ) # N個の正規分布 (平均, 分散, 出力数)
        # sample = np.random.rand( self.N_SAMPLE ) * 2 - 1 # N個の一様乱数
        # sample = np.random.beta( 4, 2, self.N_SAMPLE )   # β分布、非対称でやってみよう
        # sample = np.random.exponential( scale=1.0, size=self.N_SAMPLE )    # 指数分布
        # sample = np.random.standard_cauchy( size=self.N_SAMPLE )    # コーシー分布は平均が無い
        
        real_mean = 0   # 理論上の平均
        # real_mean = 4/(4+2) # a/(a+b) ベータ分布の平均はこうなる
        # real_mean = 1   # 指数分布の平均はscale
        
        # t分布による区間推定
        mean = np.mean(sample)  # 標本平均
        t_dist = stats.t( loc = mean,
                 scale = np.sqrt( np.var(sample, ddof=1) / self.N_SAMPLE ), # √(不偏分散/サンプルサイズ)
                 df= self.N_SAMPLE - 1 )    # 自由度
        down, up = t_dist.interval(alpha=0.95)
        # print( '95% interval: {:.3f} < x < {:.3f}'.format(down, up) )
        
        # 区間推定の結果が当たっているかどうか
        if up < real_mean :
            self.up_miss += 1
        elif down > real_mean :
            self.down_miss += 1
        
        return mean # 標本平均を返す

if __name__ == '__main__':
    me = IntervalCheck()
    me.run()


今年はブログを書くヒマがほとんど無かった。来年はもっと時間のゆとりが欲しい。
それでは、よいお年を!

エクセル近似曲線の罠

エクセル近似曲線の「指数近似」「累乗近似」は、いわゆる非線形最小二乗法ではない。
 ・エクセルで用いているのは、データの対数に直線をあてはめるという方法。
 ・いわゆる非線形最小二乗法とは、残差(誤差)の2乗和を最小にする方法。

詳しいことは以下のページで尽きているのだが、エクセルの近似曲線は便利だと、
私自身多くの人に勧めている手前、注意を忘れないよう記載しておく。

* Excelグラフ累乗,指数,多項式近似の論文記載の注意(生物科学研究所 井口研究室)
>> https://biolab.sakura.ne.jp/excel-graph.html

* 指数近似

赤い点で描かれているのが元になるデータ、
 Y = 0.8 * exp( 0.3 * X ) に、Xに比例する大きさで正規乱数を加えたもの。
下側の青い線がエクセルの指数近似曲線。
上側の緑の線はR言語による非線形回帰。

* 累乗近似

赤い点で描かれているのが元になるデータ、
 Y = 0.3 * X ^ 1.5 に正規乱数を加えたもの。
直線に近い青い線がエクセルの累乗近似曲線。
よりカーブのきつい、緑の線はR言語による非線形回帰。

ぱっと見に、あてはまりが良さそうなのはR言語の方に思える。
だからといってエクセルが“間違っている”というわけでもない。
エクセルはR言語と異なる方法で線を引いているのである。
R言語がもともとの回帰式の上で誤差を最小化するのに対し、
エクセルは両辺対数をとった式の上で誤差を最小化する。

* 指数近似
  Y = a * X ^ b   -- もともとの回帰式
  ln(Y) = ln(a) + b ln(X)   -- 両辺対数をとった式
* 累乗近似
  Y = a * exp(b * X)   -- もともとの回帰式
  ln(Y) = ln(a) + b * X   -- 両辺対数をとった式

もともとの回帰式は曲線だが、両辺対数をとった式は直線となる。
誤差の最小化は、直線の方がずっと簡単だ。
直線に直した上での回帰は以下のようになる。

* 指数近似


* 累乗近似


この直線で求めたパラメータが、エクセル近似曲線のものに一致していることが見て取れる。
 ※ exp(0.2526)=1.2874, exp(-0.1613)=0.851
エクセル近似曲線は、直線回帰にちょっとオマケを付け足して実現した機能だったのだ。
エクセルの指数近似、累乗近似では、ゼロやマイナスを含むデータに線を引くことができない。
なぜできないかというと、ゼロやマイナスの対数がとれないからである。
また、エクセルには指数回帰を行う"LOGEST"という関数があるが、この関数で求めた数字も、近似曲線で求めた数字と同じである。
(もう1つの、直線回帰のための"LINEST"関数は特に問題無い。)



さらに、エクセルでは近似曲線に対して決定係数R2を計算することができるのだが、
決定係数R2は本来、線形回帰についてのあてはまり指標であり、非線形回帰に用いることはできない。
では、エクセルに出てくる決定係数は何なのかというと、対数をとって直線に直したときの決定係数だったのだ。
なので、この数字をうのみにして、そのまま用いるべきではない。
非線形回帰であてはまりの良さを評価するには、決定係数ではなく、
残差の2乗和、あるいは(回帰分析の)標準誤差を用いるのが妥当だ。

以上の懸念は、指数近似、累乗近似の話であって、他の種類の曲線「対数近似」「多項式近似」では心配しなくてよい。
実際試してみると、対数近似、多項式近似の結果はR言語と一致する。
なぜそうなるかというと、

 ・対数近似,多項式近似 は 線形回帰
 ・指数近似,累乗近似 は 非線形回帰

だからである。
形の上で、線形=直線、非線形=曲線、と思い込んでいる人もいるかもしれないが、
回帰分析の場合、回帰式が推定パラメータについて1次のものが線形、そうでないものが非線形である。
 ・対数近似の回帰式: y = a ln(X) + b
 ・多項式近似の回帰式: y = a X^2 + b X + c
いずれの回帰式も(Xではなく)a, b, c については1次式となっている。

で、指数近似,累乗近似がやりたかったとき、結局のところどうすれば良いのか。
・エクセルを使ったときは、エクセルの結果だと明記しておく。
・できればR言語、あるいは信頼できる統計用ソフトを使った方が良い。
・実はエクセルの「ソルバー」機能を使えば、非線形回帰も不可能ではない。
 でも、そこまでするなら素直にR言語使えば良いのではないかな。

* Excelのチャートオプション"近似曲線の追加"機能の評価
>> https://ci.nii.ac.jp/naid/110007025777

以下はR言語上での操作です。

### 累乗近似、非線形回帰
# データを用意
x <- c(1:12)
y <- c(
	1.075513124,
	2.104125184,
	2.045615008,
	3.70051273,
	2.221801362,
	2.769819828,
	6.656385201,
	7.037107115,
	7.976938733,
	8.150842204,
	11.70331758,
	13.66264136 )

nonline <- nls( y ~ a * x^b, ,start=c(a=1,b=1))  # 非線形回帰
summary(nonline)  # 結果の要約

Formula: y ~ a * x^b

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   0.2680     0.1204   2.226   0.0502 .  
b   1.5587     0.1959   7.957 1.23e-05 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.132 on 10 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 2.96e-06

# プロットしてみる
a <- coef(nonline)[1]
b <- coef(nonline)[2]
plot( y ~ x)
curve( a*x^b, add=T, col="red")

### 指数近似、非線形回帰
# データを用意
x <- c(1:12)
y <- c(
	3.283583359,
	2.754235941,
	1.674398192,
	3.556187361,
	0.976339451,
	4.652897765,
	15.21122786,
	12.09124224,
	8.109277666,
	13.59872153,
	23.046961,
	24.17609857 )

nonline <- nls( y ~ a * exp(b*x), ,start=c(a=1,b=1))  # 非線形回帰
summary(nonline)  # 結果の要約

Formula: y ~ a * exp(b * x)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  1.53717    0.67319   2.283 0.045521 *  
b  0.23269    0.04152   5.605 0.000226 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.404 on 10 degrees of freedom

Number of iterations to convergence: 15 
Achieved convergence tolerance: 3.005e-07

# プロットしてみる
a <- coef(nonline)[1]
b <- coef(nonline)[2]
plot( y ~ x)
curve(a * exp(b*x), add=T, col="red")

10個の数字の計算パズル

「日本お笑い数学協会」というTwitterに流れていた問題 >> https://twitter.com/owaraisugaku?lang=ja



気になったのでメモを書き残して仕事に出かけて、帰ってきたら家族がものすごく計算してた!


>> https://ameblo.jp/rikunora/entry-12420444729.html

で、私と言えばオヤジのメンツを保つため、こっそりパソコンの力で解いたのはここだけの秘密です。
※ 以前紹介した Coprisで解けます(下を参照)>> [id:rikunora:20181025]
※(答えるときはリプライは避け、引用RTなどでお願いします!)とあったので、答は載せません。

同じTwitterに、こんな問題もありました。

やはりパソコンの力を借りて答を列挙したところ、
合計 17 となる答が 2 通り、
合計 19 となる答が 4 通り、
合計 20 となる答が 6 通り、
合計 21 となる答が 4 通り、
合計 23 となる答が 2 通り、
全部で18通りの答がありました。

いきなり上の三角形に挑むのは少々しんどいので、もう1段簡単にした問題で考えてみました。

赤、青、緑の線の足し算が同じになるように、○の中に1〜6の数を1つづつ入れよう!

答は以下の4通り。

よく見ると並び方に規則があります。

これをもとにすれば、先の1段階大きな三角形でも何とかなると思う。

せっかくなので、私も1問考えてみました。

赤、青、緑の線の足し算と、赤、青、黄の菱形の中の足し算が全て同じになるように、○の中に0〜9の数を1つづつ入れよう!

答は2通りある。
ヒント:合計は 19 になる。


/**
* □□÷□=□□÷□=□×□=□□を解く Scala + Copris スクリプト.
*/
import jp.kobe_u.copris._
import jp.kobe_u.copris.dsl._

// 10個の変数を用意する
for( n <- 0 until 10 ){
	int('x(n), 0, 9)	// 0〜9の整数
}

// 10個は全て異なる
var lst = List[Term]()
for( n <- 0 until 10 ){
	lst :+= 'x(n)
}
add( Alldifferent( lst ) )

// 問題の式
var last_term = 'x(8) * 10 + 'x(9)
add( ('x(0) * 10 + 'x(1)) === 'x(2) * last_term )	// 割り算は直接書けないので
add( ('x(3) * 10 + 'x(4)) === 'x(5) * last_term )
add( 'x(6) * 'x(7) === last_term )

/* 10の位は0にならない。0で割り算しない。この条件は無くても可 **
'x(0) !== 0
'x(2) !== 0
'x(3) !== 0
'x(5) !== 0
'x(8) !== 0
**/

find
printf( "%d%d ÷ %d = %d%d ÷ %d = %d × %d = %d%d\n",
	solution('x(0)), solution('x(1)), solution('x(2)), solution('x(3)), solution('x(4)),
	solution('x(5)), solution('x(6)), solution('x(7)), solution('x(8)), solution('x(9))
)