Python: SciPyを使った仮説検定

前回はSciPyを使った推定をとりあげましたが、今回は正規分布に関する仮説検定をSciPyを使って行います。

サンプルには前回同様、Irisデータセットを使いますので、UCIのWebサイトからダウンロードして品種ごと(setosa, virginica, versicolor)に分けておきます。それぞれの標本数は50です。

import pandas as pd
import math
from scipy import stats

iris = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', 
                 header=None, names=['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'class'])

setosa = iris[iris['class'] == 'Iris-setosa'] # Iris setosa(ヒオウギアヤメ)
virginica = iris[iris['class'] == 'Iris-virginica'] # Iris virginica
versicolor = iris[iris['class'] == 'Iris-versicolor'] # Iris versicolor

また、赤本の練習問題からも一部抜粋しております。

統計学入門 (基礎統計学?)

仮説検定のプロセス

最初に仮説検定のプロセスを確認します。

1. 帰無仮説と対立仮説を立てる
   例えば、setosaのがく片の長さ(sepal_width)は母平均は3.5であるという仮設(帰無仮説)を立てます。ここでは とします。
   対して、平均3.5ではないという仮設(対立仮説)を立てます。同じく とします。2つの仮設で大小は関係ない(ただ異なっていることを示す)ため、両側検定となります。

2. 有意水準を決める
   帰無仮説を棄却するに足るための水準を決めます。
   例えば、両側検定で有意水準95%とした場合、[0.025, 0.975]が を支持する区間(採択域)となり、この区間から外れた場合に が棄却(つまり が支持)される。

3. 統計量を計算する
   例えば、正規分布の平均であればz値を計算することで、標準正規分布上にプロットできるようになる。
   そうすると「何%の確率で が得られるのか」というのが確率密度が計算できる。

4. 有意水準と照らして判定する
   得られた確率が採択域から外れていれば を棄却( を支持)、採択域に入っていれば を棄却しないと判定する。

母平均の仮説検定

大標本(n>30)の場合は正規分布、小標本(n≦30)の場合はt分布に従うと仮定します。ただし、nが大きくなればt分布は正規分布に近づくため、nに関わらず小標本の求め方でOKだと思います。

大標本 (n>30)

setosaのがく片の長さ(sepal_width)の母平均について、帰無仮説 、対立仮説 とし、有意水準95%で仮説検定します。

z統計量を計算します。未知の母分散 については、大標本のため標本不偏分散 で推定しています。

$$ z=\frac{\bar{X} - \mu_0}{\sqrt{\frac{\sigma^2}{n}}}=\frac{\bar{X} - \mu_0}{\sqrt{\frac{s^2}{n}}} $$

z値として-1.52が得られましたので、標準正規分布における累積確率([-∞, -1.52]の区間に入る確率)を累積分布関数cdfで計算します。
結果は0.0640でしたので、95%有意水準の採択域[0.025, 0.975]に入っているため、 は棄却できない、と判定されます。

mu_0 = 3.5
z = (setosa['sepal_width'].mean() - mu_0) / math.sqrt(setosa['sepal_width'].var() / len(setosa))
print(z) # -1.5217596660094042

stats.norm.cdf(z)
# 0.06403465560380407

小標本 (n≦30)

小標本の場合、正規分布ではなくt分布で推定される。
そのため、まずは以下でt統計量を計算する。そして、自由度n-1のt分布上で累積確率を計算する。

$$ t=\frac{\bar{X} - \mu_0}{\sqrt{\frac{s^2}{n}}} $$

n=10として、同じくsetosaのがく片の長さの母平均について、帰無仮説 、対立仮説 において、有意水準95%で仮説検定した実装が以下です。

SciPyではttest_1samp()に、データと帰無仮説(popmean)を渡すと、1発で計算できます。
返り値はタプルで、1つ目がt統計量、2つ目がp値となっています。p値は、帰無仮説が成り立つと仮定した場合に、t統計量以上の値(ここでは|-1.96|=1.96以上の値)が得られる確率を表します。p値は0.082となっていますので、95%有意水準に照らして は棄却できないという判定になります。

n = 10
mu_0 = 3.5

stats.ttest_1samp(setosa.iloc[:n]['sepal_width'], popmean=mu_0)
# Ttest_1sampResult(statistic=-1.9562349357055535, pvalue=0.08213989920574083)

母平均の差の仮説検定

2つの標本XとYの平均の差について仮説検定します。等分散を仮定できるかどうかで、求め方が異なるので注意が必要。

ここでは赤本の練習問題12.2のデータを使います。同年代・同学歴・同職種の男女各10人の賃金(万円)を表すを架空の例です。
帰無仮説 、対立仮説 とし、有意水準95%で検定します。

payment_male = [15.4, 18.3, 16.5, 17.4, 18.9, 17.2, 15.0, 15.7, 17.9, 16.5]
payment_female = [14.2, 15.9, 16.0, 14.0, 17.0, 13.8, 15.2, 14.5, 15.0, 14.4]

等分散が仮定できる

2つの標本間で母分散が等しいと仮定できる場合、2つの標本の平均差 の分散を以下で計算して、2標本のt統計量を求めて仮説検定します。

$$ s^2=\frac{(n_X-1)s_X^2+(n_Y-1)s_Y^2}{n_X+n_Y-2} $$

$$ t=\frac{\bar{X}-\bar{Y}}{s\sqrt{\frac{1}{n_X}+\frac{1}{n_Y}}} $$

ナイーブに上の式でも計算できますが、ttest_ind()を使うと、2標本のデータを渡すだけでt統計量とp値を計算できます。

t統計量は3.61、p値は0.0020で有意水準(0.05/2=0.025)を大きく下回るため、帰無仮説 は棄却されます。

stats.ttest_ind(payment_male, payment_female)
# Ttest_indResult(statistic=3.606503775109937, pvalue=0.002017792820428875)

等分散を仮定できない

母分散が等しいと仮定できない場合、ウェルチの近似法でt統計量を近似します。

等分散でない平均の差は、以下の計算で得られた に最も近い整数 を自由度とするt分布に従うことを利用して検定します(ウェルチのt検定と言います)。

$$ \upsilon=\frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2})^2}{n_2-1}} $$

