風中一匹狼

#屬性資料的分析

#擲兩個骰子，兩個骰子的點數一樣的機率為多少?兩個骰子總共點數不足5的機率為多少
dbinom(2, size = 2, prob = 1/6)*6
first<- (1:6) ;second<- (1:6)
results<- c()
i<-1
while (length(results) < 100000) {
  results[i] <- (sample(first , 1)+ sample(second,1)) <5
  i<-i+1
}
sum(results)/100000

#為了檢定一個虛擬假說，進行四次試驗(互相獨立)。如果虛擬假說是正確的，至少得到一次顯著的結果 (p< 0.05)的機率有多少?
1-dbinom(4,size = 4, prob = 0.95)

#一群乳牛乳脂率低於4%佔60%，如果隨機取10頭牛，其中剛好6頭牛的乳脂率低於4%的機率有多少?
dbinom(6, size = 10, prob = 0.6) #取樣10次，其中6次成功(< 4%)的機率
#至少有4頭牛之乳脂率大於或等於4%的機率 (也就是1-成功三次以內的機率)
1- pbinom(3, size= 10, prob = 0.4) #>=4%的機率為0.4，累積成功3次的機率，100%扣除該機率

#動物藥品聲稱對某疾病有90%治癒率，使用於10頭病畜，治癒8頭。此藥是90%治癒率嗎? 以二項分布與常態近似法測定，結果是否不同? 哪種方法才是正確
#二項分布
1- dbinom(10,10,0.9)- dbinom(9, 10, 0.9) #至少治癒8頭機率為0.26，並沒有足夠證據拒絕治癒率低於90%
#常態近似法
p<- 0.9
pnorm(8, mean= 9, sd= sqrt(10*p*(1-p))/0.5) #標準偏差除以0.5校正連續性。
#常態近似法與二項分布算出的機率有些差異是因為P在接近1 (90%)且樣本數少的情況下，二項分布是有些傾斜不像常態分不一樣對稱。在該情境下二項分布會比較適用。
#http://episte.math.ntu.edu.tw/articles/sm/sm_16_07_1/index.html poisson分布介紹

#同上題，如果此藥使用於100頭病畜，治癒80頭，此藥是否有90%治癒率? 該藥治癒率的95%信賴區間為何?
pnorm(80, mean = 90, sd=sqrt(100*0.9*0.1))
xbar <-0.8 ;n<- 100; sd<- sqrt(0.8*0.2/n)
alpha = .05
zscore <- qnorm(1-alpha/2)
c(xbar-zscore*sd, xbar+zscore*sd)

#新開發果汁與老牌果汁經消費者測試，測試人數總共115人，30人兩種都喜歡，20人兩種都不喜歡，25人較喜歡老牌，40人較喜歡新的。你認為消費是否比較喜歡新的果汁?
1- pchisq((abs(40-25)-1)^2/(40+25), df= 1)
(abs(40-25)-1)^2/(40+25) > qchisq(1-0.05, 1) #新舊喜好沒有顯著差異
(abs(40-25)-1)/sqrt(40+25) > qnorm(1-0.05/2) #同上結果

#某地區乳房炎發生的情形如下，乳牛患乳房炎是否影響女牛患乳房炎的機率?
cow<- matrix(c(T,T,45,
              T,F,25,
              F,T,125,
              F,F,155),nrow= 4 ,ncol= 3,byrow=T, 
              dimnames= list(c(1:4),c("親代","子代","病例數目")))
cow<- as.data.frame(cow)
total<- sum(cow[ ,3])
p<- (45+25)/total; q<- (125+155)/total #p/q代表母牛患乳房炎機率
u<- (45+125)/total; w<- (25+155)/total #u/W代表女牛患乳房炎機率
cow$expected_number <- c(p*u,p*w,q*u,q*w)*total
(45-34)^2*sum(1/ cow$expected_number) > qchisq(1-0.05, 1) #母牛患乳房炎會顯著影響女牛
#或是使用2*2 contingency table 亦可算出Chi-Squared value
total* (45*155-125*25)^2 / ((34+136)*(36+144)*(34+36)*(136+144)) > qchisq(1-0.05, 1)

