วันพฤหัสบดีที่ 20 กุมภาพันธ์ พ.ศ. 2563

โปรแกรม R บทที่ 5


บทที่ 5
Probability and Distributions

5.1 Randon Sampling
            ในยุคต้นๆของทฤษฎีความน่าจะเป็น เรื่องที่ศึกษามักเป็นเรื่องเกี่ยวกับเกมส์และการพนัน ดังนั้นตัวอย่างสุ่มจึงมักเกี่ยวข้องกับไพ่ที่หยิบขึ้นมาจากสำรับ หรือการหยิบลูกบอลจากกล่อง เป็นต้นในR เราสามารถจำลองสถานการณ์ (simulation) เหล่านี้ได้ด้วยฟังก์ชันsample() ซึ่งมีรูปแบบดังนี้
sample(x, size, replace = FALSE, prob = NULL)
            เมื่อx คือเวคเตอร์ข้อมูลที่ประกอบด้วยค่าของข้อมูลที่ เป็นไปได้ทั้งหมด,size คือจำนวนตัวอย่างสุ่มที่ต้องการ,  replace=เป็น argument ที่ใช้เพื่อกำหนดว่าการสุ่มตัวอย่างที่กำลังจะทำนั้นเป็ นการสุ่มตัวอย่างแบบใส่คืน (with replacement) หรือไม่ใส่คืน (without replacement) และprob= เป็ argument ที่ใช้เพื่อกำหนดค่าความน่าจะเป็นของข้อมู ลแต่ละค่าในเวคเตอร์ x   หากเราต้องการสุ่มตัวเลข 5 ตัว จากตัวเลขชุดหนึ่งที่มีค่าตั้งแต่ 1 ถึง 40 สามารถทำได้ดังนี้
> sample(1:40, 5)
[1] 20 13 23 18 14
argument แรก1:40 หมายถึงเวกเตอร์ของตัวเลข 1 ถึง 40 จำนวน 40 ค่า ซึ่งถือเป็นประชากรของข้อมูล ส่วนargument ที่สองคือ 5 หมายถึงขนาดของตัวอย่างสุ่มที่เราต้องการ ในความเป็นจริงเราสามารถใช้ sample(40,5)ก็ได้ซึ่งจะให้ผลเหมือนกัน เนื่องจากตัวเลข 40 หมายถึงช่วงของเลขจำนวนเต็มที่มีค่าตั้งแต่ 1 ถึง 40 ให้สังเกตว่า default ของฟังก์ชันsample() จะเป็ นการสุ่มแบบไม่ใส่คืน ดังนั้นหากต้องการให้เป็ นการสุ่มแบบใส่คืนสามารถเปลี่ยนได้โดยเปลี่ยน argumentreplace=TRUE เช่น การจำลองข้อมูลของผลของการโยนเหรียญ 1เหรียญ 10 ครั้ง
> sample(c("H", "T"), 10, replace=TRUE)
[1] "T" "T" "T" "H" "H" "H" "T" "H" "H" "T"
            การจำลองข้อมูลในตัวอย่างข้างต้นเป็นกรณีที่ ความน่าจะเป็นในการเกิดแต่ละผลลัพธ์เท่ากัน แต่เราสามารถจำลองข้อมูลที่ความน่าจะเป็นของการเกิดแต่ละผลลัพธ์ไม่เท่ากันได้ด้วย เช่น การจำลองข้อมูลผลการผ่าตัด ซึ่งผลของการผ่าตัดมีสองแบบคือ สำเร็จ หรือล้มเหลว โดยที่ความน่าจะเป็ นที่การผ่าตัดจะประสบความสำเร็จคือ 90 เปอร์เซ็นต์ ทำได้ดังนี้