ttest_ind()のオプション引数equal_varをFalseとすることで、ウェルチのt検定が行われます(デフォルトではTrueで、これは上の通り等分散を仮定したt検定となります)。
賃金のデータに適用すると、t統計量は変わらず3.61でした。実際に を計算してみると17.17で、 は17となり、上で求めた自由度18のt統計量と異なるはずですが、そうなってないです(四捨五入ではなく切り上げで補正されているのでしょうか)。
一方で、t統計量は同じにもかかわらず、p値が少し異なっているのも気になります(何か誤解しているかもしれませんので、注意してください)。

stats.ttest_ind(payment_1, payment_2, equal_var=False)
# Ttest_indResult(statistic=3.606503775109937, pvalue=0.0021488395200926777)

標本間で対応がある場合の平均差の検定

同じ個体で比較する場合、t統計量を以下で求めます。例えば、高血圧の人にある薬を投薬してその前後の血圧の平均の差を検定する、などです。このt統計量を使った検定を、対応ありのt検定と言います。
、 とすると以下で計算されます。

$$ s^2=\frac{1}{n-1}\sum_i^n d_i^2 $$

$$ t=\frac{\bar{d}-\mu_0}{\sqrt{\frac{s^2}{n}}} $$

ここでは降圧剤を想定した架空の例として、投薬前の血圧値(before)と投薬後の血圧値(after)をデータに対して検定してみます。

SciPyでは、ttest_ind()の代わりに、ttest_rel()を使います。
t統計量1.56、p値0.14が得られました。

