pythonのpandasとnotebookとstatsmodelとscipyを使ってデータ分析(8)

いよいよ最後ですね。指数近似曲線です。
さて、ここからが問題です。
これを最後に持ってきたのは、問題を先送りに・・・してました。

再度、こちらを参考にさせて頂いて
http://daisukekobayashi.com/blog/python-least-square-method-exponential-approximation/
とりあえずコードは、こんな感じになりました。

%matplotlib inline
import pandas
import numpy
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import scipy.optimize
df = pandas.DataFrame({
        'x': [6.559112404564946264e-01,
                 6.013845740818111185e-01,
                 4.449591514877473397e-01,
                 3.557250387126167923e-01,
                 3.798882550532960423e-01,
                 3.206955701106445344e-01,
                 2.600880460776140990e-01,
                 2.245379618606005157e-01],
         'y': [1.397354195522357567e-01,
                 1.001406990711011247e-01,
                 5.173231204524778720e-02,
                 3.445520251689743879e-02,
                 3.801366557283047953e-02,
                 2.856782588754304408e-02,
                 2.036328213585812327e-02,
                 1.566228252276009869e-02]})
axes = []
fig = plt.figure(figsize=(15,6))
for i, title in enumerate(["ols fit", "Linear exp fit", "log fit", "parabora fit"]):
    axes.append(fig.add_subplot(2, 2, i+1))
    axes[i].sharey = axes[0]
    axes[i].plot(df['x'], df['y'], "ro")
    axes[i].set_title(title)
    axes[i].grid()

str = "Sample Data"
str1 = "Fitted Function:\n"

def r2_string(r):
    return "$R^{2} = %0.5f$" % (r)

def ols_string(a, b, c):
    str2 = "$y = %0.4f x %s %0.4f$\n" % (a, "+" if b >= 0 else "", b)
    return [str, str1 + str2 + r2_string(c)]

def exp_string(a, b, c):
    str2 = "$y = %0.4f e^{ %0.4f x}$\n" % (a, b)
    return [str, str1 + str2 + r2_string(c)]

def log_string(a, b, c):
    str2 = "$y = %0.4f ln(x) %s %0.4f$\n" % (a, "+" if b >= 0 else "", b)
    return [str, str1 + str2 + r2_string(c)]

def pab_string(a, b, c, d):
    str2 = "$y = %0.4f x^{2} %s %0.4f x %s %0.4f$\n" % (a, "+" if b >= 0 else "", b, "+" if c >= 0 else "", c)
    return [str, str1 + str2 + r2_string(d)]

def func_ols(a, b, x):
    return a * x + b

def func_exp(a, b, x):
    return a * numpy.exp(b * x)

def func_log(a, b, x):
    return a * numpy.log(x) + b

def func_pab(a, b, c, x):
    return a * x**2 + b * x + c

def square(x):
    return x ** 2

def fit_exp_linear(parameter, x, y):
    a = parameter[0]
    b = parameter[1]
    residual = numpy.log(y) - (numpy.log(a) + b * x)
    return residual

xx = numpy.arange(0.7, 0.2, -0.01)
parameter0 = [1, 1]

model1 = pandas.ols(y=df['y'], x=df['x'], intercept=True)
model2 = scipy.optimize.leastsq(fit_exp_linear, parameter0, args=(df['x'], df['y']))
model3 = smf.ols(formula='y ~ numpy.log(x)', data=df).fit()
model4 = smf.ols(formula='y ~ square(x) + x', data=df).fit()
#print model1
#print model3.summary()
#print model4.summary()

axes[0].plot(xx, func_ols(model1.beta.x, model1.beta.intercept, xx))
axes[1].plot(xx, func_exp(model2[0][0], model2[0][1], xx))
axes[2].plot(xx, func_log(model3.params['numpy.log(x)'], model3.params['Intercept'], xx))
axes[3].plot(xx, func_pab(model4.params['square(x)'], model4.params['x'], model4.params['Intercept'], xx))

axes[0].legend(ols_string(model1.beta.x, model1.beta.intercept, model1.r2), loc="upper left")
axes[1].legend(exp_string(model2[0][0], model2[0][1], 0), loc="upper left")
axes[2].legend(log_string(model3.params['numpy.log(x)'], model3.params['Intercept'], model3.rsquared), loc="upper left")
axes[3].legend(pab_string(model4.params['square(x)'], model4.params['x'], model4.params['Intercept'], model4.rsquared), loc="upper left")
plt.show()

こんな感じになります。
グラフ的には、エクセルと同じのようです。

今回は、「scipy.optimize.leastsq」を新たに使っています。
こちらの方が、「formula」に自由度が高いと思います。
(statsmodelの使い方が分かっていないだけかも知れませんが・・・)
「formula」の自由度から言うと、
「pandas.ols」<「smf.ols」<「scipy.optimize.leastsq」だと思います。


そして、R2は・・・
「R2-!!!」
と叫びたくなります。
出ていないですね・・・


検証に関しては、「pandas.ols」<「smf.ols」
という感じで、「scipy.optimize.leastsq」は?となります。

ちなみに、こうやると・・・

model2 = smf.ols(formula='numpy.log(y) ~ x', data=df).fit()
print model2.summary()

「R-squared」的には、エクセルと同じ値が出たりします。
でも、「x」や「Intercept」は違うのでダメですね・・・


誰か教えて・・・