> sample(c("succ", "fail"), 10, replace=TRUE, prob = c(0.90, 0.10))
[1] "succ" "succ" "succ" "succ" "succ" "succ" "succ" "succ" "succ" "succ"
            จากที่เราทราบดีอยูแล้วว่า เราสามารถอธิบายลักษณะของตัวเลขสุ่มได้ด้วยการแจกแจง (distribution) ซึ่งการแจกแจงก็คือฟังก์ชันที่ ใช้ในการหาค่าความน่าจะเป็นเมื่อตัวแปรสุ่มมีค่าใดค่าหนึ่ง หรือมีค่าอยู่ ในช่วงใดๆ เช่นP(a 67 X b)http://htmlimg4.scribdassets.com/4g6kn2nghsma7fi/images/68-b3d4a48837/000.jpg เป่นต้นซึ่งในกรณีของตัวแปรสุ่มแบบต่อเนื่องฟังก์ชันการแจกแจงของตัวแปรสุ่มจะอยู่ในรูปของฟังก์ชัน
ความหนาแน่น (
density function) แต่ถ้าตัวแปรสุ่มเป็นแบบไม่ต่อเนื่องฟังก์ชันจะอยู่ในรูปของฟังก์ชันการแจกแจง(distribution function) เช่นP(X= k) เป็ นต้น ซึ่งR มีคำสั่งที่ใช้ในการสุ่มค่าสังเกตหรือข้อมูลที่มีการแจกแจงแบบต่างๆได้ ต่อไปนี้เป็นตัวอย่างของการใช้ R ในการสุ่มตัวอย่างจากการแจกแจงแบบต่างๆ
Uniform.รูปแบบของคำสั่งที่ใช้ในการสุ่มตัวอย่างจากการแจกแจงแบบ uniform คือ
runif(n, min=0, max=1)
ตัวอย่างเช่น หากเวลาในการรอสัญญานไฟจราจรมีการแจกแจงแบบเอกรูป (uniform distribution) ในช่วง [0,2]
> runif(1,0,2)                                                     # time at light
[1] 0.7325975
> runif(5,0,2) # time at 5 lights
[1] 0.5658379 0.1512934 1.3920678 1.6725483 0.9159859

> runif(5)                                                          # 5 random number in [0,1]
[1] 0.3441802 0.2687535 0.5044450 0.8939212 0.7787676
เพื่อแสดงให้เห็นถึงลักษณะการแจกแจงแบบ uniform ที่มีค่าต่ำสุดเท่ากับ 0 และค่าสูงสุดเท่ากับ 1 ทำได้ดังนี้
> x = runif(100)
> hist(x,prob=TRUE, col=gray(.9), main="Uniform on [0,1]")
> curve(dunif(x,0,1),add=TRUE)
           
Normal.การแจกแจงแบบปกตินั้นมีพารามิเตอร์ที่เกี่ยวข้องอยู่ 2 ตัว คือµ และσ เช่น IQ มีการแจกแจงปกติด้วยแจงแบบปกติโดยมีค่าเฉลี่ ยเท่ากับ 280 และค่าเบี่ ยงเบนมาตรฐานเท่ากับ 10 เราสามารถทำให้ค่าของตัวแปรสุ่มปกตินี้เป็ นค่ามาตรฐานZ= (X µ)/σได้ ซึ่งค่ามาตรฐานนี้จะมีการแจกแจงแบบปกติมาตรฐานที่มีค่าเฉลี่ยเท่ากับ 0 และค่าเบี่ยงเบนมาตรฐานเท่ากับ
rnorm(n, mean=0, sd=1)
ตัวอย่างต่อไปนี้แสดงการสุ่มตัวอย่างจากการแจกแจงแบบปกติ
> rnorm(1,100,16)                                              # an IQ score
[1] 122.7393

> rnorm(1,mean=280,sd=10)                               # how long for a baby (10 days early)
[1] 276.9401
            เพื่อให้เห็นถึงลักษณะของการแจกแจงแบบปกติที่มีค่าเฉลี่ย
0 และค่าเบี่ยงเบนมาตรฐานเท่ากับ 1 ทำได้ดังนี้
> x = rnorm(100)
> hist(x,prob=TRUE, col= "blanchedalmond", main="Normal(0,1)")
> curve(dnorm(x),add=TRUE)
ได้กราฟดังรูปที่
5.2

รูปที่
5.2: Normal(0,1)
           