before = [141, 139, 140, 148, 138, 133, 142, 134, 139, 132, 140, 138]
after = [143, 134, 136, 140, 141, 133, 145, 128, 137, 133, 140, 134]

stats.ttest_rel(before, after)
# Ttest_relResult(statistic=1.560009076442849, pvalue=0.147047493905328)

母分散の仮説検定

母分散 の仮説検定は、 検定統計量が自由度(n-1)の 分布に従うことを利用して行います。

帰無仮説 としたときの は以下で計算されます。

$$ \chi^2=\frac{(n-1)s^2}{\sigma_0^2} $$

有意水準95%(両側検定)の場合、この 統計量が 分布の[ ]に含まれるかどうかは計算し、含まれない場合は帰無仮説を棄却します。
例えば、上の例の男性の賃金において、帰無仮説を とした時の 統計量は2.49と計算されます。自由度9における確率は0.019のため、帰無仮説は棄却されます。

sigma2_0 = 6.0

chi2_value = (len(payment_male) - 1) * pd.Series(payment_male).var() / sigma2_0
print(chi2_value) # 2.485999999999999

stats.chi2.cdf(chi2_value, len(payment_male) - 1)
# 0.018742372504659505

この 検定は、適合度の検定や独立性の検定などに発展されます。

適合度の検定

想定していた確率分布から得られる度数(期待度数)と、実際に得られた度数(観測度数)が、一致していると言えるのかどうかを検定するのが適合度の検定です。
例えば、サイコロを50回振って出た目をカウントし、そのサイコロが歪んでないか(均等に各目が1/6の確率で現れているのか)を確かめる、といった問いに答える方法です。

k個のカテゴリ(サイコロであればk=6)について、理論確率 から得られる各カテゴリの期待度数を 、観測度数を とすると、以下の式で 統計量が計算できます。この統計量は自由度k-1の 分布に従います。

$$ \chi^2=\sum_i=1^k \frac{(f_i-np_i)^2}{np_i} $$

サイコロを50回振った時の観測度数が以下の表で得られたとします。公平なサイコロであれば、期待度数は全て50/6です。このときに理論確率(全ての目が1/6の確率で出る)に対する帰無仮説 [\tex: H_0: p_1=1/6, ... , p_6=1/6 ] を立て、有意水準95%で検定します。

サイコロの目	1	2	3	4	5	6
観測度数	10	7	8	11	6	8
期待度数	50/6	50/6	50/6	50/6	50/6	50/6

SciPyでは、chisquare()の引数f_obsに観測度数、f_expに期待度数を与えることで、 統計値とp値を計算できます。
p値は0.84( は2.08)のため、帰無仮説は棄却できません。

stats.chisquare([10, 7, 8, 11, 6, 8], 
                f_exp=[50/6, 50/6, 50/6, 50/6, 50/6, 50/6])
# Power_divergenceResult(statistic=2.08, pvalue=0.8379684711477735)

独立性の検定

代数学と解析学の成績の関連など、2つの属性間の関連の有無についても [\tex: \chi^2 ] 検定で確認できます。独立性の検定と呼ばれます。

k個のカテゴリからなる属性Aと、m個のカテゴリからなる属性Bについて、観測度数を列挙します。

		...		計
		...
...	...	...	...	...
		...
計		...		n

帰無仮説 を立て、以下で 統計量を計算します。この値は自由度(k-1)(m-1)の 分布に従います。

$$ \chi^2=\sum_i \sum_j \frac{(nf_{ij}-f_{i \cdot}f_{\cdot j})^2}{nf_{i \cdot}f_{\cdot j}} $$

ここで、喫煙習慣の有無と6年後の生存/死亡を調べた表を用いて、2つが独立であることを帰無仮説として、有意水準95%で検定します(赤本の練習問題12.6からの抜粋)。

	非喫煙者(non-smoker)	喫煙者(smoker)	計
生存(alive)	950	348	1298
死亡(dead)	117	54	171
計	1067	402	1469

