PyLearn ST Ex 02 の回答例¶
回答は1つではありませんので、あくまで例です。
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
dat = np.zeros([20, 5])
# 販売数
dat[:, 0] = [10, 20, 30, 40, 50,
60, 70, 80, 90, 100,
110, 120, 130, 140, 150,
160, 170, 180, 190, 200]
# 値段
dat[:, 1] = [5810, 2280, 4820, 1780, 4200,
6890, 3510, 5720, 3290, 3100,
2970, 1680, 3550, 2510, 3220,
1900, 1850, 5300, 2450, 1000]
# 広告費
dat[:, 2] = [650, 200, 490, 130, 520,
140, 390, 290, 590, 90,
520, 210, 650, 80, 430,
140, 550, 200, 450, 100]
# ネガティブな口コミ数
dat[:, 3] = [24, 31, 19, 22, 33,
17, 21, 16, 19, 18,
11, 13, 32, 14, 12,
9, 8, 14, 3, 1]
# ポジティブな口コミ数
dat[:, 4] = [3, 4, 8, 5, 2,
12, 15, 17, 14, 9,
19, 18, 14, 17, 20,
22, 19, 18, 8, 23]
問題1.¶
In [2]:
# 販売数と値段
plt.scatter(dat[:, 1], dat[:, 0], s=30, c='orange')
plt.xlabel("nedan")
plt.ylabel("hanbai suu") #軸名
plt.grid(True)
plt.show()
In [3]:
# 販売数と広告費
plt.scatter(dat[:, 2], dat[:, 0], s=30, c='k')
plt.xlabel("koukoku")
plt.ylabel("hanbai suu") #軸名
plt.grid(True)
plt.show()
In [4]:
# 販売数とネガティブな口コミ数
plt.scatter(dat[:, 3], dat[:, 0], s=30, c='b')
plt.xlabel("negative")
plt.ylabel("hanbai suu") #軸名
plt.grid(True)
plt.show()
In [5]:
# 販売数とポジティブな口コミ数
plt.scatter(dat[:, 4], dat[:, 0], s=30, c='r')
plt.xlabel("positive")
plt.ylabel("hanbai suu") #軸名
plt.grid(True)
plt.show()
考察としては「販売数と関係がありそうな要因は、口コミである。価格はちょっと弱いけど、関係ありそうかな？広告費は関係がないかも？」あたりの察しがつけばokです。散布図は、subplotを使って、同時に4つ描いてもいいです。
問題2.¶
まずは、プログラムのコーディング。以下が書ければ及第点です。
In [6]:
# 及第点の回答
cor01 = st.pearsonr(dat[:,0], dat[:,1])
cor02 = st.pearsonr(dat[:,0], dat[:,2])
cor03 = st.pearsonr(dat[:,0], dat[:,3])
cor04 = st.pearsonr(dat[:,0], dat[:,4])
print("販売数と値段の関係数: ", cor01[0], " p値: ", cor01[1])
print("販売数と広告費の相関係数: ", cor02[0], " p値: ", cor02[1])
print("販売数とN口コミの相関係数: ", cor03[0], " p値: ", cor03[1])
print("販売数とP口コミの相関係数: ", cor04[0], " p値: ", cor04[1])
小数点以下の細かい情報には、正直そんなに意味がないので、数桁で丸めます。ここまでできるとちょっとだけ加点。
In [7]:
# ちょっと加点される回答
cor01 = st.pearsonr(dat[:,0], dat[:,1])
cor02 = st.pearsonr(dat[:,0], dat[:,2])
cor03 = st.pearsonr(dat[:,0], dat[:,3])
cor04 = st.pearsonr(dat[:,0], dat[:,4])
digit1 = 3
digit2 = 5
r1 = np.round(cor01[0], digit1)
r2 = np.round(cor02[0], digit1)
r3 = np.round(cor03[0], digit1)
r4 = np.round(cor04[0], digit1)
p1 = np.round(cor01[1], digit2)
p2 = np.round(cor02[1], digit2)
p3 = np.round(cor03[1], digit2)
p4 = np.round(cor04[1], digit2)
print("販売数と値段の相関係数: ", r1, " p値: ", p1)
print("販売数と広告費の相関係数: ", r2, " p値: ", p2)
print("販売数とN口コミの相関係数: ", r3, " p値: ", p3)
print("販売数とP口コミの相関係数: ", r4, " p値: ", p4)
分析に慣れている人は、以下のように書きます。何も考えずにここまでかけるとベストです。ポイントは、for文を使うということと、有意判定を自動で行う部分です。for文を使わないと、要因数が数十、数百になったとき、コーディングが著しくめんどくさくなります。
In [8]:
# 大変良い回答
head = ["値段", "広告費", "N口コミ数", "P口コミ数"]
r, p = [], []
digit1, digit2 = 3, 5
for i in range(1, len(dat[0, :])):
temp = st.pearsonr(dat[:,0], dat[:,i])
r.append(np.round(temp[0], digit1))
p.append(np.round(temp[1], digit2))
for i in range(len(p)):
if p[i] < 0.01:
res_str = ": 1%水準有意"
elif p[i] < 0.05:
res_str = ": 5%水準有意"
elif p[i] < 0.10:
res_str = ": 10%水準有意"
else:
res_str = "n.s."
print("販売数と", head[i], "の相関係数: ", r[i], " p値: ", p[i], " ", res_str)
結果・考察としては、販売数と値段の相関係数は-0.448（5%水準有意）、販売数と広告費との相関係数は-0.170（n.s.）、販売数とネガティブな口コミ数との相関係数は-0.745（1%水準有意）、販売数とポジティブな口コミ数との相関係数は0.740（1%水準有意）であった。したがって、値段が安いほど、ネガティブな口コミ数が少ないほど、そしてポジティブな口コミ数が多いほど、そのイヤホンは売れることになる。ただし、値段の影響は、口コミほど大きくはない。また、広告費と販売数の関連は認められなかった。 などが書ければokです。捕捉ですが、有意ではない場合、英語で「not significant」略して「n.s.」などと言います。
問題3.¶
In [9]:
Y = dat[:, 0]
# 行サイズはデータ数20、列サイズは「値段、N口コミ、P口コミ」3
# 広告費は使わないことに注意する
X = np.zeros([20, 3])
# データの代入
for i in range(len(X[:, 0])):
X[i, 0] = dat[i, 1]
X[i, 1] = dat[i, 3]
X[i, 2] = dat[i, 4]
# for文ではなく、これでもok
# X[:, 0] = dat[:, 1]
# X[:, 1] = dat[:, 3]
# X[:, 2] = dat[:, 4]
# 切片あり
ModelLR = LinearRegression() # 線形回帰モデルのセット
ModelLR.fit(X, Y) # パラメータ獲得
print('推定モデル（切片あり）: y^ = ', round(ModelLR.coef_[0], 3), ' x1 + ',
round(ModelLR.coef_[1], 3) , ' x2 + ',
round(ModelLR.coef_[2], 3) , ' x3 + ',
round(ModelLR.intercept_, 3))
# 切片なし
ModelLR_nb = LinearRegression(fit_intercept=False) # 線形回帰モデルのセット
ModelLR_nb.fit(X, Y) # パラメータ獲得
print('推定モデル（切片なし）: y^ = ', round(ModelLR_nb.coef_[0], 3), ' x1 + ',
round(ModelLR_nb.coef_[1], 3) , ' x2 + ',
round(ModelLR_nb.coef_[2], 3) , ' x3 + ',
round(ModelLR_nb.intercept_, 3))
# 切片ありモデルの誤差評価
Y_pred_01 = ModelLR.predict(X)
Mse_01 = mean_squared_error(Y, Y_pred_01)
Rmse_01 = np.sqrt(Mse_01) # 平方根を算出
print('切片ありモデル RMSE = ', np.round(Rmse_01, 3))
# 切片なしモデルの誤差評価
Y_pred_02 = ModelLR_nb.predict(X)
Mse_02 = mean_squared_error(Y, Y_pred_02)
Rmse_02 = np.sqrt(Mse_02) # 平方根を算出
print('切片なしモデル RMSE = ', np.round(Rmse_02, 3))
RMSEが小さい方、「切片ありモデル」の方が高性能と言えます。
問題4.¶
性能の良いモデルで、製品5の推定をします。
In [10]:
# 製品5の予測
x1 = dat[5, 1] # 製品3の値段
x2 = dat[5, 3] # 製品3のネガ口コミ
x3 = dat[5, 4] # 製品3のポジ口コミ
pred_y = ModelLR.coef_[0] * x1 + ModelLR.coef_[1] * x2 + ModelLR.coef_[2] * x3 + ModelLR.intercept_
print("予測値: ", pred_y)
print("実測値: ", dat[5, 0])
print("差", pred_y-dat[5, 0])
余裕があったら、小数点以下を3桁で丸めましょう。
まずは、製品5を予測できるだけでokです。でも、最終的には、製品5だけと言わず、全部推定できるようになりましょう。以下、for文を利用してすべてを推定するコードです。
In [11]:
e=[]
for i in range(len(dat[:, 0])):
x1 = dat[i, 1] # 値段
x2 = dat[i, 3] # ネガ口コミ
x3 = dat[i, 4] # ポジ口コミ
pred_y = ModelLR.coef_[0] * x1 + ModelLR.coef_[1] * x2 + ModelLR.coef_[2] * x3 + ModelLR.intercept_
print("\n *** 製品 ", i, "***")
print("予測値: ", np.round(pred_y, 2))
print("実測値: ", dat[i, 0])
print("差", np.round(pred_y-dat[i, 0], 2))
e.append( (pred_y - dat[i, 0])**2)
rmse = np.sqrt(np.sum(e)/len(dat[:,0]))
print("\n 実際に計算してみたRMSE: ", np.round(rmse, 3))
まあ、当たったり、外れたりですね。でも、 売れるものは売れる、売れないものは売れないと、それくらいは予測できているようです。
x1, x2, x3に好きな数字を入れれば、好きな予測ができます。
In [12]:
# 好きな数字で予測
x1 = 1500 # 値段は1500円
x2 = 0 # ネガ口コミはまったくない！
x3 = 25 # ポジ口コミはたくさんある！
pred_y = ModelLR.coef_[0] * x1 + ModelLR.coef_[1] * x2 + ModelLR.coef_[2] * x3 + ModelLR.intercept_
print("予測値: ", np.round(pred_y, 3))
200個も売れることになりました。実際、集めたデータの中では、もっとも多いくらいに売れている数ですね。
問題5.¶
In [13]:
# まずは、ネガティブな口コミ数の平均値を求めます。
M = np.mean(dat[:, 3])
# 以下、2つのリストを定義します。
# H_list: 高批判製品群の販売数を格納するリスト
# L_list: 低批判製品群の販売数を格納するリスト
H_list = []
L_list = []
# ネガティブな口コミ数がM以上の販売数をH_listに、
# ネガティブな口コミ数がM未満の販売数をL_listに、
# それぞれ代入します。
for i in range(len(dat[:, 0])):
if dat[i, 3] >= M:
H_list.append(dat[i, 0])
elif dat[i, 3] < M:
L_list.append(dat[i, 0])
# H_listとL_listに格納した販売数の平均値を求め、HとLに代入します。
H = np.mean(H_list)
L = np.mean(L_list)
print("高批判製品群の販売数の平均値: ", H, " 個")
print("低批判製品群の販売数の平均値: ", L, " 個")
print("2群間の差分: ", L-H)
# 2群間の販売数をt検定で比較します。
# 「批判が少ない製品が、批判が多い製品よりも多く売れる」
# という、単に差があるだけでは無く、方向を持つ差の検定なので、
# 片側t検定を採用します。
# また、問題の指示通り、「分散が等しくないt検定」を用います。
res_ttest = st.ttest_ind(H_list, L_list, equal_var = False)
pval = np.round(res_ttest[1], 5)/2 # <- 片側検定なので2で割る
print("t検定のp値: ", pval)
if pval < 0.01:
print("1%水準有意")
elif pval < 0.05:
print("5%水準有意")
elif pval < 0.10:
print("10%水準有意")
else:
print("有意差はありませんでした。")
考察としては、「ネガティブな口コミ数が平均値以上の群（高批判製品群）に属する製品の販売数は60個なのに、そうではない群（低批判製品群）の場合は150個も売れていることがわかりました。この差は90個です。また、その差は1%水準で有意でした。ネガティブな口コミが多いとまずいことがわかります」あたりが理解できればokです。なお、2.5e-05というのは「2.5 x 10の-5乗」、つまり、0.000025です。p値はかなり小さいわけです。
問題6.¶
以下が書ければ、最低限としてokです。
In [14]:
m0 = np.mean(dat[:, 0])
m1 = np.mean(dat[:, 1])
m2 = np.mean(dat[:, 2])
m3 = np.mean(dat[:, 3])
m4 = np.mean(dat[:, 4])
s0 = np.std(dat[:, 0])
s1 = np.std(dat[:, 1])
s2 = np.std(dat[:, 2])
s3 = np.std(dat[:, 3])
s4 = np.std(dat[:, 4])
print("販売数 / 平均値: ", m0, ", 標準偏差: ", s0)
print("値段 / 平均値: ", m1, ", 標準偏差: ", s1)
print("広告費 / 平均値: ", m2, ", 標準偏差: ", s2)
print("N口コミ数 / 平均値: ", m3, ", 標準偏差: ", s3)
print("P口コミ数 / 平均値: ", m4, ", 標準偏差: ", s4)
慣れてきたら、以下のようにかけるようになりましょう。そうしないと、要因数が数十、数百になったとき、大変です。
In [15]:
# 大変良い回答
m, s= [], []
head = ["販売数", "値段", "広告費", "N口コミ数", "P口コミ数"]
for i in range(len(dat[0, :])):
m.append(np.mean(dat[:, i]))
s.append(np.std(dat[:, i]))
print(head[i], "/ 平均値: ", np.round(m[i], 2), ", 標準偏差: ", np.round(s[i], 3))
続いて、約70%と約95%の区間です。集めたデータが正規性を持つと仮定すれば、
- 平均 - 1標準偏差 〜 平均 + 1標準偏差: 約70%（正確には68.27%）のデータが集中
- 平均 - 2標準偏差 〜 平均 + 2標準偏差: 約95%（正確には95.45%）のデータが集中
という傾向を持つことになりますので、値段における上記の区間は、以下のようになります。
In [16]:
print("値段は、", m1 - s1, " から ", m1 + s1, " の区間に約70%存在する。")
print("値段は、", m1 - 2*s1, " から ", m1 + 2*s1, " の区間に約95%存在する。")
もちろん、可能であれば、roundで丸めましょう。
In [17]:
print("値段は、", np.round(m1 - s1, 2), " から ", np.round(m1 + s1, 2), " の区間に約70%存在する。")
print("値段は、", np.round(m1 - 2*s1, 2), " から ", np.round(m1 + 2*s1, 2), " の区間に約95%存在する。")
これは、今回集めた20個の製品についての話ではありません。20個の製品に与えられたイヤホンの値段から、それ以外にある100個、1000個のイヤホンの値段が、だいたいこの範囲に収まることを主張するものです。もちろん、あくまで概算値なので、必ずそうなるわけではありません。詳しく知りたい人は正規分布を勉強しましょう。
おわりに¶
当たり前ですが、コピペでできても意味がないですし、見てなんとなくわかっただけでも、同様に意味がありません。
本当にまっさらな状態で、データだけが与えられたとき、同じように分析できるかを念頭に、問題に取り組んでください。
何もわからないという人は、
- if文やfor文がわからないのか、
- リストやappendがわからないのか、
- numpy配列の扱いがわからないのか、
- 相関・重回帰分析・検定・散布図がわからないのか、
まずはこれを特定し、復習しましょう。