Binomial.ตัวแปรสุ่มทวินามเป็ นตัวแปรสุ่มแบบไม่ต่อเนื่อง โดยที่ค่าของตัวแปรสุ่มทวินามคือจำนวนของความสำเร็จที่เกิดขึ้นในการทดลองที่ เป็นอิสระกัน n ครั้ง เราสามารถสุ่มตัวอย่างจากการแจกแจงทวินามได้โดยใช้คำสั่ง  rbinom(n, size, prob) ในการสร้างตัวเลขสุ่ มทวินามทำได้ดังนี้
> n = 10; p=.5
> rbinom(1,n,p)                                                 # 6 successes in 10 trials
[1] 4
> rbinom(5,n,p)
                                                 # 5 binomial numbers
[1] 5 5 4 5 6
            อย่างไรก็ตามเมื่อ
n มีขนาดใหญ่ตัวแปรสุ่มทวินามจะมีลักษณะเข้าใกล้ตัวแปรสุ่มปกติ ซึ่งเป็ นผลมาจากทฤษฎี ขีดจำกัดกลาง (Central limit theorem) โดยจะได้ว่า
             

            จะมีการแจกแจงประมาณได้ด้วยการแจกแจงแบบปกติ
รูปที่
5.3 แสดงการแจกแจงของตัวแปรสุ่มทวินามที่ n= 5,15,50 เมื่อp= 0.25 ซึ่งจากกราฟจะเห็นได้ว่าเมื่อ
n มีขนาดเพิ่มขึ้นรูปร่างของการแจกแจงทวินามจะใกล้เคียงกับลักษณะของการแจกแจงแบบปกตินั่นคือมี
ลักษณะเป็นรูประฆังคว่ำ สมมาตรซึ่งกราฟเหล่านี้ได้จากคำสั่ง
par(mfrow=c(1,3))
n=5;p=.25
x=rbinom(100,n,p)

hist(x,prob=TRUE,col="lightblue", main="n=5, p=.25")
xvals=0:n
points(xvals,dbinom(xvals,n,p),type="h",lwd=2)
points(xvals,dbinom(xvals,n,p),type="p",lwd=2)
#++++++++++++++++++++++++++++++++++++
n=15;p=.25
x=rbinom(100,n,p)

hist(x,prob=TRUE,col="lightgreen", main="n=15, p=.25")
xvals=0:n
points(xvals,dbinom(xvals,n,p),type="h",lwd=2)
points(xvals,dbinom(xvals,n,p),type="p",lwd=2)

#++++++++++++++++++++++++++++++++++++
n=50;p=.25
x=rbinom(100,n,p)

hist(x,prob=TRUE,col="brown1", main="n=50, p=.25")
xvals=0:n
points(xvals,dbinom(xvals,n,p),type="h",lwd=2)
points(xvals,dbinom(xvals,n,p),type="p",lwd=2)


Exponential.การแจกแจงแบบ exponentail เป็นการแจกแจงที่สำคัญอีกการแจกแจงหนึ่ง ซึ่งมักใช้ในการอธิบาย
ถึงลักษณะการแจกแจงของอายุการใช้งานของอุปกรณ์ไฟฟ้ า เช่น ถ้าค่าเฉลี่ยของอายุการใช้งานของหลอดไฟเท่าากับ
2500 ชั่วโมง เราอาจจะคิดได้ว่าอายุการใช้งานนี้มีการแจกแจงแบบ exponentail ที่มีค่าเฉลี่ยเท่ากับ 2500 พารามิเตอร์สำหรับการแจกแจงนี้คือrate = 1/mean รูปแบบคำสั่งสำหรับการสุ่มตัวอย่างที่มีการแจกแจงแบบ exponentailคือ
rexp(n, rate = 1)
            ตัวอย่างต่อไปนี้เป็นการสุ่มตัวอย่างข้อมูลที่มีการแจกแจงแบบ
exponentail ที่มี rate เท่ากับ 1/2500 และได้กราฟของการแจกแจงของข้อมูลชุดนี้ดังรูปที่ 5.4
> x = rexp(100,1/2500)
> hist(x, prob=TRUE, col="yellow3", main="Exponentail mean=2500")
> curve(dexp(x,1/2500),add=TRUE)