SciPyではchi2_contingency()に表で渡すことで計算できます。
返り値は4つの要素を持つタプルとなっており、( 統計量, p値, 自由度, 完全に独立の場合の期待度数)が格納されています。p値が0.22のため、帰無仮説が棄却できず、2つの属性が独立していないと判定できません。

smoking_df = pd.DataFrame({
    'non-smoker': [950, 117],
    'smoker': [348, 54]
}, index=['alive', 'dead'])

stats.chi2_contingency(smoking_df)
# (1.496886561580958,
#  0.22115103741056918,
#  1,
#  array([[942.79509871, 355.20490129],
#         [124.20490129,  46.79509871]]))

まとめ

上で紹介した以外にも、母集団の比率の検定や二項検定など、様々な検定がありますので随時触れていきたいと思います。

Python: SciPyを使った統計的推定

SciPyを使って統計的推定を行います。

Irisデータセットを使いますので、UCIのWebサイトからダウンロードし、品種ごと(setosa, virginica, versicolor)に分けておきます。それぞれの標本数は50です。

import pandas as pd
import math
from scipy import stats

iris = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', 
                 header=None, names=['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'class'])

setosa = iris[iris['class'] == 'Iris-setosa'] # Iris setosa(ヒオウギアヤメ)
virginica = iris[iris['class'] == 'Iris-virginica'] # Iris virginica
versicolor = iris[iris['class'] == 'Iris-versicolor'] # Iris versicolor

点推定

母平均 と母分散 を以下の式で点推定し、推定量 と を得ます。

• 母平均は、標本平均で推定
• 母分散は、不偏分散で推定

$$ \hat{\mu}=\bar{X}=\frac{1}{n} \sum_{i=1}^n X_i $$

$$ \hat{\sigma^2}=s^2=\frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X}) $$

PandasのSeriesでは、mean()とvar()でそれぞれ計算できます。

setosa['sepal_width'].mean()
# 3.418

setosa['sepal_width'].var()
# 0.1451795918367347

区間推定

平均

中心極限定理によって、nが大きくなると標本平均 は正規分布 に近づきます。
これを利用し、平均 の100(1-α)%の信頼区間は以下で計算します。zは標準正規分布N(0,1)に従う値で、95%信頼区間の場合、[ ]=[ ]となります。

$$ \bar{X}-z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma^2}{n}} < \mu < \bar{X}+z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma^2}{n}} $$

ただし、通常では分散 が未知のため、上の式そのままでは求められず、工夫が必要となります。

大標本

標本数が多い大標本(目安としてはn>30)場合、中心極限定理により とみなせるため、上の式の をsで置き換えればOKです。

$$ \bar{X}-z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} < \mu < \bar{X}+z_{\frac{\alpha}{2}}\sqrt{\frac{s^2}{n}} $$

scipy.statsを使うと、こんな感じで実装できます。
stats.norm.interval()で、平均(loc)・標準偏差(scale)の正規分布でalpha×100%となる範囲を、中央値を中心として取得できます。例えば、stats.norm.interval(alpha=0.95, loc=0, scale=1)なら[-1.96, 1.96]となります。

stats.norm.interval(alpha=0.95, 
                    loc=setosa['sepal_width'].mean(), 
                    scale=math.sqrt(setosa['sepal_width'].var() / len(setosa)))
# (3.3123873659408125, 3.5236126340591887)

小標本

標本数が少ない(n≦30)場合、 は正規分布 よりも、さらに裾野が広い、自由度n-1のt分布の方が当てはまりが良いことが知られています。

$$ \bar{X}-t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} < \mu < \bar{X}+t_{\frac{\alpha}{2}}(n-1)\sqrt{\frac{s^2}{n}} $$

正規分布と同じように、stats.t.interval()を使うことで、自由度(df)・平均(loc)・標準偏差(scale)となる範囲を計算できます。

