統計学でパン屋のオヤジの不正を見破れ !

こんにちは!
趣味は犬の散歩です!
今回は、シリーズ「コンピューターを使って統計学を学ぶ ?」の第1回です!

Webサービスを提供している弊社では、日々多様なデータを扱っています。
データを元に現状を把握・分析することができ、
更には、中に潜む意味ある情報を顕にすることで、将来に向けて有益な事実・性質を発見することもできます。

まずはデータを「見る」ことから始まりますが、
その道具として、統計という手段を使うことがあります。
ただ、その理屈をきちんと理解することはなかなか難しいと思います。

ここではまず、「案外使えるな、面白いな」と思ってもらえることを目標に、 身近な例を使いつつ、統計の面白さを紹介できたらと思っています。

今回は、統計によって、 パン屋のオヤジがパンの重量を誤魔化していることを暴きたいと思います。

パン屋と私

事例

i社 SK部所属のSさんは、とあるパン屋さんで定期的に「●●パン 1,000グラム」なる商品を購入しています。
そして、Sさんは統計学の勉強のため、毎回重量を計って記録していました。

ある程度記録した後に分析してみたところ、
「ムムムッ ! 何かおかしいで…。平均は950グラムやし、…やし…」と思いました。*1
パン屋のオヤジに問うてみると、「製造機械が故障していた、すぐ直す」とのことのでした。

その後、Sさんは「やれやれ」と思いつつも、とてもおいしいパンなので同じ商品を購入し続けました。
そして暇なSさんは引き続き記録を残し、また分析を行いました。
今回は 平均が 1,000グラムになりました。

しかし、Sさんはパン屋のオヤジが依然として分量をちょろまかしていることを見破り、 再度オヤジに注意しました。
オヤジは驚き、その後不正は無くなりました。

2回目のパンの重量の記録が下になります。
皆さんも下のデータを元に、どこがおかしいのか一度考えてみてください。

Excelで分析できるかと思いますが、以下ではR言語を使って解説をさせて頂きます。

2回目の記録

コンマ区切りで 300個のデータ (小数点以下第4位、単位はグラム) になります。

(平均は、ほぼ 1000になるのですが…)