รูปที่
5.3: Random binomial data with the theoretical distribution

รูปที่
5.4: Random exponential data with the theoretical density
5.2 Probability Calculation and Combinatorics  
            ในกรณีที่การสุ่มตัวอย่างเป็ นแบบไม่ใส่คืน เช่นsample(1:40), 5 เราทราบว่าความน่าจะเป็ นของการที่จะได้ค่าใดค่าหนึ่งเป็นค่าแรกคือ1/40 ค่าที่ สองเป็ น1/39 ตามลำดับ ดังนั้นในตัวอย่างนี้ความน่าจะเป็ นที่ จะได้ตัวอย่างสุ่มชุดใดชุดหนึ่งจะเท่ากับ1/(40× 39× 38××37× 36) หรือในR เราสามารถใช้ฟังก์ชันprod() ในการคำนวณหาค่าความน่าจะเป็นนี้ได้ดังนี้
> 1/prod(40:36)
[1] 1.266449e-08
อย่างไรก็ตามจะสังเกตเห็นได้ว่าค่าความน่าจะเป็นที่จะได้หมายเลขชุดใดๆนั้นเกิดขึ้นตามลำดับที่กำหนด (
nPr) ถ้าในกรณี ที่ ลำดับของหมายเลขไม่ สำคัญหมายความว่าตัวเลขชุดเดียวกันถึงแม้ว่าจะมีลำดับการเรียงที่แตกต่างกันก็ ยังถือว่าเป็นตัวเลขชุดเดียวกัน (nCr) เราสามารถหาความน่าจะเป็นได้ดังนี้
> prod(5:1)/prod(40:36)
[1] 1.519738e-06
หรือใช้ฟังก์ชัน
choose() ใน R คำนวณคานี้ได้ดังนี้
> 1/choose(40,5)
[1] 1.519738e-06
5.3 Buil-In Distributions inR
            ถูกสร้างเป็นฟังก์ชันไว้ในR เรียบร้อยแล้ว เราจึงสามารถใช้ฟังก์ชันเหล่านี้ในการหาค่าความน่าจะเป็นหรือค่าวิกฤตสำหรับการทดสอบสมมติฐานได้ โดยใช้แทนตารางสถิติแบบดั้งเดิม ในที่นี้เราจะพิจารณาการแจกแจงแบบปกติ และการแจกแจงแบบทวินาม อย่างไรก็ตามสำหรับการแจกแจงแบบอื่นๆก็มีรูปแบบเดียวกัน
ต่อไปนี้เป็ นสิ่งที่เราสามารถหาได้จากการแจกแจงทางสถิติ
• Density หรือ point probability (ขึ้นอยู่กับเป็ นการแจกแจงไม่ต่อเนื่องหรือต่อเนื่อง)
• Cumulated probability, distribution function
• Quantiles
• Pseudo-random numbers
            Rมีฟังก์ชันสำหรับหาค่าเหล่านี้ได้ เช่น สำหรับการแจกแจงปกติจะมีฟังก์ชัน dnorm(), pnorm(), qnorm()และ rnorm()ตามลำดับ (density,probability,quantile, และrandom)
           
5.3.1 Densities
            Density สำหรับการแจกแจงแบบต่อเนื่องคือค่าของความน่าจะเป็ นที่ ตัวแปรสุ่มจะมีค่าเข้าใกล้ค่า x
ซึ่ งความน่าจะเป็นที่ตัวแปรสุ่มจะมีค่าอยู่ในช่วงที่กำหนดนี้คือพื้นที่ใต้กราฟที่สอดคล้องกับการแจกแจงของตัวแปรสุ่มนี้  สำหรับการแจกแจงแบบไม่ต่อเนื่อง คำว่า
density หมายถึงความน่าจะเป็ นแบบจุด ซึ่งคือความน่าจะเป็นที่ตัวแปรสุ่มจะมีค่าเท่ากับ x density function เป็นการแจกแจงปกติ สามารถทำได้ดังนี้
> x = seq(-4,4,0.1)
> plot(x,dnorm(x), type="l")
            ไดกราฟดังรูปที่