n = 10
stats.t.interval(alpha=0.95,
                 df=n-1,
                 loc=setosa.iloc[:n]['sepal_width'].mean(), 
                 scale=math.sqrt(setosa.iloc[:n]['sepal_width'].var() / n))
# (3.0902871970662513, 3.529712802933748)

n=10において、正規分布では[3.12, 3.50]で区間推定される一方、t分布では[3.09, 3.53]が得られており、t分布の方が少し裾野が広くなっていることがわかります。
nが大きくなるほど、t分布は正規分布に近づくため、区間推定で得られる範囲はほぼ一致します。そのため、nに関わらず、最初からt分布を使うことが多いようです。

平均の差の区間推定

setosaとverginicaのがく片幅の平均値の差の区間推定について考えます。

平均の差は以下で推定できます(2つの正規分布の差 の分布より導出できます)。こちらも通常は分散 が未知のため、工夫が必要となります。

$$ \bar{X_1}-\bar{X_2}-z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}} < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+z_{\frac{\alpha}{2}}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}} $$

大標本

大標本(n>30)であれば、 として、以下の式で求められます。

$$ \bar{X_1}-\bar{X_2}-z_{\frac{\alpha}{2}}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}} < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+z_{\frac{\alpha}{2}}\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}} $$

標本数n=50のsetosaとvirginicaのがく片幅の平均差を、95%信頼区間で推定します。

stats.norm.interval(alpha=0.95,
                   loc=setosa['sepal_width'].mean()-virginica['sepal_width'].mean(),
                   scale=math.sqrt(setosa['sepal_width'].var()/len(setosa) + virginica['sepal_width'].var()/len(virginica)))
# (0.30563607258957026, 0.5823639274104314)

小標本 (等分散を仮定できる)

小標本(n≦30)の場合、さらに等分散を仮定できる場合と仮定できない場合で、求め方が異なります。

等分散を仮定できるなら、t分布により計算できます。

$$ \bar{X_1}-\bar{X_2}-t_{\frac{\alpha}{2}}(n_1+n_2-2)\sqrt{ \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}(\frac{1}{n_1}+\frac{1}{n_2}) } < \mu_1-\mu_2 < \bar{X_1}-\bar{X_2}+t_{\frac{\alpha}{2}}(n_1+n_2-2)\sqrt{ \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}(\frac{1}{n_1}+\frac{1}{n_2}) } $$

setosa(標本数 、標本平均 )とvirginica(標本数 、標本平均 )のがく片の幅の平均差について、母平均の95%信頼区間を求めます。
不偏分散は概ね同じ(0.94≒0.96)ため、等分散を仮定し、上の式で計算して求めます。

n_1 = 10 # setosaの標本数
n_2 = 15 # virginicaの標本数

s_1 = setosa.iloc[:n_1]['sepal_width'].var()
s_2 = virginica.iloc[:n_2]['sepal_width'].var()

# 分散が近そうなことを確認
print(s_1) # 0.09433333333333332
print(s_2) # 0.09638095238095239

stats.t.interval(alpha=0.95,
                 df=n_1+n_2-2,
                 loc=setosa.iloc[:n_1]['sepal_width'].mean()-virginica.iloc[:n_2]['sepal_width'].mean(),
                 scale=math.sqrt((((n_1-1)*s_1 + (n_2-1)*s_2)/(n_1+n_2-2)) * (1/n_1 + 1/n_2)))
# (0.14223996284479556, 0.6644267038218702)

小標本 (等分散を仮定できない)

等分散でない場合、t統計量を近似するウェルチの近似法を使います。

ウェルチの近似法では、t統計量が で近似します。
下式で に最も近い整数 をパラメータとする に従うことを利用する。
あとは上の式で信頼区間を推定できます。

$$ \upsilon=\frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2})^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2})^2}{n_2-1}} $$

