Pythonでの回帰分析

Question

前提・実現したいこと

Python初心者です．
あるデータを対数正規分布の累積分布関数で回帰分析を行いたいと思っています．
近似線を書いたのですが，範囲が0~∞にならず，負の領域になってしまいます．
大変恐縮ですがご教授お願いします．

該当のソースコード

import numpy as np
from scipy import stats
from scipy.optimize import curve_fit

x = np.array([  0,   4,   9,  14,  19,  24,  29,  34,  39,  44,  49,  54,  59,
        64,  69,  74,  79,  84,  89,  94,  99, 104, 109, 114, 119, 124,
       129, 134, 139, 144, 149, 154, 159, 164])

y = np.array([ 0.        ,  0.17304493,  0.28618968,  0.50083195,  0.55407654,
        0.65058236,  0.73044925,  0.83527454,  0.87687188,  0.92845258,
        0.93510815,  0.95174709,  0.96006656,  0.9750416 ,  0.97670549,
        0.98169717,  0.98169717,  0.9843594 ,  0.98868552,  0.9906822 ,
        0.9906822 ,  0.99234609,  0.99234609,  0.99567388,  0.99567388,
        0.99567388,  0.99567388,  0.99567388,  0.99567388,  0.99567388,
        0.99733777,  1.        ,  1.        ,  1.        ])
def cdf(x, a, b, c):
    return stats.lognorm.cdf(x, a, b, c)
param, cov = curve_fit(cdf, x, y)
y_pred = cdf(x, *param)
# param =[  0.58354202,  -9.49552953,  25.29707224]くらい
# y_pred[0] = 0.046559241691627652となり，0になりません

Answer 1

MSE を見る限り、それなりに近似できているように見えますが、もっとパラメータを真値に近い値にしたいということでしょうか？あくまで近似なので、完全に一致することはないと思います。

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from sklearn.metrics import mean_squared_error
from scipy.optimize import curve_fit

x = np.array([0, 4, 9, 14, 19, 24, 29, 34, 39, 44, 49, 54, 59,
              64, 69, 74, 79, 84, 89, 94, 99, 104, 109, 114, 119, 124,
              129, 134, 139, 144, 149, 154, 159, 164])

y = np.array([0. , 0.17304493, 0.28618968, 0.50083195, 0.55407654,
              0.65058236, 0.73044925, 0.83527454, 0.87687188, 0.92845258,
              0.93510815, 0.95174709, 0.96006656, 0.9750416 , 0.97670549,
              0.98169717, 0.98169717, 0.9843594 , 0.98868552, 0.9906822 ,
              0.9906822 , 0.99234609, 0.99234609, 0.99567388, 0.99567388,
              0.99567388, 0.99567388, 0.99567388, 0.99567388, 0.99567388,
              0.99733777, 1. , 1. , 1.])


def cdf(x, a, b, c):
    return stats.lognorm.cdf(x, a, b, c)
param, cov = curve_fit(cdf, x, y)
y_pred = cdf(x, *param)

# 描画する。
plt.plot(x, y, linestyle='--', marker='o', color='b', ms=2, label='data')
plt.plot(x, y_pred, linestyle='--', marker='o', color='g', ms=2, label='prediction')
plt.legend()
plt.show()

mse = mean_squared_error(y, y_pred)
print(mse)  # 0.00030439892191327277

追記1

対数正規分布の累積分布関数としているのに，xが0の時にyが0とならないのでどうにかしたいのです．

対数正規分布の定義域は 0 < x < ∞ なのに、その累積分布関数で cdf(0) = 0 とならないのはおかしいということですね。
cdf(0) = 0 とならない理由は scipy.stats.lognorm 関数に loc というシフトするパラメータを含んでいるからです。リファレンス を参考にしてください。

なので、このまま fit() すると、loc も推定対象なので、近似したものでは cdf(0) = 0 とはなりません。scipy.stats.lognorm にも fit() という関数があるので、こちらを使うと、loc=0 と固定した状態で残りのパラメータを推定できます。

import seaborn as sns

from scipy import stats
from sklearn.metrics import mean_squared_error
from scipy.optimize import curve_fit