942.6196,1009.8084,1042.9022,1009.2511,1028.6365,969.0221,964.1312,1061.9494,1007.1534, 1015.3306,1045.6238,1077.3554,1035.5813,1016.0476,991.7477,994.9836,1023.6000,971.2351, 1017.9573,988.2233,
1023.0658,961.0094,1018.1547,1047.1131,1095.7430,974.5152,974.3017, 1010.5197,1013.6388,1011.2781,1021.6128,964.3344,969.2688,1050.7094,1004.0404,1013.0937 ,980.8748,919.9666,978.8745,1066.9316,
1012.0491,1034.5944,964.7270,1031.7082,1024.8805 ,948.5610,1010.3529,988.3572,1028.0611,1016.4201,1034.5434,980.9697,977.4229,986.6875, 1024.0293,969.7522,1025.5132,987.4987,986.9957,1030.4608,
987.0861,941.2553,999.9153 ,969.5143,939.6738,1006.6746,1052.8287,982.1435,1020.0974,1039.4144,1012.4396,1008.6194 ,993.6562,996.6083,1035.5999,946.7522,978.1967,960.5838,1005.8291,993.0979,
1023.6000, 1000.2802,1023.5553,1012.2793,1016.1940,971.6040,1013.7196,1022.1846,985.8477,1008.6194 ,962.4956,969.4809,1021.7654,993.6392,1021.1901,1041.2004,1030.5544,1005.8150,992.4146, 1026.5620,
1027.8279,979.6216,957.6936,956.6530,990.5143,998.9499,977.4040,945.7144 ,939.1462,1010.5328,980.7811,961.4427,1029.8466,1011.5566,985.9380,981.1647,957.9942, 1033.8574,988.8894,1030.7943,
983.6360,1057.8485,1042.8581,989.1975,978.9598,1015.3505, 1007.5134,934.5193,1012.8271,994.7286,995.4682,936.4363,994.7452,1005.0297,982.2436, 1001.0642,1012.4901,1010.8445,999.1350,1026.2960,
958.7274,1018.1503,1029.0578,1034.8865 ,938.6173,1023.4427,1033.3289,1038.6833,1051.8372,998.7876,1016.6356,1013.8618,959.5765, 1006.5952,1048.3189,938.8169,1029.3245,956.7513,1013.7648,984.9830,
1000.9276,967.2816 ,979.6819,970.3903,956.9860,1002.2706,1059.9072,990.0690,997.3581,1049.2533,1018.6856, 1061.3912,1009.2963,1091.6018,1034.3929,994.7286,978.9090,1037.4274,1081.7007,1000.9361,
1033.3289,972.2284,1040.5924,1007.6892,1019.7736,1008.6262,1059.8326,974.8430,1012.0339, 1024.0473,937.4193,996.6093,990.7152,1029.1215,1026.2960,1027.8279,1050.9055,976.6762, 1004.6035,1027.8009,
1004.2662,986.6875,1002.6342,1057.9404,964.0399,1101.6676,1014.5224, 1064.8283,988.7543,976.5099,997.8609,1019.8313,973.7088,971.2765,1062.1122,963.8089, 1024.7650,987.6943,1013.0937,984.5705,
964.9949,966.1785,997.8238,1014.5462,1012.5004, 1010.8434,996.5197,977.6954,991.6123,1005.8179,943.2477,988.3649,1009.3081,964.7396, 1021.8199,1002.0320,964.7757,988.2771,1059.7388,1003.7594,
1048.2966,1023.0502,1025.3499 ,975.0784,1011.2385,1043.8224,1043.1437,976.7518,985.7092,965.2652,1003.1360,1012.5692 ,995.8157,1007.0538,1011.8803,969.8088,1086.9752,963.0390,966.5357,1019.8943,
1014.3999, 1013.3269,934.2373,894.3765,1005.8179,947.7668,1022.9102,994.3349,1052.3795,1028.4340 ,982.0136,968.0434,1024.0293,946.1655,967.5313,1005.4424,949.6137,959.5577,1047.2260, 1072.1457,
1029.0578,1019.6781,985.3175,952.5064,1008.9496,1104.7787,1076.5669,984.6077, 1023.0815,965.4801,1002.2706,986.2137,928.7853,987.1063,1022.2172,941.5854,997.5557 ,997.3581,971.9854,983.1708

解説

正規分布

データを集計する際、ある項目(例えば、重量)について、
ある値(例えば、980グラム)をとる個数(例えば、30個)をグラフにすることが良くあります。
こうした集計結果(分布)について、何かしらの特徴に対し、名前が付けられていることがあります。

今回は、その中で「正規分布」と呼ばれる分布がカギになります。

今回の場合のように、
ある製品をある目標値に向かって製造するなら、誤差の分布たる正規分布に従うと仮定することができます。*2

この分布はグラフで見たとき、左右対称の山形をしています。
そしてその形は、その分布を構成する集団の平均値と標準偏差の値の2つのみによって一意に決まります。

今回は便宜上、本来の正しい製造方法によれば、平均値1,000、標準偏差50の正規分布に従うものとして下さい。

平均を1,000グラム、標準偏差は50グラムとしておくと、下のような形状になります。
f:id:i-plug-develop:20170219194133p:plain
即ち、1,000個製造した場合、重量が1,000グラム辺りのパンは150個程度出来上がり、900グラム辺りのパンも15個程度は存在しそうです。

上の300個のパンの重量データについて分布をグラフにしてみましょう。
(Rによるコード例はこちら)

f:id:i-plug-develop:20170219220853p:plain

んー…
上のヒストグラムに、平均1,000グラム、標準偏差50グラムの10,000個のパンから、300個を無作為抽出した場合のヒストグラムを重ねてみます。
(Rによるコード例はこちら)

f:id:i-plug-develop:20170219221457p:plain

形が少し違います。
水色のヒストグラムの1000グラム付近の個数が、ピンク色に比べて多くなっているし、 左右対称であるはずの分布にしては少し右に偏っています(これはちょっとわかりにくいですが…)

また、今回の正しい製造方法では標準偏差は50辺りになるものとしましたが、
水色の、パン屋から購入したデータでは、その標準偏差は 34.1程度です。
(Rによるコード例はこちら)

これは、50と大きく異なり、 標準偏差(データ集団のバラツキの指標)が小さいということで、 作為的な何かが働いたことを思わせます。

シミュレーション

