C# クエリストリング(?var=hoge&...)を作る

C#でクエリストリングを作る方法の備忘録です。

クエリストリングは、URLのパスの後ろに?変数名1=値1&変数名2=値2&...といった形で任意の値が渡される文字列です。例えば以下のような文字列です。

http://ohke.hateblo.jp/search?q=Python&page=1504879200

LINQで文字列を連結することもできますが、' '(スペース)や'&'(アンパサンド)などの特殊記号をエスケープする必要があったりするため、自前での実装は避けたほうがいいです。

かわりにSystem.Web.HttpUtility.ParseQueryString.aspx)メソッドを使うと、安全かつ容易にクエリストリングを生成できます。

  • ' 'や'&'などはパーセントエンコーディングされています
var query = System.Web.HttpUtility.ParseQueryString("");

query.Add("name", "Tanaka, Ichiro");
query.Add("childs", "jiro&saburo");

var queryString = query.ToString();

Console.WriteLine(queryString);
// name=Tanaka%2c+Ichiro&childs=jiro%26saburo

生成されたクエリストリングをSystem.UriBuilder.aspx)クラスに渡せば、URIも一発生成できます。

var uriBuilder = new System.UriBuilder("example.host.name") {
    Query = query.ToString()
};

Console.WriteLine(uriBuilder.Uri);
// http://example.host.name/?name=Tanaka%2c+Ichiro&childs=jiro%26saburo

Python: 回帰モデルで市場反応分析

これまで緑本などで学んできた統計モデルを、マーケティングに応用するための勉強を行っています。

今回は市場反応分析を線形回帰モデルとポアソン回帰モデルで行います。
市場反応分析に関する理論や使用するデータは、マーケティングの統計モデル (統計解析スタンダード)の第3章「消費者の市場反応のモデル化」を参考にさせていただいてます。

マーケティングの統計モデル (統計解析スタンダード)

実装はこちらです。

statistical_modeling_for_marketing/chapter3.ipynb at master · ohke/statistical_modeling_for_marketing · GitHub

市場反応分析

マーケティング施策は、顧客の行動(購入や予約など)の変化を促すために行われます。

  • 例えば、ある商品を値下げすると、その商品の販売個数が増える
  • 例えば、ある商品を山積み陳列すると、その商品の販売個数は増えるが、類似商品の販売個数は減る ...など

施策を効率的に実施するためには、上のように施策の効果を定性的に捉えるだけではなく、定量化する必要があります。

  • ある商品をX%値下げすると、その商品の販売個数がY%増える
  • ある商品を山積み陳列すると、その商品の販売個数はX%増えるが、類似商品の販売個数はY%減る

このように、施策を受けた顧客の行動変化を定量的にモデル化することを市場反応分析といいます。

価格と陳列によるPI値の変化を分析

それでは、擬似的なデータを使って、具体的に市場反応分析をしていきます。データはこちらからダウンロードできます。

最初の事例として、2つのブランドの醤油(醤油Aと醤油Bとする)を販売するケースを考えます。

両方の商品について「価格を変更する」「山積み陳列する」というそれぞれ2つの施策が組み合わせて行われています。
上の施策を実施した結果、醤油AのPI値の対数がどのように変化するのかをモデル化します。つまり醤油AのPI値の対数が目的変数となります。

データを概観する

まずはデータを概観します。

  • LogPI_Aが醤油AのPI値の対数
    • PI値とは、その商品の販売点数÷客数で計算される値で、その商品の人気度・支持度を表す指標で、
    • LogPI_Aは、  \log \frac{\mbox{販売点数}\times \mbox{客数}}{1000} で計算されています
  • LogPriceIndex_AとLogPriceIndex_Bはそれぞれ、  \log \frac{\mbox{商品Aの売価}}{\mbox{商品Aの期間最大売価}}  \log \frac{\mbox{商品Bの売価}}{\mbox{商品Bの期間最大売価}} で計算される値で、ここでは価格掛率と呼びます
    • 値が小さいほど安売りされていることを表す
  • Display_AとDisplay_Bは、それぞれの商品が山積みされていれば1となるダミー変数

f:id:ohke:20180324132509p:plain

醤油Aのと価格掛率の相関を見てみます。

  • 醤油Aが安売りされている(価格掛率LogPriceIndex_Aが小さい)と、PI値が大きくなる
  • 醤油Bが安売りされている(価格掛率LogPriceIndex_Bが小さい)と、PI値が大きくなる
    • 醤油Aと醤油Bが競合関係にあることがわかります