#根據先前報告，70%的甲種雞的社會地位比乙種雞高，當兩種雞相處時，他們每小時每隻雞會鬥爭1.5次。如果我們將甲乙兩種雞各10隻關在一起，需要多少時間的觀察才可測出甲雞種的社會地位比乙種高? (要求檢定的低一類錯誤機率alpha= 0.05，第二類錯誤機率 Beta= 0.1)
p<- 0.7; q<- 1-p
p0<- 0.5 #虛無假設H0: P=0.5
alpha<- 0.05;beta<- 0.1; zscoreII<- qnorm(1-alpha/2); zscoreI<- qnorm(1-beta)
n<- (zscoreI + zscoreII)^2*p*q/ (p-p0)^2
ceiling(n) #至少需要鬥爭56次
fight<- 1.5*20/2 #每小時每隻雞鬥爭1.5次，共有20隻雞，但兩隻雞才是唯一次鬥爭
fight_real<- fight*10/(20-1) #實際跟不同種雞鬥爭的次數/小時，10隻不同種，9隻同種
ceiling(n)/fight_real

#A種雞勝C種雞的機率為0.6，B種雞勝C種雞的機率為0.4。這次將A、C種各10隻同圈飼養，B、C也各10隻同圈飼養，兩群雞各需多少時間的觀察才可測出A種雞比B種雞會打架? (鬥爭頻率、alpha、Beta值同上)
p1<- 0.6; q1<- 1-p1; p2<- 0.4; q2<- 1-p2 #虛無假設H0: P1=p2
alpha<- 0.05;beta<- 0.1; zscoreII<- qnorm(1-alpha/2); zscoreI<- qnorm(1-beta)
n<- (zscoreI + zscoreII)^2*(p1*q1 + p2*q2)/ (p1-p2)^2
fight_real<- fight*10/(20-1)
ceiling(n/fight_real)

#希望藉由民意測驗了解候選人可得到的選票，需要多少的樣本才可保證95%的結過不偏離實際結果5%?
L<- 0.05 ;p<-0.5;q<- 1-p #p和q為0.5時，pq值最大
alpha<- 0.05; zscore<- qnorm(1-alpha/2)
n<- zscore^2*p*q/L^2
ceiling(n)

#對於某法令，抽樣調查兩政黨人士的意見如下，兩政黨人士對於此法令支持程度是否有差異?
law<- matrix(c("G",T,49,
              "G",F,10,
              "B",T,39,
              "B",F,30),nrow= 4 ,ncol= 3,byrow=T, 
       dimnames= list(c(1:4),c("political party","opinion","number")))
law<- as.data.frame(law)
n<- sum(as.numeric(law$number))
(n*(49*30-10*39)^2) / ((49+39)*(10+30)*(10+39)*(39+30)) > qchisq(1-0.005, 1) #極顯著差異 (p < 0.005)

#雞場蛋雞患馬立克氏病的機率是0.5%，抽驗1000隻蛋雞，而不曾發現此病症的機率有多大?
dbinom(0, 1000, 0.005)#二項式分布
#p很小，n很大的情況也可以使用possion 分布
np<- 1000*0.005
dpois(0,np)

#假設台北市每星期因車禍死亡平均2.5人，連續兩周沒有人因車禍死亡的機率又多大?
#符合n很大，p很小，可以假設該變數符合poisson distribution
dpois(0, 2.5)*dpois(0, 2.5)

#在一野生動物族群內，帶有某一隱性突變性狀的個體(具有兩個隱性突變基因)佔有4%。根據遺傳學哈溫定律，如果此族群個體對性狀而言都以逢機交配，則沒有帶有該隱性基因的個體應佔有(1-sqrt(0.04))^2= 0.64，而僅帶有一個隱性基因沒有表現此隱性特質的個體應佔有32%。
#在此族群中此性狀表現正常之個體有多少比例帶有隱性基因?
p<- 0.96 ; q<- 1-p
aa<- q
AA<- (1-sqrt(q))^2
Aa<- 1-AA- aa
Aa/(Aa+AA)
#若父母此性狀均正常但帶有隱性基因，子代表現隱性特質的機率為0.25。如果父母之一未帶有隱性基因，則其子代不會表現隱性特質。在此族群中，逢機取兩個正常表現的個體交配，其子代會表現隱性特質的機率是多少?
(Aa/(Aa+AA))^2*0.25#取到正常且帶有隱性基因之機率，且父母都要是該情況，且子代表現機率是0.25

#肉牛中有些外表矮小的現象 (dwarfism)是由於遺傳上的隱性基因所造成的，由父母個得到一個矮小基因的個體就會有矮小的體型，僅有一個矮小基因仍有正常的體型，但一般稱為擔體 (carrier)。兩個單體交配產生的子代有0.75的機率仍是正常的。如果有一頭體型正常的優良公牛，將之與一群屬於擔體的母牛交配，需要生產多少頭正常體型的仔牛，才有99%的把握來認為此公牛並不是一頭擔體?
n<- log(0.01)/log(0.75)

