pandasで移動平均や高値線、安値線を計算する

前回がローソク足だったので今回も市場データでよく使われるテクニックから移動平均を取り上げてみたいと思います。
ついでにHLバンド(ドンチャンチャンネル/高値線,安値線)も同様にもとまるので紹介します。

技術としては、window関数と呼ばれる種類の関数を使って算出します。

ドキュメントはこの辺り。
Window
pandas.DataFrame.rolling
pandas.Series.rolling

DataFrameとSeries両方に実装されていて、同じように使うことができます。
rolling() で 指定期間ごとに区切ったデータを作り、それに対して、 meanやmax,minなどの関数を適用して
平均や最大値、最小値を算出して配列として返します。

実際に見た方が早いと思うので、乱数でランダムウォークデータを生成し、
計算して可視化してみましょう。


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データ作成
data = pd.DataFrame(np.random.normal(0, 100, 200).cumsum() + 20000)
# 移動平均、期間の高値/安値線の計算
ma = data.rolling(10).mean()
high_band = data.rolling(20).max()
low_band = data.rolling(20).min()
# 可視化
fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(1, 1, 1)
ax.plot(data, label="元データ")
ax.plot(ma, label="移動平均")
ax.plot(high_band, label="高値線")
ax.plot(low_band, label="安値線")
plt.legend()
plt.show()

結果がこちら。

期間の最初の方はデータ不足により線が途切れています。
この辺りの制御は min_periods などの引数で細かく調整できるので、
データ量やその時の目的に応じて調整して使いましょう。

pythonでローソク足を描く

以前(pythonを勉強し始めた頃)は、matplotlibでローソク足をかけたはずなのですが、最近は方法が変わってしまったようなのでそのメモです。
なお、ここでサンプルに使うデータはすでにcsvフィアルか何かで保存されているものとします。

以前は matplotlib.finance というのをimport でき、これを使ってかけたのですが、version 2.0 からなくなってしまったようです。
matplotlib.finance
This module is deprecated in 2.0 and has been moved to a module called mpl_finance.

そしてさらに良くないことに、移動先の mpl_finance ですが、あまりしっかり保守されてない様子。

githubのリポジトリに下記の文言があります。
The code is provided as is and is basically un-maintained.

ただ、一応動くようなので動かしてみましょう。
anacondaには含まれていないようなので、インストールから必要です。
pip install mpl_finance
これで、
mpl-finance==0.10.0
が入りました。

さて、使い方ですがun-maintainedの宣言通り、 mpl-finance の公式ドキュメントらしきものは見当たらず、
上の、matplotlib.finance時代のドキュメントを読んで使わないといけないようです。