f:id:ohke:20180324134021p:plain

山積み陳列の有無による醤油AのPI値の分布を箱ひげ図で見てみます。

  • 醤油Aが山積み陳列される(Display_A=1)と、PI値が大きくなる
  • 醤油Bが山積み陳列される(Display_B=1)と、PI値が小さくなる

f:id:ohke:20180324134518p:plain

2標本t検定で醤油Bの山積み陳列の効果を確認

上の箱ひげ図では、醤油Aの山積み陳列によるPI値の増加ははっきり現れていますが、醤油Bの山積み陳列による影響はそれより弱く、誤差という見方ができるかもしれません。

モデル化の前に、醤油Bの山積み陳列によるPI値の変化が誤差ではない(起きにくい)ことを、念のため検定で確認しておきます。ここでは、価格掛率は一旦無視します。なおこの判断は誤りであることが後にわかります。線形回帰モデルを作る時にそれも見ていきます。

なお、検定では統計学入門 (基礎統計学?)の第11・12章を参考にしています。

平均値の差が有意かどうかをt検定で確認します。2つのグループで分散が異なっている(Display_B=1では0.765、Display_B=0では0.910)ため、ウェルチのt検定を行います。

まずはt統計量を計算します。  s が不偏分散、  n がデータ数です。下付き添字が1ならDisplay_B=1、0ならDisplay_B=0を表します。

$$ t=\frac{\overline{PI_1} - \overline{PI_0}}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0}}} $$

t分布の自由度  \nu は以下の式で近似できます。 \nu から最も近い整数  \nu^{*} を選び、自由度  \nu^{*} のt分布の5%の下側信頼限界の  t_{0.05}(\nu^{*}) 値を算出します。上で求めた値がこの値よりも小さければ(つまり  t \mbox{<} t_{0.05}(\nu^{*}) )、有意差があると判定します。これは左片側検定です。

$$ \nu=\frac{(\frac{s_1^2}{n_1}+\frac{s_0^2}{n_0})^2}{\frac{(s_1^2 / n_1)^2}{n_1-1}+\frac{(s_0^2 / n_0)^2}{n_0-1}} $$

Pythonではこんな感じで計算できます。

m = len(treatment_A)
n = len(control_A)
s_1 = treatment_A.var()
s_2 = control_A.var()

t = (treatment_A.mean() - control_A.mean()) / math.sqrt(s_1**2 / m + s_2**2 / n)
print('t:', t)

# 自由度νを近似的に計算
nu = (s_1 ** 2 / m + s_2**2 / n)**2 / ((s_1**2 / m)**2 / (m-1) + (s_2**2 / n)**2 / (n-1))
nu_star = round(nu) # 整数ν*に丸める

print('-t_0.95(ν*):', -stats.t.ppf(0.95, nu_star))

 t\mbox{<}-t_{0.05}(\nu^{*}) \iff -8.47\mbox{<}-1.65 となるため、グループ間の平均値に有意差があることも確認できました。

線形回帰でモデル化

PI値は連続変数なので、正規分布を仮定した線形回帰でモデル化します。

$$ LogPIA_t=\beta_0+LogPriceIndexA_t\beta_1+LogPriceIndexB_t\beta_2+DisplayA_t\beta_3+DisplayB_t\beta_4+\epsilon_t $$

import statsmodels.api as sm

model = sm.OLS(data['LogPI_A'],
               sm.add_constant(data[['LogPriceIndex_A', 'LogPriceIndex_B', 'Display_A', 'Display_B']]))
result = model.fit()

結果は以下のとおりです。いずれも各施策の効果が係数となっています。

  • LogPriceIndex_Aの係数が負数(-5.02)、Display_Aの係数が正(0.774)となっており、醤油Aの販売にプラスに働いていることが数値化されています
  • LogPriceIndex_Bの傾きは正数(1.93)で、競合商品の値引きの影響を受けることが確認できます

f:id:ohke:20180324152547p:plain