そうです、パン屋のオヤジは、うるさい客であるSさんには、作ったパンの中から適当に大きいものを渡していたのです。

水色のデータ(パン屋から購入した300個)の値は、
実は、上のオヤジの操作を下のようにシミュレートして作成しました
(Rによるコード例はこちら)

平均950グラム、標準偏差50グラムの正規分布に従う集団の中からn個をランダムに抽出し、 そのn個のうちの最大値を取る、という操作を300回繰り返す

このnの値を1,2,3,…と順番に試し、抽出した300個の平均値が1,000くらいになるnを調べると、
nの値が増えるにつれて、平均値も増加し、
n = 4 のときにだいたい1,000くらいになりました。

このように、最大のものを作為的に抽出した場合、その300個のヒストグラムは左右対称ではなく、右に偏った形になります。
nの値が大きくなると、ある程度まで、平均値とともにこの偏りも大きくなり、一目で正規分布ではないことがわかります。
(今回はこの偏りが少しわかりにくかったですね…)

最後に

まとめ

今回の事例は、フランスの数学者 H.ポアンカレに関する逸話から取り上げました。
この話は有名で、ネット検索しても何件か引っかかります。
ただ、具体的な考察をしているものが無いなーと思っていたので、シミュレーションを試みました。
正規分布標準偏差などの説明をもっとしたかったのですが、興味を持たれた方は是非参考図書などを紐解いてみて下さい !

R言語によるコード例

今回、統計のお話に多くの紙面を使ってしまったので、
R言語とは何か?、その特徴、などはまた別の機会にできたらと思います。
(このページの一番下に簡単なインストール手順を載せました)

● 2回目の購入記録からその平均値と標準偏差を求める。

# 300個のデータをベクトル(Rのデータ構造の一つで、値の集まりを表現するもの)
bread_300 <- c(
942.6196,1009.8084,1042.9022,1009.2511,1028.6365,969.0221,964.1312,1061.9494,1007.1534, 1015.3306,1045.6238,1077.3554,1035.5813,1016.0476,991.7477,994.9836,1023.6000,971.2351, 1017.9573,988.2233, 1023.0658,961.0094,1018.1547,1047.1131,1095.7430,974.5152,974.3017, 1010.5197,1013.6388,1011.2781,1021.6128,964.3344,969.2688,1050.7094,1004.0404,1013.0937 ,980.8748,919.9666,978.8745,1066.9316, 1012.0491,1034.5944,964.7270,1031.7082,1024.8805 ,948.5610,1010.3529,988.3572,1028.0611,1016.4201,1034.5434,980.9697,977.4229,986.6875, 1024.0293,969.7522,1025.5132,987.4987,986.9957,1030.4608, 987.0861,941.2553,999.9153 ,969.5143,939.6738,1006.6746,1052.8287,982.1435,1020.0974,1039.4144,1012.4396,1008.6194 ,993.6562,996.6083,1035.5999,946.7522,978.1967,960.5838,1005.8291,993.0979, 1023.6000, 1000.2802,1023.5553,1012.2793,1016.1940,971.6040,1013.7196,1022.1846,985.8477,1008.6194 ,962.4956,969.4809,1021.7654,993.6392,1021.1901,1041.2004,1030.5544,1005.8150,992.4146, 1026.5620, 1027.8279,979.6216,957.6936,956.6530,990.5143,998.9499,977.4040,945.7144 ,939.1462,1010.5328,980.7811,961.4427,1029.8466,1011.5566,985.9380,981.1647,957.9942, 1033.8574,988.8894,1030.7943, 983.6360,1057.8485,1042.8581,989.1975,978.9598,1015.3505, 1007.5134,934.5193,1012.8271,994.7286,995.4682,936.4363,994.7452,1005.0297,982.2436, 1001.0642,1012.4901,1010.8445,999.1350,1026.2960, 958.7274,1018.1503,1029.0578,1034.8865 ,938.6173,1023.4427,1033.3289,1038.6833,1051.8372,998.7876,1016.6356,1013.8618,959.5765, 1006.5952,1048.3189,938.8169,1029.3245,956.7513,1013.7648,984.9830, 1000.9276,967.2816 ,979.6819,970.3903,956.9860,1002.2706,1059.9072,990.0690,997.3581,1049.2533,1018.6856, 1061.3912,1009.2963,1091.6018,1034.3929,994.7286,978.9090,1037.4274,1081.7007,1000.9361, 1033.3289,972.2284,1040.5924,1007.6892,1019.7736,1008.6262,1059.8326,974.8430,1012.0339, 1024.0473,937.4193,996.6093,990.7152,1029.1215,1026.2960,1027.8279,1050.9055,976.6762, 1004.6035,1027.8009, 1004.2662,986.6875,1002.6342,1057.9404,964.0399,1101.6676,1014.5224, 1064.8283,988.7543,976.5099,997.8609,1019.8313,973.7088,971.2765,1062.1122,963.8089, 1024.7650,987.6943,1013.0937,984.5705, 964.9949,966.1785,997.8238,1014.5462,1012.5004, 1010.8434,996.5197,977.6954,991.6123,1005.8179,943.2477,988.3649,1009.3081,964.7396, 1021.8199,1002.0320,964.7757,988.2771,1059.7388,1003.7594, 1048.2966,1023.0502,1025.3499 ,975.0784,1011.2385,1043.8224,1043.1437,976.7518,985.7092,965.2652,1003.1360,1012.5692 ,995.8157,1007.0538,1011.8803,969.8088,1086.9752,963.0390,966.5357,1019.8943, 1014.3999, 1013.3269,934.2373,894.3765,1005.8179,947.7668,1022.9102,994.3349,1052.3795,1028.4340 ,982.0136,968.0434,1024.0293,946.1655,967.5313,1005.4424,949.6137,959.5577,1047.2260, 1072.1457, 1029.0578,1019.6781,985.3175,952.5064,1008.9496,1104.7787,1076.5669,984.6077, 1023.0815,965.4801,1002.2706,986.2137,928.7853,987.1063,1022.2172,941.5854,997.5557 ,997.3581,971.9854,983.1708 )