5.5 ซึ่ง argument ค่า -4 ถึง 4 และ 0.1 หมายถึงเราสร้างให้เวกเตอร์x ประกอบด้วยข้อมูลที่มีค่าตั้งแต่ -4 ถึง 4 โดยที่แต่ละค่าห่างกัน 0.1 นั่นคือเวกเตอร์ x ประกอบด้วย (4.0, 3.9, 3.8,..., 3.8, 3.9, 4.0) ในฟังก์ชันplot() เราใช้ argumenttype = "l" เพื่อให้มีการสร้างกราฟเป็นอีกวิธีหนึ่ งที่ สามารถสร้างกราฟในแบบเดียวกันนี้ได้คือใช้ฟงก์ชันcurve() ดังนี้

รูปที่
5.5: Density of normal distribution
> curve(dnorm(x), from=-4, to=4)
คำสั่งนี้จะสะดวกมากกว่าในการสร้างกราฟรูปแบบนี้ แต่ ฟังก์ชันcurve() ต้องการค่าของy ในรูปฟังก์ชันของx
สำหรับการแจกแจงแบบไม่ต่อเนื่ อง ค่าของตัวแปรสุ่มจะเป็
นค่ามีจุดตัดที่ ชัดเจน จึงนิยมที่ จะสร้างกราฟแบบเส้นตรงตั้งฉากกับแกนx เช่นเราต้องการกราฟของการแจกแจงทวินามที่มีค่าn= 50 และp= 0.33 ทำได้ดังนี้
> x = 0:50
> plot(x,dbinom(x,size=50,prob=0.33), type="h")
            ได้กราฟดังรูปที่
5.6 ในฟังก์ชันdnorm() นี้เราใช้ argument 3 arguments โดย argument แรกคือx ซึ่งเป็น
ค่าของตัวแปรสุ่ม
argument ที่สองคือsize=50 เป็นการกำหนดขนาดของการทดลอง และ argument ที่สามคือ
prob=0.33เป็ นการกำหนดความน่าจะเป็ นของการเกิดความสำเร็จในการทดลองแต่ละครั้ง ส่วน argument type
= "h"ในฟังก์ชัน plot()เป็ นการสร้างกราฟเส้นตรงตั้งฉากกับแกนx
           
5.3.2 Cumulative Distribution Functions
            ฟังก์ชันการแจกแจงความนาจะเป็ นสะสม (cumulative distribution function) หมายถึงค่าความน่าจะเป็นที่ตัวแปรสุ่มมีค่าเท่ากับหรือน้อยกว่าค่าที่กำหนด ฟังก์ชันในR ที่ใช้ในการคำนวณค่าความน่าจะเป็ นสะสม
คือฟังก์ชันที่ขึ้นต้นด้วย
"p" เช่นpnorm และpbinom เป็ นต้น
            หากเราทราบว่าตัวแปรสุ่มตัวหนึ่งมีการแจกแจงแบบปกติ มีค่าเฉลี่ยเท่ากับ
132 และค่าเบี่ยงเบนมาตรฐานเท่ากับ 13 และอยากทราบว่าความน่าจะเป็ นที่ตัวแปรสุ่มนี้จะมีค่ามากกว่า 160 จะเท่ากับเท่าไร ทำได้ดังนี้

> 1 - pnorm(160, mean = 132, sd=13)
[1] 0.01562612