# 対数正規分布に従うサンプルを生成する。
shape, scale = 0.5, 1.
sample = stats.lognorm(s=0.5, loc=0, scale=1.).rvs(size=2000)
sns.distplot(sample, norm_hist=True, kde=False)
plt.show()

# loc は固定して、パラメータを推定する。
shape_pred, loc_pred, scale_pred = stats.lognorm.fit(sample, floc=0)
print('shape={}, loc={}, scale={}'.format(shape, loc, scale))
# shape=0.5094946936562328, loc=0, scale=1.008179589277641

x = np.linspace(0, 10, 100)
y = stats.lognorm.cdf(x, s=shape, loc=0, scale=scale)
y_pred = stats.lognorm.cdf(x, s=shape_pred, loc=loc_pred,  scale=scale_pred)

# 描画する。
plt.plot(x, y, linestyle='--', marker='o', color='b', ms=2, label='true')
plt.plot(x, y_pred, linestyle='--', marker='o', color='g', ms=2, label='prediction')
plt.legend()
plt.show()

mse = mean_squared_error(y, y_pred)
print(mse)  # 7.72357454555545e-06
print(y_pred[0])  # 0.00000000e+00 cdf(0) = 0 となっている。

追記2

curve_fit() でも次のようにすることでパラメータを固定できました。
stats.lognorm.fit() のほうでは、fitting のアルゴリズムが異なるのか、データ数が少ない影響等により精度よく近似できませんでした。

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from sklearn.metrics import mean_squared_error
from scipy.optimize import curve_fit

x = np.array([0, 4, 9, 14, 19, 24, 29, 34, 39, 44, 49, 54, 59,
              64, 69, 74, 79, 84, 89, 94, 99, 104, 109, 114, 119, 124,
              129, 134, 139, 144, 149, 154, 159, 164])

y = np.array([0., 0.17304493, 0.28618968, 0.50083195, 0.55407654,
              0.65058236, 0.73044925, 0.83527454, 0.87687188, 0.92845258,
              0.93510815, 0.95174709, 0.96006656, 0.9750416, 0.97670549,
              0.98169717, 0.98169717, 0.9843594, 0.98868552, 0.9906822,
              0.9906822, 0.99234609, 0.99234609, 0.99567388, 0.99567388,
              0.99567388, 0.99567388, 0.99567388, 0.99567388, 0.99567388,
              0.99733777, 1., 1., 1.])

# scipy.optimize.curve_fit を使うやり方
######################################################
def cdf(x, a, b):
    return stats.lognorm.cdf(x, s=a, loc=0, scale=b)
[s, scale], cov = curve_fit(cdf, x, y)
print('s={}, scale={}'.format(s, scale))  # s=0.8933602211719341, scale=14.750787612138023
# 近似した関数の結果
y_pred = cdf(x, s, scale)

# 描画する。
plt.plot(x, y, linestyle='--', marker='o', color='b', ms=2, label='data')
plt.plot(x, y_pred, linestyle='--', marker='o', color='g', ms=2, label='prediction')
plt.legend()
plt.show()

print(mean_squared_error(y, y_pred))  # 0.000730706490266028
print(y_pred[0])  # 0.0

Comments

[kouheichild]

tiitoi　様
ご回答ありがとうございます．
実は，こちらで算出した式に従ってデータを作成したいのですが
対数正規分布の累積分布関数としているのに，xが0の時にyが0とならないのでどうにかしたいのです．
伝え方が悪くて申し訳ございません．

[tiitoi]

回答を追記しました。この認識であっていますでしょうか。

[kouheichild]

稚拙な説明で申し訳ございません．
そちらの認識であっています．ありがとうございます．
ちなみに，ここまでしていただいたのに大変恐縮なのですが，サンプルを生成せずに，（ヒストグラムのみ既知）である場合に，loc=0として他のパラメタを推定することはできるのでしょうか．

[tiitoi]

質問のように x と y のペアが与えられたときに loc=0 としてパラメータを求めるということでしょうか？
一応、追記2という形で記載しました。

[kouheichild]

tiitoi　様

お時間を割いていただきまして，大変ありがとうございます．
まさに，やりたかったことそのものです．本当にありがとうございます．感無量です．
今後もteratailを利用させていただく予定ですので，どうかお時間があれば，再びご教授いただけたら幸いです．
ありがとうございました．