# 平均値を求める。
mean(bread_300)    # 1002.307
# 標準偏差を求める。
sd(bread_300)      # 34.14505

● 2回目の購入記録からヒストグラムを作成する (上のコードを実行してからこちらを実行して下さい)。

hist( bread_300, breaks = seq(800, 1200, 10), 
        xlab = "重量(グラム)",
        ylab = "個数",
        ylim = c(0, 80),
        col = "#00ffff40" )

● 平均値1,000、標準偏差50の正規分布に従う集団から300個を抽出したデータのヒストグラムを重ねる (上の2つのコードを実行してからこちらを実行して下さい)

sample_n <- c( 994.1369,987.0971,1047.5415,1002.8592,1054.1306,1015.3721,988.1536,970.0364,912.6958, 966.0189,1013.2729,1035.1735,924.7157,1050.6771,1003.3077,1007.2228,1069.9694,1032.3997, 975.1960,926.9996,1064.8699,941.4907,1049.3500,991.1917,961.6760,999.0000,1061.5649, 929.0185,1006.2891,962.1312,1013.9594,993.7704,942.8584,1029.7913,991.1096,1032.9023, 915.4206,1009.3299,975.8709,1045.4671,947.8475,947.1513,1097.7692,1005.0567,994.5444, 1042.1183,858.8239,893.8004,982.8945,986.1262,995.9118,1041.5223,936.5344,952.7476, 937.0836,965.2329,984.6119,1054.4379,923.8032,983.4004,945.3204,969.5557,895.8705, 972.8970,1128.0027,959.3284,929.1028,945.3997,1117.6408,1007.8059,1085.6142,1070.4184, 895.5513,979.4731,1066.8736,966.2401,1016.6104,984.8446,1052.6221,945.9118,904.4617, 929.0646,980.5280,933.8230,1015.3361,987.4836,996.1149,1088.0353,1020.9122,997.2176, 1040.6960,943.6930,1106.0469,970.2570,985.7500,1007.8669,1052.8253,1002.1506,1033.2751, 980.9386,1017.8374,942.3575,1081.6865,974.1010,987.0210,929.4259,959.6827,950.5818, 1014.5022,949.9070,921.3852,1028.3031,1052.4466,1029.6142,997.9244,884.4263,1029.7378, 1003.8702,1037.6065,935.8152,972.5690,948.4661,1034.8465,981.5814,984.8942,972.2062, 1047.9757,1029.7674,1052.0063,1016.0905,1060.6002,1009.9170,974.7477,961.4822,939.9263, 985.1878,1007.9722,1032.6484,924.1573,1013.9437,971.5437,1012.7835,1010.9390,983.4363, 967.7916,1031.2695,952.8575,925.3614,1046.8850,1074.3239,1090.2499,1008.4107,1020.5373, 921.3591,991.3792,1003.1841,1023.3568,1083.7206,1008.9094,1059.4036,995.5736,973.0823, 1034.1843,939.6046,908.3522,995.4870,1034.0152,987.2099,953.9585,996.4512,1021.1059, 937.8940,1005.3128,991.6352,1004.8901,1049.5250,945.8929,867.6495,930.1818,958.9819, 1009.8446,1078.9360,926.5778,1002.4438,981.5163,976.6745,989.9523,975.8811,852.5694, 1022.4677,930.2307,1046.7965,986.1886,1042.5594,959.0563,969.2619,958.6423,952.9141, 937.4160,912.2349,1049.6675,1017.1742,1005.1197,974.8972,968.2147,1035.2391,1019.6467, 1001.2841,959.7172,984.6958,1036.2826,1037.7327,1008.5656,1020.4331,1031.3752,1027.6935, 1029.8042,987.3990,1083.1286,1065.5840,962.5940,968.6915,966.8113,967.3567,988.4064, 1005.8973,954.8120,984.7683,1012.1604,1008.2556,1017.7931,1007.5473,991.6473,971.2482, 972.6246,1034.5663,956.4438,939.9172,957.0331,1003.9330,998.5091,1064.3668,1016.1435, 931.6898,981.2179,1031.6605,957.8420,1010.5162,1056.3974,1042.2632,1098.9695,1147.1012, 1056.5529,1029.3007,978.3719,952.8251,1086.3560,986.9735,992.6212,997.2452,974.3695, 1102.3596,984.8489,974.1963,1008.1070,1047.4390,1025.6372,907.9092,959.1172,1040.0858, 940.7205,1040.2311,962.6010,960.2923,923.8196,1040.9454,971.9379,1089.2588,1010.5561, 956.8893,966.9814,943.1299,1122.0712,978.1521,1073.7456,1006.8336,948.0014,1056.6535, 958.7573,1012.4727,951.4484,999.2419,967.4731,976.1112,983.4778,992.8126,972.0378, 1028.9414,1031.4803,1015.0866 )