รูปที่
5.6: Point probabilities in binom(50,0.33)
            จากค่าความน่าจะเป็ นที่ได้หมายความว่ามีเพียง 1.5 เปอร์เซ็นต์ ของประชากรเทานั้นที่มีค่ามากกว่า 160 ซึ่งฟังก์ชัน pnorm()จะให้ค่าความน่าจะเป็ นที่ตัวแปรสุ่มจะมีค่าน้อยกว่าหรือเท่ากับค่าที่กำหนดไว้ใน argument แรกในฟังก์ชันด้วยค่าเฉลี่ยและค่าเบี่ยงเบนมาตรฐานที่กำหนด  โดยปกติเราจะนำการหาค่าความน่าจะเป็ นนี้ไปใช้ในการทดสอบสมมุ ติฐานทางสถิติต่างๆ เช่น ใน Simple sign test:มีคนไข้ 20 คน แต่ละคนได้รับทรีทเมนต 2 ชนิด และผูทดลองได้ถามคนไข้แต่ละคนว่าชอบทรีทเมนต์ใดมากกว่ากันปรากฏว่ามีคนไข้ 16 คน ที่ชอบทรีทเมนต์ A มากกว่า จากข้อมูลที่ได้เราต้องการทดสอบว่าทรีทเมนต์ A ดีกว่าทรีทเมนต์ B อย่างมีนัยสำคัญทางสถิติหรือไม่ ซึ่งถ้าทรีทเมนต์ทั้งสองดีพอๆกัน เราย่อมคาดหมายว่าจำนวนคนไข้ชอบทรีทเมนต์ A จะมีการแจกแจงแบบทวินามด้วยค่าp= 0.5 และn= 20 ดังนั้นความน่าจะเป็ นที่ จำนวนคนไข้ที่ ชอบทรีทเมนตA จะน้อยกว่าหรือเทากับ 16 คน คือ
> pbinom(16, size=20,prob=0.5)
[1] 0.9987116
            เราต้องการค่าความน่าจะเป็ นที่
X16 ถ้าเรานำค่าความน่าจะเป็ นที่ ได้นี้ไปลบออกจาก 1 เราจะไม่ได้ค่าความน่าจะเป็ น ที่ถูกต้อง เราควรทำดังนี้
> 1-pbinom(15,size=20,prob=0.5)
[1] 0.005908966
            หากต้องการหาค่าความน่าจะเป็ นของการทดสอบสมมติฐานแบบสองทาง เช่น ความน่าจะเป็
นที่ มีจำนวนคนไข้ที่ ชอบทรีทเมนต์ A มากกว่าหรือเท่ากับ 16 คน หรือน้อยกว่าเท่ากับ 4 คน
> 1-pbinom(15,size=20,prob=0.5) + pbinom(4,size=20,prob=0.5)
[1] 0.01181793
            ซึ่งเท่ากับสองเท่าของความน่าจะเป็ นทางเดียว
           
5.3.3 Quantiles
           
Quantile function เป็ นสวนกลับของฟังก์ชันความน่าจะเป็ นสะสม นั่นคือp-quantile คือค่าที่มีคุณสมบัติที่ทำให้ค่าความน่าจะเป็ นที่ ตัวแปรสุ่มจะมีค่าเท่ากับหรือน้อยกว่าค่าๆนี้เท่ากับp เช่น ตามนิยามแล้วค่ามัธยฐานจะเป็ นค่าควอนไทล์ที่ 50%
            โดยทั่วไปตารางของการแจกแจงทางสถิติมักจะให้ค่าที่อยู่ในรูปของควอนไทล์ สำหรับชุดของค่าความน่าจะเป็นที่กำหนด ตารางสถิติจะแสดงขอบเขตที่ค่าสถิติทดสอบต้องมีค่ามากกว่าค่านี้เพื่อการปฏิเสธสมมติฐานหลักที่ระดับนัยสำคัญที่ กำหนด  ในทางทฤษฎีมักใช้ค่าควอนไทล์ในการคำนวณหาช่วงความเชื่ อมั่นของพารามิเตอร์
เช่น ต้องการหาช่วงความเชื่ อมั่น 95%ของµ จาก
x+ σ/n× N0.025 µ¯x+ σ/n× N0.975
            เมื่อN0.025 คือควอนไทล์ที่2.5% ในการแจกแจงปกติ ถ้าσ= 12 และn= 5 และเราหาค่าเฉลี่ย
x= 83 เราสามารถหาช่วงความเชื่ อมั่นได้ดังนี้
> xbar = 83
> sigma = 12
> n = 5
> sem = sigma / sqrt(n)
> sem
[1] 5.366563
> xbar + sem * qnorm(0.025)
[1] 72.48173
> xbar + sem * qnorm(0.975)
[1] 93.51827