がく片長(sepal_length[cm])について、setosa(標本数 , 標本平均 )とvirginica(標本数 , 標本平均 )の平均の差を推定します。分散はそれぞれ0.0849、0.505のため、等分散を仮定できないものとなっています。
t統計量を で近似して95%信頼区間を計算すると、[-2.07, -1.13]となります。

n_1 = 10 # setosaの標本数
n_2 = 15 # virginicaの標本数

s_1 = setosa.iloc[:n_1]['sepal_length'].var()
s_2 = virginica.iloc[:n_2]['sepal_length'].var()

# 分散が近いとは言えない
print(s_1) # 0.08488888888888892
print(s_2) # 0.5054285714285711

upsilon = (s_1/n_1 + s_2/n_2) / ((s_1/n_1)**2/(n_1-1) + (s_2/n_2)**2/(n_2-1))
t = round(upsilon)

stats.t.interval(alpha=0.95,
                 df=t,
                 loc=setosa.iloc[:n_1]['sepal_length'].mean()-virginica.iloc[:n_2]['sepal_length'].mean(),
                 scale=math.sqrt((((n_1-1)*s_1 + (n_2-1)*s_2)/(n_1+n_2-2)) * (1/n_1 + 1/n_2)))
# (-2.0683594356953647, -1.131640564304636)

分散

母分散 の区間推定では、カイ二乗分布を使います。

母分散 と不偏分散 自由度n-1のカイ二乗分布に従います。

$$ \frac{1}{\sigma^2}\sum_i^n (X_i-\bar{X})^2 = \frac{(n-1)s^2}{\sigma^2} \sim \chi^2(n-1) $$

これを について整理すると、100(1-α)%の信頼区間を求める以下の式が導出できます。

$$ \frac{(n-1)s^2}{\chi^2_\frac{\alpha}{2}(n-1)} < \sigma^2 < \frac{(n-1)s^2}{\chi^2_{1-\frac{\alpha}{2}}(n-1)} $$

setosaのがく片幅の95%信頼区間を推定します。
正規分布やt分布同様、stats.chi2.interval()で、自由度(df)のカイ2乗分布でalpha×100%となる範囲を取得できますので、あとは上の式に当てはめると計算できます。
結果は[0.101, 0.223]となりました。

df = len(setosa) - 1 # 自由度n-1
chi2_025, chi2_975 = stats.chi2.interval(alpha=0.95, df=df)
(df*setosa['sepal_width'].var() / chi2_975, df*setosa['sepal_width'].var() / chi2_025)
# (0.10130383788745644, 0.225441889805552)

時系列データで使うPandas小技集

時系列データを扱うにあたって役に立った、Pandasのテクニックを紹介します。

• 文字列型のSeriesから日時型・日付型のSeriesへ変換する
• 日付に欠測値を含むデータを日毎に集計する
• 累積和を計算する

今回の例に使う時系列データは以下です。
ある商品の4/1〜4/3の購入履歴をイメージしてください。ユーザ(user_id)の購入日時(timestamp)と購入数(item_count)が入ったDataFrameとなっています。

import pandas as pd

df = pd.DataFrame(...)

print(df.info())
# <class 'pandas.core.frame.DataFrame'>
# RangeIndex: 5 entries, 0 to 4
# Data columns (total 3 columns):
# timestamp     5 non-null object
# user_id       5 non-null object
# item_count    5 non-null int64
# dtypes: int64(1), object(2)
# memory usage: 200.0+ bytes

日時型や日付型への変換

timestampカラムは文字列になっていますので、最初に扱いやすい日時型へ変換する必要があります。

Seriesの型変換にはastype()がしばしば使われますが、文字列から日時へ変換することはできません。

代わりに、pandas.to_datetime()を使うと、Series型の各要素を日時型へ一括変換できます。

変換後は、Pythonのdatetime型ではなく、Pandasの提供するTimestamp型となります。このTimestamp型は、いろいろなフォーマットの日時文字列に対応しているため、上のような'/'区切りの日付でも問題ありません。