# 正しく製造した場合に、無作為に300個抽出した場合のヒストグラムを上のヒストグラムに重ねる。
sample_n <- sample(rnorm(10000, mean = 1000, sd = 50), 300, replace = F)
hist( sample_n,
        breaks = seq(800, 1200, 10),
        xlab = "重量[g]",
        ylab = "個数",
        ylim = c(0, 80),
        col = "#ff00ff40",
        add = T )             # 重ねて表示させる。

● オヤジの不正のシミュレーション

# 平均値950、標準偏差50の正規分布から10000個の標本を抽出する。
samples_950 <- rnorm(10000, mean = 950, sd = 50)

# (平均950gの)サンプル(全10,000個)の中から i個をランダムに抽出し、
# そのi個のうちの最大値の平均を眺めるコード (iは 1〜10)
res_list <- list()
for (i in 1:10) {
    s <- c()
    for (j in 1:300) {
        s <- append( s, max(sample(samples_950, i, replace = F)) )
    }
    l = list(s)
    res_list <- append(res_list, l)
}

for (i in 1:10) {
  print( mean(res_list[[i]]) )
}

参考図書

R、RStudioのインストール

RStudioには、デスクトップ版とサーバー版がありますが、今回は手元のパーソナルコンピューターにインストールして使用したいと思いますので、デスクトップ版をインストールします。 http://www.rstudio.com/ide/download/desktop から
Mac版は、「RStudio ●.●.● - Mac OS X 10.6+ (64-bit) 」を、
Windows版は、「RStudio ●.●.● - Windows Vista/7/8/10 」をダウンロードし、
そのファイルをダブルクリックしてインストール作業をして下さい。

f:id:i-plug-develop:20170215145819p:plain

*1:この「おかしい」の根拠も統計学の良い使用例になるのですが、ここでは扱いません

*2:正規分布は、数学者C.F.ガウスの測定誤差の研究の中で発見されました。多くの現象に当てはまり、理論上も応用上も非常に重要な確率分布です