ดังนั้นช่วงความเชื่อมั่น
95% ของµ คือ 72.48 ถึง 93.52
โดยปกติ ค่าควอนไทล์มักเขียนอยู่ ในรูป
Φ−1 (0.975) เมื่อΦ คือสัญลักษณ์มาตรฐานของฟังก์ชันการแจกแจง
ความน่ าจะเป็
นสะสมของการแจกแจงแบบปกติ
5.4 Simulations
            หากเราสามารถจำลองข้อมูลสุ่มชนิดต่างๆได้จะทำให้เราสามารถทำการทดลองเพื่อตอบคำถามเกี่ยวกับการแจกแจง ของขอมูลได้เนื่องจากมีฟังก์ชันมากมายที่สามารถใช้ในการสร้างตัวเลขสุ่มได้ เราจะสามารถทราบถึงลักษณะการแจกแจงของค่าตัวเลขสุ่มเหล่านี้ได้โดยพิจารณาจาก histogram หรือวิธีการอื่นๆ ในที่นี้เราต้องการสร้างตัวเลขสุ่มชนิดใหม่ขึ้นมาและตรวจสอบว่ าการแจกแจงของตัวเลขสุ่ มนี้เป็ นการแจกแจงแบบใด
           
5.4.1 The Central Limit Theorem
           
เราเริ่มต้นด้วยตัวอย่างของ central limit theorem (CLT) ซึ่งเป็ นทฤษฎีที่สำคัญมากในทางสถิติ ซึ่งทฤษฎีนี้ได้กล่าวไว้ว่า หากXi ถูกสุ่มมาอย่างเป็ นอิสระกันจากประชากรที่มีค่าเฉลี่ยµ และค่าเบี่ยงเบนมาตรฐานσ ที่ทราบค่าแล้วค่ามาตรฐานของค่าเฉลี่ยของX คือ
                                   

            จะมีการแจกแจงเข้าใกล้การแจกแจงแบบปกติที่มีค่าเฉลี่ยเท่ากับ
0 และค่าเบี่ยงเบนมาตรฐานเท่ากับ 1 นั่นคือหากn มีขนากใหญ่พอแล้ว ค่าเฉลี่ยของตัวอย่างจะมีการแจกแจงประมาณได้ด้วยการแจกแจงแบบปกติ โดยมีค่าเฉลี่ยเท่ากับ
                                   

จะมีการแจกแจงประมาณได้ด้วย
N(0,1) และเราจะใช้R ในการตรวจสอบทฤษฎีนี้ดังต่ อไปนี้
> n = 10
> p=.25
> s=rbinom(100,n,p)
> x = (s - n*p)/sqrt(n*p*(1-p))
> hist(x,prob=TRUE,col="gold2")

> curve(dnorm(x), col="blue", add=TRUE)
            จากคำสั่งในข้างต้นจะได้กราฟดังรูปที่
5.7 ซึ่งจะเห็นได้ว่าการแจกแจงของข้อมูลนั้นใกล้เคียงกับN(0,1)

รูปที่
5.7: Scaled binomial data is approximatelyN(0,1)
5.4.2 For Loops
โดยทั่วไปกระบวนการในการสร้างตัวเลขสุ่มจำนวน 100 ค่า อาจทำได้ไม่ง่ายนัก ดังนั้นเราอาจจำเป็นต้องที่จะสร้าง
ตัวเลขสุ่มนี้ที่ละค่า ซึ่งทำได้โดยการใช้ "for loop" ตัวอย่างต่อไปนี้แสดงการสร้างตัวเลขสุ่มจำนวน 100 ค่า ที่ละค่า
> results = numeric(0)
# a place to store the results
> for (i in 1:100){
# the for loop
+ s = rbinom(1,n,p)
# just 1 this time
+ results[i] = (s - n*p)/sqrt(n*p*(1-p))
# store the answer
+ }
> hist(results)
จากคำสั่งเหล่านี้เราได้สร้างตัวแปรresults ไว้สำหรับเก็บข้อมูลที่ เราสร้างขึ้น เราสามารถดูลักษณะของข้อมูลที่ เรา
สร้างขึ้นได้โดยใช้คำสั่ง hist(results) ซึ่งได้ดังรูปที่ 5.8
CLT with normal data.CLT สามารถใช้กับตัวแปรสุ่มปกติได้ด้วย พิจารณาตัวอย่างต่อไปนี้ ให้Xi เป็นตัวเลขสุ่
ปกติที่มีค่าเฉลี่ยµ= 5 และσ= 5
> results=c()
> mu=0;sigma=1
> for(i in 1:200) {
+ x=rnorm(100,mu,sigma)

รูปที่
5.8: Simulation of CLT with binomial data
+ results[i]=(mean(x) - mu)/(sigma/sqrt(100))
+ }
> hist(results,prob=TRUE,ylim=c(0,.4), col="green")
> curve(dnorm(x), add=TRUE)
จากคำสั่งข้างบนได้กราฟดังรูปที่ 5.9 โดย histogram ที่ได้ชี้ให้เห็นว่าข้อมูลมีการแจกแจงประมาณไดด้วยการแจกแจง
ปกติ
Histogram of results

