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])
販売数と値段の関係数:  -0.44772120028952445  p値:  0.04775386187686266
販売数と広告費の相関係数:  -0.17001265285902242  p値:  0.47361711529507955
販売数とN口コミの相関係数:  -0.7452767272271268  p値:  0.0001626733205946405
販売数とP口コミの相関係数:  0.7398212794024966  p値:  0.0001925770784567566

小数点以下の細かい情報には、正直そんなに意味がないので、数桁で丸めます。ここまでできるとちょっとだけ加点。

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)
販売数と値段の相関係数:  -0.448  p値:  0.04775
販売数と広告費の相関係数:  -0.17  p値:  0.47362
販売数とN口コミの相関係数:  -0.745  p値:  0.00016
販売数とP口コミの相関係数:  0.74  p値:  0.00019

分析に慣れている人は、以下のように書きます。何も考えずにここまでかけるとベストです。ポイントは、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  p値:  0.04775   : 5%水準有意
販売数と 広告費 の相関係数:  -0.17  p値:  0.47362   n.s.
販売数と N口コミ数 の相関係数:  -0.745  p値:  0.00016   : 1%水準有意
販売数と P口コミ数 の相関係数:  0.74  p値:  0.00019   : 1%水準有意

結果・考察としては、販売数と値段の相関係数は-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))
推定モデル(切片あり): y^ =  -0.007  x1 +  -2.691  x2 +  3.751  x3 +  123.673
推定モデル(切片なし): y^ =  -0.001  x1 +  -0.006  x2 +  7.906  x3 +  0.0
切片ありモデル RMSE =  32.102
切片なしモデル RMSE =  39.276

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])
予測値:  75.39181758565994
実測値:  60.0
差 15.391817585659936

余裕があったら、小数点以下を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))
 *** 製品  0 ***
予測値:  30.25
実測値:  10.0
差 20.25

 *** 製品  1 ***
予測値:  39.52
実測値:  20.0
差 19.52

 *** 製品  2 ***
予測値:  69.29
実測値:  30.0
差 39.29

 *** 製品  3 ***
予測値:  70.94
実測値:  40.0
差 30.94

 *** 製品  4 ***
予測値:  13.39
実測値:  50.0
差 -36.61

 *** 製品  5 ***
予測値:  75.39
実測値:  60.0
差 15.39

 *** 製品  6 ***
予測値:  99.2
実測値:  70.0
差 29.2

 *** 製品  7 ***
予測値:  104.91
実測値:  80.0
差 24.91

 *** 製品  8 ***
予測値:  102.35
実測値:  90.0
差 12.35

 *** 製品  9 ***
予測値:  87.6
実測値:  100.0
差 -12.4

 *** 製品  10 ***
予測値:  144.84
実測値:  110.0
差 34.84

 *** 製品  11 ***
予測値:  144.61
実測値:  120.0
差 24.61

 *** 製品  12 ***
予測値:  65.58
実測値:  130.0
差 -64.42

 *** 製品  13 ***
予測値:  132.44
実測値:  140.0
差 -7.56

 *** 製品  14 ***
予測値:  144.18
実測値:  150.0
差 -5.82

 *** 製品  15 ***
予測値:  168.86
実測値:  160.0
差 8.86

 *** 製品  16 ***
予測値:  160.64
実測値:  170.0
差 -9.36

 *** 製品  17 ***
予測値:  116.94
実測値:  180.0
差 -63.06

 *** 製品  18 ***
予測値:  128.7
実測値:  190.0
差 -61.3

 *** 製品  19 ***
予測値:  200.35
実測値:  200.0
差 0.35

 実際に計算してみたRMSE:  32.102

まあ、当たったり、外れたりですね。でも、 売れるものは売れる、売れないものは売れないと、それくらいは予測できているようです。

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))
予測値:  207.09

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.0  個
低批判製品群の販売数の平均値:  150.0  個
2群間の差分:  90.0
t検定のp値:  2.5e-05
1%水準有意

考察としては、「ネガティブな口コミ数が平均値以上の群(高批判製品群)に属する製品の販売数は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)
販売数 / 平均値:  105.0 , 標準偏差:  57.66281297335398
値段 / 平均値:  3391.5 , 標準偏差:  1564.1907652201505
広告費 / 平均値:  341.0 , 標準偏差:  197.38034349954913
N口コミ数 / 平均値:  16.85 , 標準偏差:  8.568984770671493
P口コミ数 / 平均値:  13.35 , 標準偏差:  6.405271266699015

慣れてきたら、以下のようにかけるようになりましょう。そうしないと、要因数が数十、数百になったとき、大変です。

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)) 
販売数 / 平均値:  105.0 , 標準偏差:  57.663
値段 / 平均値:  3391.5 , 標準偏差:  1564.191
広告費 / 平均値:  341.0 , 標準偏差:  197.38
N口コミ数 / 平均値:  16.85 , 標準偏差:  8.569
P口コミ数 / 平均値:  13.35 , 標準偏差:  6.405

続いて、約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%存在する。")
値段は、 1827.3092347798495  から  4955.690765220151  の区間に約70%存在する。
値段は、 263.11846955969895  から  6519.8815304403015  の区間に約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%存在する。")
値段は、 1827.31  から  4955.69  の区間に約70%存在する。
値段は、 263.12  から  6519.88  の区間に約95%存在する。

これは、今回集めた20個の製品についての話ではありません。20個の製品に与えられたイヤホンの値段から、それ以外にある100個、1000個のイヤホンの値段が、だいたいこの範囲に収まることを主張するものです。もちろん、あくまで概算値なので、必ずそうなるわけではありません。詳しく知りたい人は正規分布を勉強しましょう。

おわりに

当たり前ですが、コピペでできても意味がないですし、見てなんとなくわかっただけでも、同様に意味がありません。

本当にまっさらな状態で、データだけが与えられたとき、同じように分析できるかを念頭に、問題に取り組んでください。

何もわからないという人は、

  • if文やfor文がわからないのか、
  • リストやappendがわからないのか、
  • numpy配列の扱いがわからないのか、
  • 相関・重回帰分析・検定・散布図がわからないのか、

まずはこれを特定し、復習しましょう。