本節書摘來自華章計算機《數學模組化：基于r》一書中的第1章，第1.7節,作者：薛 毅 更多章節内容可以通路雲栖社群“華章計算機”公衆号檢視。
1.7.1 問題的提出 該題選自2013年“深圳杯”數學模組化夏令營a題.
“民以食為天”,食品安全關系到千家萬戶的生活與健康.随着人們對生活品質的追求和安全意識的提高,食品安全已成為社會關注的熱點,也是政府民生工程的一個主題.城市食品的來源越來越廣泛,人們消費加工好的食品的比例也越來越高,是以除食材的生産收獲外,食品的運輸、加工、包裝、貯存、銷售以及餐飲等每一個環節皆可能影響食品的品質與安全.另一方面,食品品質與安全又是一個專業性很強的問題,其标準的制定和抽樣檢測及評價都需要科學有效的方法.深圳是食品抽檢、監督最統一、最規範、最公開的城市之一.請下載下傳2010年、2011年和2012年深圳市的食品抽檢資料（注意蔬菜、魚類、雞鴨等抽檢資料的擷取） 資料下載下傳網站：www.szaic.gov.cn（深圳市市場監督管理局網站）.并根據這些資料來讨論：
1） 如何評價深圳市這三年各主要食品領域微生物、重金屬、添加劑含量等安全情況的變化趨勢.
2） 從這些資料中能否找出某些規律性的東西：如食品産地與食品品質的關系；食品銷售地點（即抽檢地點）與食品品質的關系；季節因素等.
3） 能否改進食品抽檢的辦法,使之更科學更有效地反映食品品質狀況且不過分增加監管成本（食品抽檢是需要費用的）,例如對于抽檢結果穩定且抽檢頻次過高的食品領域該作怎樣的調整？
解決上述問題的最大困難是資料整理,因為這三年中的資料格式不統一（三年的excel表的格式不同）,檔案不統一（有excel表和word文檔）,有的表隻有不合格食品,沒有說明本次抽檢共抽檢了多少食品.這些都給資料分析帶來困難,需要花大量的時間來分析整理 本節的資料是根據2013年“深圳杯”數學模組化夏令營的幾篇優秀論文整理得到的..
1.7.2 問題1：三年各主要食品領域安全情況的變化趨勢
對深圳市2010年至2012年的所有能夠得到的不合格資料進行彙總,并将檢測名額大緻分為以下四類：微生物、重金屬、食品添加劑和其他類,彙總結果如表1.18所示.
表1.18 2010—2012年各類主要食品抽檢情況彙總表
年度抽檢總數抽檢不合格種類微生物重金屬添加劑其他類總數
2010 2261 60 7 31 40 138
2011 7814 250 264 239 79 832
2012 10607 116 35 29 78 258
食品抽檢隻有兩種結果,一是合格,二是不合格.是以,選擇二項分布的方法進行分析,這裡采用prop.test()函數,pairwise.prop.test()函數和prop.trend.test()函數完成安全情況分析.
用prop.test()函數檢驗這三年抽檢不合格的比率是否相同,即檢驗h0:p1=p2=p3, h1:p1,p2,p3不全相同如果拒絕原假設,則表明這三年的不合格比率有顯著差異.
用pairwise.prop.test()函數作多重比率檢驗,即檢驗誰與誰有顯著差異,該函數的使用格式為pairwise.prop.test(x, n,
p.adjust.method = p.adjust.methods, ...)參數x為整數向量,表示試驗成功的次數,或者為一個2列的矩陣,第1清單示成功的次數,第2清單示失敗的次數.
n為整數向量,表示試驗的次數,當x為矩陣時,該值無效.
p.adjust.method為p值的調整方法(細節見2.2.2節).
...為附加參數.
用prop.trend.test()函數作比率趨勢性檢驗,即檢驗h0:p1=p2=…=pk, h1:p1≤p2≤…≤pk或者檢驗h0:p1=p2=…=pk, h1:p1≥p2≥…≥pk其使用格式為prop.trend.test(x, n, score = seq_along(x))參數x為事件數.n為試驗次數.score為分組得分,預設值為自然順序.
将資料按表格形式存入資料檔案(檔案名:summer_1.data),在表頭中,sampling表示抽檢總數,microorganism表示微生物方面的食品不合格數,metal表示重金屬方面的食品不合格數,additive表示添加劑含量方面的食品不合格數,others表示其他方面的食品不合格數,unqualified表示食品不合格的總數.用read.table()函數讀取表格形式的資料檔案.
(1) 關于微生物方面的食品安全情況的變化趨勢> rt <- read.table("summer_1.data")
prop.test(x = rt$microorganism, n = rt$sampling)
3-sample test for equality of proportions without
continuity correction
data: rt$microorganism out of rt$sampling
x-squared = 103.3388, df = 2, p-value < 2.2e-16
alternative hypothesis: two.sided
sample estimates:
prop 1prop 2prop 3
0.026536930.031993860.01093617拒絕原假設,說明關于微生物方面的食品安全情況這三年是有差異的.> pairwise.prop.test(x = rt$microorganism, n = rt$sampling)
pairwise comparisons using pairwise comparison of
proportions
12
20.21-
32.4e-08< 2e-16
p value adjustment method: holm 從計算結果來看,2010年與第2012年,2011年與2012年微生物方面的食品安全情況有顯著差異.> prop.trend.test(x = rt$microorganism, n = rt$sampling)
chi-squared test for trend in proportions
data: rt$microorganism out of rt$sampling,
using scores: 1 2 3
x-squared = 70.1004, df = 1, p-value < 2.2e-16從趨勢檢驗來看,拒絕原假設,說明微生物方面的食品不合格率在逐年減少.
(2) 關于重金屬方面的食品安全情況的變化趨勢> prop.test(x = rt$metal, n = rt$sampling)
data: rt$metal out of rt$sampling
x-squared = 310.7125, df = 2, p-value < 2.2e-16
0.0030959750.0337855130.003299708 拒絕原假設,說明關于重金屬方面的食品安全情況這三年是有差異的.> pairwise.prop.test(x = rt$metal, n = rt$sampling)
27.1e-15 -
31< 2e-16
p value adjustment method: holm 從計算結果來看,2010年與第2011年,2011年與2012年重金屬方面的食品安全情況有顯著差異.> prop.trend.test(x = rt$metal, n = rt$sampling)
chi-squared test for trend in proportions
data: rt$metal out of rt$sampling,
using scores: 1 2 3
x-squared = 65.8371, df = 1, p-value = 4.898e-16從趨勢檢驗來看,拒絕原假設,說明重金屬方面的食品不合格率在逐年減少.
(3) 關于添加劑含量方面的食品安全情況的變化趨勢> prop.test(x = rt$additive, n = rt$sampling)
data: rt$additive out of rt$sampling
x-squared = 245.0698, df = 2, p-value < 2.2e-16
0.0137107470.0305861270.002734044拒絕原假設,說明添加劑含量方面的食品安全情況這三年是有差異的.> pairwise.prop.test(x = rt$additive, n = rt$sampling)
21.7e-05-
32.3e-11< 2e-16
p value adjustment method: holm 從計算結果來看,2010年與第2011年,2010年與第2012年,2011年與2012年添加劑含量方面的食品安全情況有顯著差異.> prop.trend.test(x = rt$additive, n = rt$sampling)
data: rt$additive out of rt$sampling,
using scores: 1 2 3
x-squared = 111.1509, df = 1, p-value < 2.2e-16從趨勢檢驗來看,拒絕原假設,說明添加劑含量方面的食品不合格率在逐年減少.
1.7.3 問題2：找出某些規律性的東西
1.食品生産到銷售各環節與食品品質的關系
食品在生産到銷售的過程中,大緻可分為:餐飲、儲存、生産、使用、運輸和銷售等環節,将2010～2012年食品在各環節的抽檢情況彙總(見表1.19).
2.食品産地與食品品質的關系
不同食品産地對食品品質的影響主要展現在三個方面：運輸距離的長短,不同食品産地本身的品質和深圳市對不同産地産品流入市場時的監管力度.
在資料進行中,由于部分産地的檢測數目較小,該省份或城市的合格率會産生較大的偶然性,是以統一剔除檢測數目小于50的産地.篩選後所得的36個食品産地與合格率資訊,如表1.20所示.
考查食品産地與食品品質之間的關系,可從近距離(800km以内)和遠距離(800km以上)兩方面考慮.考查的名額是食品産地到銷售地的距離和食品抽檢的不合格率.由于這兩項名額的資料不一定滿足正态分布,是以這裡采用秩相關檢驗——spearman相關檢驗.
将資料按表格形式存入資料檔案(檔案名:summer_22.data),在表頭中,distance表示食品産地距深圳市的距離,sampling表示抽檢總數,unqualified表示食品不合格數.以下是r程式(檔案名:summer_22.r)及計算結果.
先計算近距離(800km以内)的食品安全情況：rt <- read.table("summer_22.data")
rt$ratio <- rt$unqualified / rt$sampling
cor.test(~ distance + ratio, data = rt, subset = 1:18,
method = "spearman")
spearman's rank correlation rho
data: distance and ratio
s = 804.9153, p-value = 0.5018
alternative hypothesis: true rho is not equal to 0
rho
0.169334接受原假設,說明食品産地與食品品質是獨立的.
再計算遠距離(800km以上)的食品安全情況：cor.test(~ distance + ratio, data = rt, subset = 19:36,
s = 1256.987, p-value = 0.231
-0.2972007接受原假設,說明食品産地與食品品質是獨立的.
将近距離和遠距離的食品放在一起讨論：cor.test(~ distance + ratio, data = rt,
s = 11417.17, p-value = 0.003876
-0.4693908拒絕原假設,說明食品産地與食品品質是有關的.從相關系數來看,相關系數的符号為“-”,也就是說,食品産地越遠,食品抽檢的不合格率越小,食品品質越高.這一結論與我們的生活常識似乎是沖突的,應該是産地越遠,食品抽檢的不合格率越高.
仔細分析,上述計算得到的結論是合理的,因為産地越遠,食品的包裝就會越好,運輸條件也相對近距離要好,否則到了銷售地後,食品的包裝損壞嚴重,會給生産廠家帶來不必要的損失.
3.食品銷售地點（即抽檢地點）與食品品質的關系
為了研究這個問題,将深圳市按照城區劃分,各個城區的合格數與不合格數如表1.21所示.
表1.21 2010—2012年各地區食品安全情況彙總表
抽檢總數 5773 2682 690 4899 2754 2893 734 986
不合格數 293 128 67 256 113 107 32 53
作列聯表檢驗,檢驗食品品質與抽檢地是否獨立.以下是程式(檔案名:summer_23.r)和計算結果：x <- scan()
5773 2682 690 4899 2754 2893 734 986
293 128 67 256 113 107 32 53
x <- matrix(c(x［9:16］, x［1:8］-x［9:16］), nr = 2, byrow = t)
chisq.test(x)
pearson's chi-squared test
data: x
x-squared = 49.5067, df = 7, p-value = 1.805e-08拒絕原假設,說明食品的抽檢地點與食品品質是有關系的.這個結論或許與每個地區的食品來源有關.
也可以作比率檢驗(prop.test()函數),其結果是相同的,這一檢驗留給讀者完成.
4.季節與食品品質的關系
為了研究季節與食品品質的關系,整理出飲品、果蔬、糧食、焙烤、肉蛋和水産品共6類主要食品的抽檢情況,如表1.22所示.
表1.22 六類主要食品分季度彙總表
飲品 777 15 254 7 1225 29 988 25
果蔬 861 33 275 29 518 25 1049 50
糧食 809 47 721 32 1171 16 1156 36
焙烤 476 130 205 17 989 18 491 31
肉蛋 426 62 391 36 280 16 766 76
水産 98 2 222 10 373 22 330 18
将資料以表格形式存放在文本檔案(summer_24.data)中,使用列聯表作獨立性檢驗.
(1) 關于飲品類食品的檢驗情況rt <- read.table("summer_24.data")
x <- matrix(as.numeric(rt［1, ］), nc = 4)
pearson's chi-squared test
x-squared = 0.8809, df = 3, p-value = 0.83程式中的as.numeric()函數是将其他類型資料強制轉換成數值型.計算結果是,接受原假設,即飲品類食品的安全情況與季節無關.
(2) 關于果蔬類食品的檢驗情況x <- matrix(as.numeric(rt［2, ］), nc = 4)
x-squared = 17.4588, df = 3, p-value = 0.0005687拒絕原假設,即果蔬類食品的安全情況與季節有關.
(3) 關于糧食類食品的檢驗情況x <- matrix(as.numeric(rt［3, ］), nc = 4)
x-squared = 29.5963, df = 3, p-value = 1.678e-06拒絕原假設,即糧食類食品的安全情況與季節有關.
(4) 關于焙烤類食品的檢驗情況x <- matrix(as.numeric(rt［4, ］), nc = 4)
x-squared = 197.4468, df = 3, p-value < 2.2e-16拒絕原假設,即焙烤類食品的安全情況與季節有關.
(5) 關于肉蛋類食品的檢驗情況x <- matrix(as.numeric(rt［5, ］), nc = 4)
x-squared = 12.5369, df = 3, p-value = 0.005753拒絕原假設,即肉蛋類食品的安全情況與季節有關.
(6)關于水産類食品的檢驗情況x <- matrix(as.numeric(rt［6, ］), nc = 4)
alternative hypothesis: two.sided接受原假設,即水産類食品的安全情況與季節無關.注:由于第一季度的不合格數為2,小于5,是以,使用fisher精确檢驗.
1.7.4 問題3：如何改進食品的抽檢辦法
這個問題實際上相當的複雜,改進抽檢辦法的目的就是在保證抽檢品質的前提下,減少抽檢次數,進而降低抽檢成本.這裡涉及抽檢方法、抽檢成本,以及優化等方面的知識.
本小節以簡單随機抽樣為基礎,簡單計算出抽樣的樣本容量.
設y1,y2,…,yn為随機抽樣的樣本,y=∑ni=1yi為樣本均值,n為總體的總數,f=n/n為抽樣比.由抽樣理論可知,var(y)=1－fns2=1n－1ns2(1.79)其中s2為總體的方差.由式(1.79)得到1n=1n+var(y)s2(1.80)式(1.80)雖然明确地給出了樣本容量n的計算公式,但式中的總體方差s2是未知的,而抽樣誤差var(y)也與總體方差s2有關.
通常人們用置信度1－α和絕對誤差限y－y≤d來計算抽樣的樣本容量,即考慮p{y－y≤d}=1－α（1.81）其中y為總體的均值.由式(1.81)和雙側分位點uα/2的定義分别得到py－yvar(y)≤dvar(y)=1－α（1.82）
py－yvar(y)≤uα/2=1－α（1.83）對比式(1.82)和式(1.83)得到,uα/2=d/var(y),即var(y)=d2/u2α/2（1.84）将式(1.84)代入式(1.85),得到1n=1n+d2u2α/2s2（1.85）如果考慮置信度和相對誤差py－yy≤r=1－α（1.86）類似前面的推導,可以得到1n=1n+r2y2u2α/2s2（1.87）對于食品的安全檢查,人們最關心的是食品的不合格率p,是以,用p替換y,并令s2=pq,其中q=1－p.代入式(1.87)得到1n=1n+r2pu2α/2q（1.88）在實際應用中,通常取α=0.05,用正态分布的分位點近似uα/2.如果某類食品的總數n=10000,食品的不合格率p=0.05,檢驗的相對誤差r=0.1,即10%,按式(1.88)計算得到n=4220約為總數的一半.
這個抽檢次數遠遠大于目前的抽檢次數,是以,隻使用簡單随機抽檢的方法是不能減少抽檢次數的.關于其他抽檢方法的讨論已超出本書的内容,是以不再進行讨論.
1.7.5 結論
(1) 深圳市2010年至2012年三年來的食品安全狀況逐年變好;
(2) 食品的餐飲、儲存等各環節與食品品質有關;
(3) 食品産地與食品品質有關,而且是負相關;
(4) 食品的抽檢地點與食品品質有關;
(5) 飲品和水産類食品的品質與季節無關,果蔬、糧食、焙烤和肉蛋類食品的品質與季節有關.
習題1
求學生身高的平均值、中位數、方差、極差和四分位極差.
2.假定某地區成年男子的身高(機關:cm)服從正态分布x~n(170,7.692),求該地區成年男子身高超過175cm的機率.如果希望有95%以上的人在通過門時不低頭,門的高度至少應為多少?
3.繪出習題1中學生身高資料的直方圖、核密度估計曲線.(要求:直方圖高度為頻率/間距,并将核密度估計曲線與直方圖畫在一張圖上).
4.繪出習題1中學生身高資料的正态q-q圖,說明該資料是否接近正态分布.
5.正常男子血小闆計數均值為225×109/l,今測得20名男性油漆作業勞工的血小闆計數值(機關:109/l)如下：
220 188 162 230 145 160 238 188 247 113
126 245 164 231 256 183 190 158 224 175
假定資料服從正态分布,問油漆勞工的血小闆計數與正常成年男子有無顯著差異?
6.已知某種燈泡壽命服從正态分布,在某星期所生産的該燈泡中随機抽取10隻,測得其壽命(機關:h)為
1067 919 1196 785 1126 936 918 1156 920 948
求這個星期生産出的燈泡能使用1000h以上的機率.
7.為研究國産四類新藥阿卡波糖膠囊效果,某醫院用40名Ⅱ型糖尿病病人進行同期随機對照實驗.試驗者将這些病人随機等分到試驗組(阿卡波糖膠囊組)和對照組(拜唐蘋膠囊組),分别測得試驗開始前和8周後空腹血糖,計算得到空腹血糖下降值如表1.23所示.假設資料服從正态分布,試用t檢驗(讨論方差相同和方差不同兩種情況)和成對資料的t檢驗來判斷:國産四類新藥阿卡波糖膠囊與拜唐蘋膠囊對空腹血糖的降糖效果是否相同.分析三種檢驗方法哪一種最好,并說明理由.
表1.23 實驗組與對照組空腹血糖下降值
試驗組-0.70-5.602.002.800.703.504.005.807.10-0.50
2.50-1.601.703.000.404.504.602.506.00-1.40
對照組3.706.505.005.200.800.200.603.406.60-1.10
6.003.802.001.602.002.201.203.101.70-2.00
8.一項調查顯示某城市老年人口比重為14.7％.該市老年研究協會為了檢驗該項調查是否可靠,随機抽選了400名居民,發現其中有57位是老年人.問調查結果是否支援該市老年人口比重為14.7％的看法(α=0.05).
9.研究者發現,婦女患乳腺癌可能與初次分娩時的年齡有關,表1.24給出國際衛生組織在1970年的報告.試分析初次分娩婦女各年齡段患乳腺癌的比例是否相同.
表1.24 初次分娩年齡與乳腺癌患病人數
初次分娩年齡<2020~2425~2930~34≥3
5乳腺癌數32012061011463220
調查總數1742563839041555626
10.某項調查詢問了2000名青年人.問題是:“你認為我們的生活環境比過去更好、更差,還是沒有變化?”有800人覺得“越來越好”,有720人感覺“一天不如一天”,有400人表示“沒有變化,一直如此”,還有80人說不知道.試分析:在總體中是否可以認為“我們的生活環境比過去更好”的人比“我們的生活環境比過去更差”的人多?
11.在某養魚塘中,根據過去經驗,魚的長度的中位數為14.6cm,現對魚塘中魚的長度進行一次估測,随機地從魚塘中取出10條魚，長度如下：
13.32 13.06 14.02 11.86 13.58
13.77 13.51 14.42 14.44 15.43
将它們作為一個樣本進行檢驗.試分析,該魚塘中魚的長度是在中位數之上,還是在中位數之下？(1)用符号檢驗;(2)用wilcoxon符号秩檢驗.
12.考慮例1.19若新方法與原方法得到排序結果改為表1.25所示的情形,能否說明新方法比原方法顯著提高了教學效果?
表1.25 學生數學能力排序結果
新方法467910
原方法12358
13.為了比較一種新療法對某種疾病的治療效果,将40名患者随機地分為兩組,每組20人,一組采用新療法,另一組采用原标準療法.經過一段時間的治療後,對每個患者的療效作仔細地評估,并劃分為差、較差、一般、較好和好五個等級.兩組中處于不同等級的患者人數如表1.26所示.試分析,由此結果能否認為新方法的療效顯著地優于原療法(α=0.05)?
表1.26 不同方法治療後的結果
等級差較差一般較好好
新療法組01973
原療法組221141
14.據經驗知,機床發生故障的機率服從均勻分布,某工廠中的房間在一周内統計所有車床發生故障的頻數資料如表1.27所示,試分析,這組資料是否服從均勻分布?
表1.27 某工廠中的房間在一周内所有車床發生故障的頻數資料
星期一二三四五六
頻數78391617
mendel用豌豆的兩對相對性狀進行雜交實驗,黃色圓滑種子與綠色皺縮種子的豌豆雜交後,第二代根據自由組合規律,理論分離比為黃圓∶黃皺∶綠圓∶綠皺=916∶316∶316∶116實際實驗值為:黃圓315粒,黃皺101粒,綠圓108粒,綠皺32粒,共556粒,問此結果是否符合自由組合規律?
16.在一次實驗中,每隔一定的時間觀察一次由某種鈾所放射到達計數器的α粒子數x,共觀察了100次,其資料如表1.28所示.試分析,能否認為α粒子數x服從poisson分布.
表1.28 觀察到的α粒子數
粒子數01234567891011
頻數 1516172611992121
17.一般認為長途電話通過電話總機的過程是一個随機過程,其間打進電話的時間間隔服從指數分布.某個星期下午1:00以後最先打進的10個電話的時間為
1:06 1:08 1:16 1:22 1:23 1:34 1:44 1:47 1:51 1:57
試用kolmogorov-smirnov檢驗分析打進電話的時間間隔是否服從指數分布.
18.試用shapiro-wilk檢驗判斷習題1中的身高資料是否服從正态分布.
19.在高中一年級男生中抽取300名考察其兩個屬性:b是1500米長跑,c是每天平均鍛煉時間,得到4×3列聯表,如表1.29所示.試對α=0.05,檢驗b與c是否獨立.
表1.29 300名高中學生體育鍛煉的考察結果1500米長跑記錄鍛煉時間2小時以上1~2小時1小時以下合計5″01′~5″30′451210675″31′~6″00′462028946″01′~6″30′282330816″31′~7″00′11123558合計13067103300
20.為研究分娩過程中使用胎兒電子監測儀對剖腹産率有無影響,對5824例分娩的經産婦進行回顧性調查,結果如表1.30所示,試進行分析.
表1.30 5824例經産婦回顧性調查結果剖腹産胎兒電子監測儀使用未使用合計是358229587否249227455237合計285029745824
21.為研究人腦的左右半球惡性惡性良性腫瘤的發病率是否有顯著差異,對人腦惡性惡性良性腫瘤和良性惡性良性腫瘤的發病情況作了調查,調查結果如表1.31所示.試進行分析.
表1.31 人腦左右半球惡性惡性良性腫瘤和良性惡性良性腫瘤的發病情況
左半球9312
右半球134
合計10616
22.表1.32列出某高中18名學生某門課程的聯考成績和模拟考試成績,這組資料能否說明:聯考成績與模拟考試成績是相關的？
表1.32 聯考成績和模拟考試成績
學号 聯考成績 模拟成績 學号 聯考成績 模拟成績 學号 聯考成績 模拟成績1 87 90 7 78 65 13 90 100
2 76 98 8 91 90 14 92 97
3 77 92 9 76 84 15 100 97
4 85 87 10 100 92 16 100 95
5 89 87 11 96 100 17 90 94
6 83 62 12 96 98 18 99 100
23.調查某大學學生每周學習時間(機關:h)與得分的平均等級之間關系,現抽查10名學生的資料如表1.33所示.表中等級10表示最好,1表示最差.試用秩相關檢驗(spearman檢驗和kendall檢驗)分析學習時間與學習等級有無關系.
表1.33 每周學習時間與學習等級(機關:h)
學習時間24172041522346181529
學習等級81479510326
![資料挖掘分類技術[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)