一方でDisplay_Bの傾きがわずかですが正数(0.0374)となっています。95%信頼区間では0を跨いでいるので負になる可能性もありますが、検定とは反した結果となっています。
LogPriceIndex_BとDisplay_Bの相関が高い(-0.641)にも関わらず、LogPriceIndex_Bを無視して検定したため、Display_Bの影響を過大評価してしまいました。安売りする商品を山積み陳列するというのは、相関係数的にも直感的にもよくあるパターンに見えます。

  • 他の条件が異なるグループ間での比較にt検定を使うべきではありません
    • LogPriceIndex_Bも揃えたグループ間で平均値を比較すべきです
  • こうした複数の変数がからむ場合、初めから回帰モデルに落とした方が良いかと思います

モデルにもLogPriceIndex_B×Display_Bの交互作用項を追加する必要がありそうです。
LogPriceIndex_B×Display_Bを交互作用項(LogPriceIndex_Display_B)としてパラメータに追加して線形回帰した結果が以下です。Display_Bの係数は-0.330(95%信頼区間は[-0.550, -0.110])で、Display_B単体では醤油Aの販売に負に影響していることがわかります。

f:id:ohke:20180326083747p:plain

価格と陳列による販売個数の変化を分析

同じく醤油Aと醤油Bを販売する事例について考えますが、今回は販売個数(つまりカウントデータ)のモデル化を行います。

データを概観する

同様にデータを概観します。

  • Sale_Unit_Aは醤油Aの販売個数
  • LogPriceIndex_AとLogPriceIndex_Bは、先の事例と同じ価格掛率ですが、今回は対数化されていません(0以上1以下の値域をとる)
  • Display_AとDisplay_Bも同じく山積み陳列の有無
  • Visitorsは来店者数

f:id:ohke:20180324155858p:plain

販売個数と価格掛率・来店者数の分布をプロットします。

  • 醤油Aが安くなれば販売個数は増え、醤油Bが安くなれば販売個数が減る傾向が見えます
  • 来店者数と販売個数も緩い正相関がありそう

f:id:ohke:20180324160455p:plain

山積み陳列の有無毎に、販売個数を箱ひげ図でプロットします。ここでもやはり、競合商品(醤油B)を山積み陳列した時は醤油Aの販売個数が減る傾向が伺えます。

f:id:ohke:20180324160506p:plain

ポアソン回帰でモデル化

販売個数はカウントデータなので、ポアソン回帰(+対数リンク関数)でモデル化します。ポアソン分布の平均値  \mu は以下で計算します。

  • Visitorsはオフセット(傾きを持たない定数項)として扱っています

$$ \log(\mu_t)=\log(Visitors_t)+\beta_0+PriceIndexA_t\beta_1+PriceIndexB_t\beta_2+DisplayA_t\beta_3+DisplayB_t\beta_4 $$

Pythonでの実装は以下となります。

  • offsetオプションで、オフセット値を指定できます
import statsmodels.formula.api as smf

model = smf.glm('Sale_Unit_A ~ PriceIndex_A + PriceIndex_B + Display_A + Display_B',
                data=data,
                offset=np.log(data['Visitors']),
                family=sm.families.Poisson(link=sm.families.links.log))
result = model.fit()

結果は以下です。今回はDisplay_Bの傾きが負数となっており、わずかながらですが競合商品の山積み陳列が悪影響を与えることがわかります。

f:id:ohke:20180324161921p:plain

Python: statsmodelsでガンマ回帰

これまで読み進めてきた緑本を改めて復習するために、statsmodelsで一般化線形モデル(GLM)を作り、糖尿病に関するデータの回帰分析を行います。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

実装はこちらです。

blog/statsmodels_gamma.ipynb at master · ohke/blog · GitHub

データセット

今回はscikit-learn付属の糖尿病データセットを使います。

import pandas as pd
from sklearn.datasets import load_diabetes

diabetes = load_diabetes()

data = pd.DataFrame(diabetes.data, columns=("age", "sex", "bmi", "map", "tc", "ldl", "hdl", "tch", "ltg", "glu"))
data['y'] = diabetes.target

データは442レコード有り、各レコードは11個の変数を持ちます。
1年後の糖尿病の進行状況が目的変数  y (値が大きいほど進行していることを表す)となっており、それ以外は説明変数です。説明変数は全て平均0で正規化されています。この10個の説明変数から1年後の進行状況  y を予測する、回帰分析のタスクとなります。

f:id:ohke:20180317113506p:plain

各変数の相関係数を見ます。
 y との相関係数では、 bmi (BMI)が最も高く、 ltg (おそらく中性脂肪)が続いています。