df['timestamp'] = pd.to_datetime(df['timestamp'])

print(df.info())
# <class 'pandas.core.frame.DataFrame'>
# RangeIndex: 5 entries, 0 to 4
# Data columns (total 3 columns):
# timestamp     5 non-null datetime64[ns]
# user_id       5 non-null object
# item_count    5 non-null int64
# dtypes: datetime64[ns](1), int64(1), object(1)
# memory usage: 200.0+ bytes

timestampは日時ですが、日単位で集計していきますので、日付型のカラム"date"を追加します。

Series型にはdtというアクセサが提供されており、Timestamp型を含む日時型の要素から、日付や時刻のみの要素へ一括変換できます。

df['date'] = df['timestamp'].dt.date

print(type(df.loc[0, 'date']))
# datetime.date

日付に欠測があるデータを日毎に集計する

4/1〜4/3の3日間の、各ユーザの購入数を集計したいとします。

user_idとdateでgroupbyすれば良さそうですが、それでは購入していない日の購入数(つまり0)が取得できません。

# 以下では購入していない日('0001'は4/2、'0002'は4/3)の値が取得できない
pd.DataFrame(df.groupby(['user_id', 'date']).sum()['item_count']).reset_index()

こうした場合、日付とユーザIDのみからなるDataFrameを最初に作成し、上の集計済みDataFrameと左外部結合して、欠測値を0埋める、という3段階に分けて行います。

最初に日付+ユーザIDのDataFrameの作成です。

pandas.date_range()を使うと、指定日時('2018-04-01')から一定の間隔(ここではfreq='D'のため、1日ごと)の日時型のSeriesを作ってくれます。periodsは繰り返し回数を指定し、ここでは3回となっています。

date_df = pd.DataFrame(pd.date_range('2018-04-01', periods=3, freq='D').date, columns=['date'])

ユーザIDのDataFrameも作ります。unique()で重複を除外したユーザIDの一覧を取得します。

user_id_df = pd.DataFrame(df['user_id'].unique(), columns=['user_id'])

この2つのDataFrameをクロスジョインで結合します。
merge()でサポートされるのはキーによる結合だけですので、ここでは2つのテーブルの全ての行が同じ値0となるカラムkeyを追加し、それをキーとして結合しています。これにより、一方のDataFrameの各行を、もう一方のDataFrameの全ての行と結合できます。結合後、カラムkeyは不要ですので、drop()で削除しておきます。

date_df['key'] = 0
user_id_df['key'] = 0

tmp_df = date_df.merge(user_id_df, on='key').drop('key', axis=1)

次に、上で作成したDataFrameと集計値のDataFrameをmerge()で左外部結合します。
購入していない日はNaNの値が入っていることがわかります。

sum_df = pd.DataFrame(df.groupby(['user_id', 'date']).sum()['item_count']).reset_index()
sum_df = tmp_df.merge(sum_df, on=['date', 'user_id'], how='left')

最後に、fillna()で0埋めします。

sum_df = sum_df.fillna(0)

累積和を計算する

購入数の累積和を計算する場合、cumsum()が使えます。

以下では、日毎の購入個数の累積和を計算しています。
DataFrameは累積和を計算したい順にソートする必要があります。ここでは、インデックス(日付)順にソートすることで、4/2は4/1までの購入数+4/2の購入数、4/3は4/2までの購入数+4/3の購入数となるようにしています。

daily_sum_df = sum_df.groupby('date', group_keys=False).sum()['item_count']
daily_sum_df.sort_index()

daily_sum_df.cumsum()
# date
# 2018-04-01    3.0
# 2018-04-02    6.0
# 2018-04-03    9.0
# Name: item_count, dtype: float64

さらに、groupbyと組み合わせると、ユーザごとの累積和も計算できます。

sum_df['user_cumsum'] = sum_df.groupby('user_id').cumsum()['item_count']
sum_df = sum_df.sort_values(['date', 'user_id'])