รูปที่
5.9: Simulation of CLT with normal data
5.4.3 Normal Probability Plots
กราฟอีกชนิดหนึ่งที่ สามารถนำมาใชในการพิจารณาวาขอมู ลมี การแจกแจงแบบปกติ หรือไม ไดแกnormal probability
plotซึ่งกราฟชนิดนี้ใชไดดีกวา histogram ในการตรวจสอบลักษณะการแจกแจงปกติของขอมูล แนวคิดของการสราง
กราฟรูปแบบนี้มาจากการพลอตควอนไทลของขอมูลของเรากับควอนไทลของการแจกแจงแบบปกติ ซึ่งควอนไทลที่qคือคาของขอมู ลที่ ขอมู ลจำนวนq×100% ของขอมู ลทั้งหมดที่ มีคานอยกวาคานี้ ดังนั้นควอนไทล 0.25 คือ ควอไทลที่1, ควอนไทล 0.5 คือควอไทลที่ 2 หรือ median และควอนไทล 0.75 คือควอไทลที่ 3 ควอนไทลของการแจกแจงทางทฤษฎีก็คลายคลึงกันเพียงแต แทนที่ จะเป นจำนวนของขอมูลที่ มีคานอยกวา ก็กลายเป นพื้นที่ ที่อยู ทางดานซายของคาที่กำหนด เชน median จะแบงพื้นที่ใตโคงความหนาแนน (density curve) ออกเป นสองสวนเทาๆกัน
การแปลความหมายของ normal probability plot ทำไดโดยพิจารณาการกระจายของจุดบนกราฟ หากการกระจายมีแนวโนมเป นเสนตรงแลวแสดงวาขอมูลมีการแจกแจงประมาณไดดวยการแจกแจงแบบปกติ หากการกระจายเป นเสนโคงแสดงว าการแจกแจงนั้นมีหางยาวหรือหางสั้นได
ในR มีฟงกชันที่ใชในการสราง normal probability plot ไดแกqqnorm() หรือฟงกชันทั่วไปคือqqplot()
และqqline() ที่ ใชในการลากเสนอางอิง
ตัวอยางตอไปนี้แสดงการสราง normal probability plot สำหรับขอมูลที่มีการแจกแจงแบบตางๆ
> par(mfrow=c(1,4))
> x = rnorm(100);qqnorm(x, main="Normal(0,1)");qqline(x)
> x = rnorm(100,10,15);qqnorm(x, main="Normal(10,15)");qqline(x)
> x = rexp(100,1/10);qqnorm(x, main="Exponential mu = 10");qqline(x)
> x = runif(100,0,1);qqnorm(x, main="Uniform(0,1)");qqline(x)
ไดกราฟดังรูปที่ 5.10 จะเห็นไดวากราฟ 2 รูปแรกมีลักษณะของการกระจายคลายเสนตรง ในขณะที่กราฟสองรูปหลัง
ไม

รูปที่
5.10: Some normal plots

ไม่มีความคิดเห็น:

แสดงความคิดเห็น