また、 age (年齢)は各変数と緩やかに相関があります。加齢とともに tc (総コレステロール)や ltg が増加している傾向は見て取れます。

f:id:ohke:20180317114039p:plain

 y の分布も確認しておきます。平均152、標準偏差77となっています。

f:id:ohke:20180317120129p:plainf:id:ohke:20180317120139p:plain

GLMの設計

それでは  y を予測するモデルをGLMで設計します。

GLMでは、確率分布・リンク関数・線形予測子を決める必要がありました (1/252/12の投稿も参考にしてください)。

  • 0以上で上限の無い連続値のようですので、確率分布にはガンマ分布を使います
  • ガンマ分布ですので、リンク関数には対数リンク関数を使います

ガンマ分布は、以下の確率密度関数で表される0以上の連続確率分布です。
 \Gamma(s) はガンマ関数で、 s が正の整数なら  \Gamma(s)=(s-1)! となります。ガンマ分布の平均は  s\mu 、分散は  s\mu^2 です。

$$ p(y \mid s,\mu)=\frac{\exp(-\frac{y}{\mu})}{\Gamma (s)\mu^s}y^{s-1} $$

線形予測子で指定するパラメータは以下の5パターンで比較します。比較基準にはAICを使います(参考)。

  1.  \beta_1 : パラメータなし(切片のみ)モデル
  2.  \beta_1+\beta_2bmi : bmiモデル
  3.  \beta_1+\beta_2bmi+\beta_3sex : bmi+sexモデル
  4.  \beta_1 + \beta_2ldl+\beta_3age : ldl+ageモデル
  5.  \beta_1 + \beta_2ldl+\beta_3age+\beta_4 \times ldl \times age : ldl+age+ldl×ageモデル(交互作用項有り)

まずはパラメータなしモデルで、ベースラインとなるAICを確認します。

  • AICは5041
  • このモデルの最大尤度は-2519.6 (残差逸脱度は126.8)
# pip install statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf

# y~1とすることで、切片のみモデルを作る
model0 = smf.glm('y ~ 1', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result0 = model0.fit(disp=0)

print('AIC:', result0.aic)
# AIC: 5041.186062839457

result0.summary()

f:id:ohke:20180317123522p:plain

次に、bmiモデルで最尤推定させたところ、AICは4884となりました。パラメータ無しのモデルと比較すると、-157改善されています。

model1 = smf.glm('y ~ bmi', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result1 = model1.fit(disp=0)

print('AIC:', result1.aic)
# AIC: 4883.66124171387

もう1つ、bmi+sexモデルでも最尤推定します。AICは4885となり、bmiモデルよりも大きいため、複雑度が増すだけのモデルと言えます。したがって、bmiモデルの方が良い予測を行えるモデルです。

model2 = smf.glm('y ~ bmi+sex', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result2 = model2.fit(disp=0)

print('AIC:', result2.aic)
# AIC: 4885.365291972465

ここまでの3つのモデルの予測値を比較してみます。

  • パラメータ無し(model0)は、単一の平均値(152)で推定されています
  • bmiモデル(model1)とbmi+sexモデル(model2)はほぼ同じ予測値となっています

f:id:ohke:20180317130045p:plain

次に、ldl+ageモデルとldl+age+ldl×ageモデルから、交互作用項の有無による差を見ていきます。
一般的に、ldl(LDLコレステロール)は加齢に伴って上がる検査値です(相関係数からも緩い正相関が見られます)。つまり、2つは完全に独立したパラメータではありません。そこで交互作用項ldl×ageを導入します。加齢によるLDLコレステロール値の上昇を打ち消す(つまり負の係数となる)ことを期待しています。

2つモデルのAICの差は5.5で、交互作用項の導入による予測精度の向上が期待できます。

model3 = smf.glm('y ~ ldl+age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result3 = model3.fit(disp=0)
print('AIC:', result3.aic)
# AIC: 5019.939140159227

model4 = smf.glm('y ~ ldl+age+ldl:age', data=data, family=sm.families.Gamma(link=sm.families.links.log))
result4 = model4.fit(disp=0)
print('AIC:', result4.aic)
# AIC: 5014.413561464431

ldl+age+ldl×ageの各パラメータを見ると、 ldlとageは正の係数(それぞれ1.49、1.64)となっていますが、 ldl×ageは負の係数(-28.92)となっており、期待通りの結果となっていることがわかります。

f:id:ohke:20180318105207p:plain