#由過去長久的經驗得知肉品工廠的罐頭，大約有5%是不合格。如果有天生產的罐頭逢機取5個，至少有三個不合格的機率?
p<- 0.05; n<-5
Y_3<- dbinom(3,size= n,prob = p) #3個不合格
Y_4<- dbinom(4,size= n,prob = p) #4個不合格
Y_5<- dbinom(5,size= n,prob = p) #5個不合格
sum(Y_3,Y_4,Y_5)

#三小時內能100公釐以上雨水的大雨，在某城市一年內平均平均發生次數剛好一次，且符合poisson分布。請問下一年度發生這種大雨的次數超過4次的機率有多少?
1- ppois(4, 1) #累積發生四次的機率

#某獸醫研究者認為以剖腹生產的仔豬，隔離飼養，並餵飼以熱壓消毒的飼料，應該有50%的機率不會有某種病原菌。經試驗屠宰檢查43頭豬，其中有28頭沒有這種病原菌。
#常態近似法
p<- 0.5 ; n<- 43; p_bar<- 28/n; sd<- sqrt(n*p*(1-p))
expected<- n*p
(28-expected) / sd > qnorm(1-0.05/2) #感染率顯著高於50%
#chi-squares近似法
((28-expected)^2+ ((43-28)-expected)^2) / expected > qchisq(1-0.05, df=1) #結果同上

#50隻小鼠注射疫苗，另外50隻注射鹽水，之後都讓之感染病原。結果注射疫苗者有4隻死亡，注射鹽水者有13隻死亡。
#注射疫苗是否有效，有多少信心?
n1<- 50; x1<- 4; n2<- 50; x2<- 13; df<-4-1-1
chi<- ((n1+n2)*(x1*(n2-x2)-(n1-x1)*x2)^2) / (n1*n2*(x1+x2)*((n2-x2)+(n1-x1)))
pchisq(chi,df) #有94%信心水準
#如果對照組存活率是80%，若疫苗能使存活率提高至90%，希望能有80%的機率顯著的檢定出疫苗的效果，需要測定幾隻小鼠?
p1= 0.8 ; p2<- 0.9 ;alpha<- 0.05; beta<- 0.2
zscore1<- qnorm(1-alpha/2); zscore2<- qnorm(1-beta)
n<- (zscore1+zscore2)^2*(p1*(1-p1)+p2*(1-p2)) / (p1-p2)^2

#調查240頭乳牛，分娩後胎衣滯留子宮有108頭。而其中有76頭在分娩後90天內仍未配種成功。這些資料是否足夠顯示胎衣滯留與分娩後配種成功有關連?
240*(32*52- 76*80)^2 / (32+76)*(80+52)*(32+80)*(76+52) > pchisq(1-0.05, 1)

#檢定452頭公豬配種能力，並依品種與豬場規模分類
pig<- matrix(c("P","L",47,24,
               "P","S",43,15,
               "H","L",98,28,
               "H","S",154,43), nrow=4,ncol=4,byrow=T,
             dimnames = list(c(1:4),c("species","scale","good","bad")))
#分別計算純種豬與雜交豬配種能力合格率的95%信賴區間
p1 <-(47+43)/n1 ;n1<- (47+43+24+15); sd1<- sqrt(p1*(1-p1)/n1)
alpha = .05
zscore <- qnorm(1-alpha/2)
c(p1-zscore*sd1-1/(2*n1), p1+zscore*sd1+1/(2*n1))

p2 <-(98+154)/n2 ;n2<- (98+28+154+43); sd2<- sqrt(p2*(1-p2)/n2)
c(p2-zscore*sd2-1/(2*n2), p2+zscore*sd2+1/(2*n2))
#純種豬與雜交豬配種合格率是否不同?
#卡方檢定
chi_species<- (n1+n2)*((47+43)*(28+43)-(24+15)*(98+154))^2 / ((47+43+24+15)*(98+28+154+43)*(47+43+98+154)*(24+15+28+43))
pchisq(chi_species,1) #未顯著，但有趨勢

#豬場規模與配種合格率是否有關係?
chi_scale <- (452*((47+98)*(15+43)-(43+154)*(24+28))^2) / ((47+98+43+154)*(15+28+24+43)*(47+24+98+28)*(43+15+154+43))
pchisq(chi_scale,1) # 沒有差異

提供R script參考( 格式為 CP950)