ローソク足を書く関数は次の4つあり、それぞれデータの渡し方が違います。
.candlestick2_ochl(ax, opens, closes, highs, lows, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick2_ohlc(ax, opens, highs, lows, closes, width=4, colorup=’k’, colordown=’r’, alpha=0.75)
.candlestick_ochl(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)
.candlestick_ohlc(ax, quotes, width=0.2, colorup=’k’, colordown=’r’, alpha=1.0)

今回は手元のデータと相性が良いので .candlestick_ohlc を使います。
quotes に 日付、始値、高値、安値、終値、の5列のデータがデータ件数行だけ並んだ配列を渡してあげる必要があります。
ここで面倒なのは日付の渡し方で、float型で渡す必要があります。
ドキュメントに time must be in float days format – see date2numとある通り、専用の関数があるのでそれを使います。

matplotlib.dates.date2num(d)

また、この関数は引数がdatetime型なので、元々が2019-05-07 のような文字列になっているならば、
datetime型に変換しておく必要があります。
それにはpandasの to_datetimeを使います。
pandas.to_datetime
(いつもならそれぞれ1記事使ってるようなテクニックですね。to_datetimeの方は便利なのでそのうち専用記事書くかも。)

前置きが長くなりましたが、ここまでの情報でできるので日経平均のcsvデータからローソク足を書いてみましょう。


import pandas as pd
import mpl_finance
import matplotlib.pyplot as plt
from matplotlib.dates import date2num

# データの読み込み
df = pd.read_csv("./日経平均データ.csv")
print(df.head())

'''
        date      open      high       low     close
0  2019-3-11  21062.75  21145.94  20938.00  21125.09
1  2019-3-12  21361.61  21568.48  21348.81  21503.69
2  2019-3-13  21425.77  21474.17  21198.99  21290.24
3  2019-3-14  21474.58  21522.75  21287.02  21287.02
4  2019-3-15  21376.73  21521.68  21374.85  21450.85
'''

# dateの型変換
# まずdatetime型にする
print(df["date"].dtypes)  # object
df["date"] = pd.to_datetime(df["date"])
print(df["date"].dtypes)  # datetime64[ns]
# 続いて float型へ
df["date"] = matplotlib.dates.date2num(df["date"])
print(df["date"].dtypes)  # float64

print(df.head())
'''
       date      open      high       low     close
0  737129.0  21062.75  21145.94  20938.00  21125.09
1  737130.0  21361.61  21568.48  21348.81  21503.69
2  737131.0  21425.77  21474.17  21198.99  21290.24
3  737132.0  21474.58  21522.75  21287.02  21287.02
4  737133.0  21376.73  21521.68  21374.85  21450.85
'''

# 可視化
fig = plt.figure(figsize=(13, 7))
ax = fig.add_subplot(1, 1, 1)
mpl_finance.candlestick_ohlc(ax, df.values)
plt.show()

こうして出来上がるチャートが次です。

正直これ単体では手間の割に可視化するメリットがないなーというのが正直なところです。
ただ、matplotlibの仕組みに乗っかっているので、
自分のオリジナルの指標などを追加していくことができます。

matplotlibのpcolorとpcolormesh

先日のトピックモデルの記事中で、試しにヒートマップでの可視化を試みた時、
matplotlibのpcolorって関数を使用しました。

参考:pythonでトピックモデル(LDA)

matplotlibでヒートマップを描こうと思うと少々無理やりな実装になるものも含めて、
imshowや、contourf、pcolorなど複数の方法が考えられ、結構迷いますが一番自然に書ける気がして採用しました。
しかしどうやらこのpcolor、あまり評判がよろしくないようです。

公式ドキュメントを見ても次ように記載があります。

matplotlib.axes.Axes.pcolor

Hint

pcolor() can be very slow for large arrays. In most cases you should use the similar but much faster pcolormesh instead. See there for a discussion of the differences.

要するに pcolormesh を使う方が良いようです。

ドキュメントはこちら。
matplotlib.axes.Axes.pcolormesh

Differences between pcolor() and pcolormesh()
の部分を読んでも、あまり pcolorにメリットを感じないので、
いっそのこと pcolor 自体を pcolormesh のエイリアスか何かに変えてしまっても良さそうなのですが、
戻り値の型が違うこともありますし、何か同じように使えない事情もあるようですね。

とりあえず今後はこれまでpcolorを使っていた場面ではpcolormeshを使うようにしようと思います。
ちなみに、x軸y軸が等間隔の場合はimshowの方がさらに速いそうです。

imshow
If X and Y are each equidistant, imshow can be a faster alternative.

標準化レーベンシュタイン距離

以前の記事で、レーベンシュタイン距離を計算できるライブラリを紹介しました。
参考:pythonで編集距離(レーベンシュタイン距離)を求める

どちらかというと、ライブラリよりもアルゴリズム側の説明の続きなのですが、
標準化されたレーベンシュタイン距離(normalized Levenshtein distance)というものも提案されています。
これは、二つの文字列のレーベンシュタイン距離を、文字数が多い方の文字数で割った値として定義されます。
固有名詞の名寄せなどでレーベンシュタイン距離を使う場合、
こちらを使った方がうまく行くことが多いようです(個人的な経験から。)

以前紹介した、 python-Levenshtein
実装されているんじゃないかと期待してしばらく調べていたのですが、どうやらこれは実装されてないようです。
特にLevenshtein.ratio という関数に期待したのですがこれは全然違いました。

ということで自分で実装しましょう。


import Levenshtein


def normalized_distance(text0, text1):
    dist = Levenshtein.distance(text0, text1)
    max_len = max(len(text0), len(text1))
    return dist / max_len

関数名は normalized_levenshtein_distance にしたかったが流石に長すぎるので少し短縮。
ただ、これでも長いですね。(pep8的につらい)

これによって、例えば通常のレーベンシュタイン距離では、
“バニラ” と “アイス” の組み合わせも “チョコレート” と “チョコレートアイス” もどちらも距離は3でしたが、
標準化した距離を使うことで、
前者は距離1、後者は 1/3(=0.333…)と算出されるようになり、前者の方が離れてると見なせます。

使ってみた結果がこちら。


print(Levenshtein.distance("アイス", "バニラ"))  # 3
print(Levenshtein.distance("チョコレートアイス", "チョコレート"))  # 3
print(normalized_distance("アイス", "バニラ"))  # 1.0
print(normalized_distance("チョコレートアイス", "チョコレート"))  # 0.3333333333333333

scipyで距離行列を計算する

前の記事でちらっと pdist関数が登場したので、scipyで距離行列を求める方法を紹介しておこうと思います。

距離行列の説明はwikipediaにあります。
距離行列 – Wikipedia

要するに、N個のデータに対して、(i, j)成分がi番目の要素とj番目の要素の距離になっているN*N正方行列のことです。

これは次の二つの関数を使って計算できます。
scipy.spatial.distance.pdist
scipy.spatial.distance.squareform

numpyで適当に5点取ってやってみましょう。


import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# 出力する桁数を抑える
np.set_printoptions(precision=3)
# 乱数生成
X = np.random.randint(-5, 6, size=(5, 2))
print(X)
'''
[[-2 -4]
 [-3 -4]
 [ 2 -1]
 [ 4 -2]
 [-1 -2]]
'''

y = pdist(X)
print(y)
'''
[1.    5.    6.325 2.236 5.831 7.28  2.828 2.236 3.162 5.   ]
'''

M = squareform(y)
print(M)
'''
[[0.    1.    5.    6.325 2.236]
 [1.    0.    5.831 7.28  2.828]
 [5.    5.831 0.    2.236 3.162]
 [6.325 7.28  2.236 0.    5.   ]
 [2.236 2.828 3.162 5.    0.   ]]
'''

変数Mに距離行列が入っているのが確認できますね。
(1,2)成分が [-2 -4] と [-3 -4] の距離の1,
(1,3)成分が [-2 -4] と [ 2 -1] の距離の5など、きちんと距離が入っています。

なお、 pdistの戻り値自体は、正方行列の形をしておらず、squareformで正方行列に変形する必要があるので注意です。
pdist の戻り値(コード中では変数y)には、上三角行列部分の値が、1行めから順番に入っています。
配列の長さは$(N*(N-1)/2) = 5*4/2 = 10$ です。

pdist 関数は metric という引数に様々な距離関数の値をとることができ、いろんな種類の距離行列を作ることができます。
(デフォルトはユークリッド距離)

ここまで紹介した、pdistは1つのデータセットに対して距離行列を生成しましたが、
cdistという関数を使うと、2個のデータセットからそれぞれ1個ずつ要素を抽出して距離行列を作ることもできます。

こちらはsquareformを使わなくても初めから行列の形で結果が得られます。


from scipy.spatial.distance import cdist
XA = X[:3]
XB = X[3:]
print(cdist(XA, XB))

# 出力
'''
[[6.325 2.236]
 [7.28  2.828]
 [2.236 3.162]]
'''

最初の例においても、
pdistとsquareform を使うよりも、 cdist(X, X) した方が便利かもしれません。
結果は同じになります。

pandasでカテゴリー変数を数値に変換する

pandasの関数を使ってカテゴリー変数を数値に変換する方法のメモです。
one-hot表現とは異なり、[“a”, “b”, “c”] みたいなデータを、 “a”は0,”b”は1,”c”は2みたいな感じで数字に変換します。

scikit-learnのデータセットにおける、正解ラベルのようなデータ型にする、と言ったほうがわかりやすいかもしれません。
Rにはカテゴリカル変数の型があるのでそちらの方がイメージしやすい人もいるかも。

使う関数は、pandas.factorizeです。

この関数に変換したいリストを渡すと、数値に変換したリストと、何番めの数がもともと何の値だったのかを示す配列を返してくれます。

実際にやってみましょう。
numpyでランダムに10個のデータを作って、それを数値化します。


import numpy as np
import pandas as pd

data = np.random.choice(["a", "b", "c"], 10)
print(data)
# ['b' 'c' 'a' 'c' 'c' 'c' 'c' 'a' 'b' 'a']
labels, uniques = pd.factorize(data)
print(labels)
# [0 1 2 1 1 1 1 2 0 2]
print(uniques)
# ['b' 'c' 'a']

うまくいきました。
uniques の値をアルファベット順にしたい時は、
factorize する時に sort=Trueを指定するとできます。


labels, uniques = pd.factorize(data, sort=True)
print(labels)
# [1 2 0 2 2 2 2 0 1 0]
print(uniques)
# ['a' 'b' 'c']

できました。

numpyで配列内の値を特定の値に制限する

前の記事の最後で、
異常時の前処理として、1〜99パーセンタイルでクリップするって話を少し書いたので、
それをnumpyで実現する関数の紹介です。
と言っても、わざわざ専用関数を使わなくても容易に実装できるのですが、せっかく用意されているのがあるので知っておくと便利です。

ドキュメントはこちら。
numpy.clip

np.clip(配列, 最小値, 最大値)と指定すると、
配列の値のうち、区間[最小値, 最大値]からはみ出た値を、その範囲に収まるように区切ってくれます。

ためしに、標準正規分布に従う値20個を生成して、[-1, 1]の範囲にクリッピングしてみましょう。


import numpy as np
data = np.random.randn(20)
print(data)
'''
[-1.71434343  0.33093523 -0.0648882   1.34432289 -0.15426638 -1.05988754
 -0.41423379 -0.8896041   0.12403786  1.40810052  0.61199047  1.60193951
 -0.72897283 -0.00861939 -0.38774556  0.40188148  0.08256356  1.61743754
 -0.12320721  1.45184382]
'''
data = np.clip(data, -1, 1)
print(data)
'''
[-1.          0.33093523 -0.0648882   1.         -0.15426638 -1.
 -0.41423379 -0.8896041   0.12403786  1.          0.61199047  1.
 -0.72897283 -0.00861939 -0.38774556  0.40188148  0.08256356  1.
 -0.12320721  1.        ]
'''

-1 より小さかった値は-1に、 1より大きかった値は1になりました。

ちなみに下記のコードでも同じことができます。


data[data < -1] = -1
data[data > 1] = 1

前回紹介した、percentileと組み合わせて使うことで、
nパーセンタイルからmパーセンタイルにクリップするということも簡単に実現できます。

試しに 5〜95パーセンタイルにクリップしたのが次のコードです。


data = np.random.randn(10)
print(data)
'''
[-0.41127091 -1.34043164  0.09598778 -1.19662011 -0.04607188 -0.02745831
  0.23184919  0.85601106  0.58430572  0.88205005]
'''
c_min, c_max = np.percentile(data, [5, 95])
print(c_min, c_max)
'''
-1.2757164503037743 0.8703325037861378
'''
data = np.clip(data, c_min, c_max)
'''
[-0.41127091 -1.27571645  0.09598778 -1.19662011 -0.04607188 -0.02745831
  0.23184919  0.85601106  0.58430572  0.8703325 ]
'''
print(data)

この例だけ見てもありがたみを感じないのですが、実際のデータを決定木などにかける時、
ほんの数件のデータだけ極端な外れ値になっていたりすると、
いい感じの範囲にデータを収めることができるので便利です。

また、scikit-learnなどのライブラリのコードを見てみると、
値を 0より大きく1より小さい範囲に収める目的などでも使われています。
ここなど
n以上m以下、ではなくnより大きいmより小さい、で区切る時は便宜上、eps=1e-15のような非常に小さい値を用意して、
[n+eps, m-eps]で代用するようですね。
こういう書き方も参考になります。

numpyのpercentile関数の仕様を確認する

中央値や四分位数を一般化した概念に分位数ってのがあります。
その中でも特にq/100分位数をqパーセンタイルといい、numpyに専用の関数が用意されています。
numpy.percentile

データの可視化や外れ値の除外で使うためにこれの仕様を確認したのでそのメモです。

そもそも僕が何を疑問に思ったのかを説明したほうがいいと思うので、いくつか例を紹介します。

まずわかりやすい例で50パーセンタイル。
これは、奇数個の値があればその中央の値、偶数個の値に対しては、真ん中の二つの値の中点を返します。


import numpy

# 5個の値の3番目の数を返す
data_1 = np.array([3, 12, 3, 7, 10])
print(np.percentile(data_1, 50))  # 7.0

# 6個の値の3番目の数と4番目の数の平均を返す
data_2 = np.array([3, 12, 3, 7, 10, 20])
print(np.percentile(data_2, 50))  # 8.5

同様にして、区切りのいい値がある時のパーセンタイルは非常にわかりやすい。
11個の値があれば、それぞれ順番に 0パーセンタイル, 10パーセンタイル, … 90パーセンタイル, 100パーセンタイルです。


data_3 = np.random.randint(0, 2000, 11)
print(data_3)
# 出力
# [1306  183 1323  266  998 1263 1503 1986  250  305 1397]
for p in range(0, 101, 10):
    print(p, "パーセンタイル・・・", np.percentile(data_3, p))
# 出力
'''
0 パーセンタイル・・・ 183.0
10 パーセンタイル・・・ 250.0
20 パーセンタイル・・・ 266.0
30 パーセンタイル・・・ 305.0
40 パーセンタイル・・・ 998.0
50 パーセンタイル・・・ 1263.0
60 パーセンタイル・・・ 1306.0
70 パーセンタイル・・・ 1323.0
80 パーセンタイル・・・ 1397.0
90 パーセンタイル・・・ 1503.0
100 パーセンタイル・・・ 1986.0
'''

ここまではわかりやすいのですが、自分が疑問に思ったのは、
もっと中途半端なパーセンタイルです。

(例)この出力の40.16ってどうやって算出された?


data_4 = np.array([15, 52, 100, 73, 102])
print(np.percentile(data_4, 17))
#  出力
# 40.16

この疑問放置したままなのが気持ち悪かったので、
これまでパーセンタイルや四分位数、そしてこれらを使う箱ひげ図などを使わなかったのですが、
とあるタスクの中で箱ひげ図を使いたくなったのでこの機会に仕様を確認しました。

といっても、numpyの該当ページにもNote.として記されていますし、
wikipediaにも普通に載ってます。
分位数
あと、pを1刻みで動かして適当なデータに対してパーセンタイル算出してプロットしたら明快にわかりました。

要は、中途半端な値に対しては、隣接の2つの値を線形補完して求めるそうです。
上の例で言えば、
15が0パーセンタイル、52が25パーセンタイルなので、17パーセンタイルは
$(52-15)*17/25+15=40.16$ と計算されています。
仕様がわかったのでこれからはバシバシ使おう。

機械学習を行う時、異常時の前処理として、1〜99パーセンタイルでクリップすると有効なことがあるという話を最近聞いたので、
それも試してみたいです。

pandasでクロス集計

テーブルデータを分析する時、2つの列の値の関係を調べたいことが良くあります。
対象の列がどちらも連続値をとるのであれば、散布図を書きますが、
これが連続値ではなく性別や都道府県、世代、アンケートのn段階評価などの時は、散布図はあまり適しません。

このような場合は、クロス集計を使うと便利です。
(wikipediaでは分割表と呼んでいるようです。)

pandasにクロス集計を行う専用の関数があるのでそれを紹介します。
(職場ではcsvに書き出してtableauに読み込ませてやることも多いのですが、pythonで完結できればそれはそれ便利です。)

ドキュメントはこちら。
pandas.crosstab

aggfuncに値を渡せば、数え上げ以外にも色々できることがあり、pivot_table的な使い方もできるのですが、
とりあえずダミーデータを用意してデータの数え上げでクロス表を作ってみましょう。
ABテストなどで頻繁に使いますからね。


import pandas as pd
import numpy as np
# ダミーデータ生成
df = pd.DataFrame(
        {
            "data1": np.random.randint(1, 6, 100),
            "data2": np.random.choice(["foo", "bar"], 100),
        }
    )
# サンプル(5件)
print(df.sample(5))
# 出力
'''
    data1 data2
71      2   bar
13      1   foo
7       5   bar
16      1   bar
56      2   foo
'''

# クロス集計
cross_df = pd.crosstab(df.data2, df.data1)
print(cross_df)
# 出力
'''
data1   1   2   3   4  5
data2                   
bar    10  13   8  12  8
foo    12   9  11   9  8
'''

とても楽にクロス集計ができました。
SQLでこれを実装しようとするとかなりの手間なので、地味に重宝しています。

完全に蛇足ですが、crosstabを知らなかった頃はこういうコードを書いていました。
indexとcolumnsだけ指定して値が全部0のデータフレームを作って、
元のデータの値を順番に数え上げています。
今思えばcrosstab知らなくても、もう少しまともな書き方がありそうです。


cross_df_2 = pd.DataFrame(
        0,
        index=set(df.data2),
        columns=set(df.data1),
    )
for _, row in df.iterrows():
    cross_df_2.loc[row.data2, row.data1] += 1    
print(cross_df_2)
# 出力
'''
     1   2   3   4  5
foo  12   9  11   9  8
bar  10  13   8  12  8
'''

pythonで一般化線形モデル(GLM) – ポアソン回帰

以前の記事で、statsmodelsを使って線形回帰をやったので、今回はポアソン回帰をやってみます。
参考:statsmodelsで重回帰分析

データは久保拓弥先生の、データ解析のための統計モデリング入門 (通称緑本)の第3章から拝借し、
本に載っているのと同じ結果を得ることを目指します。
ちなみに本のコードはRで実装されています。
Rは勉強中なので写経はしましたがそれはそれとして、
今回はいつも使っているpythonのstatsmodelsでやってみます。

データはこちらのサポートページからダウンロードできます。
生態学のデータ解析 – 本/データ解析のための統計モデリング入門
3章の data3a.csv を保存しておきましょう。

このデータは、架空の植物の種子の数に関するもので、
種子の数が$y_i$, 植物のサイズが$x_i$, 施肥処理を行ったかどうかが$f_i$列に格納されています。

そして、種子の個数$y_i$が期待値 $\lambda_i=\exp(\beta_1+\beta_2x_i+\beta_3f)$のポアソン分布に従うと仮定した場合の尤度を最大化するように係数を決定します。
($f$は施肥処理を行ったら1, 行ってない場合は0)

早速やってみましょう。
ドキュメントはこのあたりが参考になります。
statsmodels.genmod.generalized_linear_model.GLM
Model Family に Poisson を指定すると、リンク関数は自動的にlogが選ばれます。


import pandas as pd
import statsmodels.api as sm

# データの読み込み
# ファイルの取得元
# http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html
df = pd.read_csv("./data3a.csv")

# f 列の値を Tならば 1 , その他(Cのみ)ならば0に符号化
df["fT"] = (df.f == "T").astype(int)

y = df.y
X = df[["x", "fT"]]
# 定数列(1)を作成
X = sm.add_constant(X)

# モデル生成と学習
model = sm.GLM(y, X, family=sm.families.Poisson())
result = model.fit()

# 結果出力
print(result.summary())

# 以下出力結果
                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                Sun, 21 Apr 2019   Deviance:                       84.808
Time:                        16:36:24   Pearson chi2:                     83.8
No. Iterations:                     4   Covariance Type:             nonrobust
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.2631      0.370      3.417      0.001       0.539       1.988
x              0.0801      0.037      2.162      0.031       0.007       0.153
fT            -0.0320      0.074     -0.430      0.667      -0.178       0.114
==============================================================================

本のP58に載っているのと全く同じ係数を得ることができました。

ついでに実際のデータとこのモデルから予測される種子数をプロットして可視化してみます。


import matplotlib.pyplot as plt
import numpy as np

# プロット用のデータ作成
xx = np.linspace(7, 13, 101)
yy0 = np.exp(result.params["const"] + result.params["x"]*xx)
yy1 = np.exp(result.params["const"] + result.params["x"]*xx + result.params["fT"])

# 可視化
fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel("個体のサイズ")
ax.set_ylabel("種子の個数")
ax.scatter(X[X.fT == 0]["x"], y[X.fT == 0].values, marker="x", label="肥料無し")
ax.scatter(X[X.fT == 1]["x"], y[X.fT == 1].values, alpha=0.7, label="肥料有り")
ax.plot(xx, yy0, linestyle="--",  label="肥料無し")
ax.plot(xx, yy1, label="肥料有り")
plt.legend()
plt.show()

出力結果がこちら。