[关闭]
@Macux 2015-12-01T06:49:06.000000Z 字数 6760 阅读 2284

R语言_概率篇

R语言_学习笔记

Part 1:随机数

1、生成服从某一分布随机数

  1. > Fox <- runif(100,min=-19,max=80)
  2. # 生成服从均匀分布的100个随机数，且均匀分布的左边界为-19,右边界为80
  1. > Fox <- rnorm(200,mean=18,sd=76)
  2. # 生成服从正太分布的200个随机数，且正态分布的均值为18,标准差为76

上面两个是常用的两个生成随机数的方法。当然也可以玩点小花样出来：

  1. > sds <- sample(c(0,1,2), 100, replace=TRUE) 
  2. # replace=TRUE表示是放回抽样。
  3. > means <- rchisq(100,df=10)
  4. > Fox <- rnorm(1000,mean=means,sd=sds)

如果只想生成正数或负数，调整一下标准差或均值即可。

  1. > Fox <- rwilcox(100,m=10,n=20)

当想生成整数随机数列，此分布是最佳选择。

  1. > Fox <- rpois(100,lambda=300)  #lambda表示均值。
  1. > Fox <- rchisq(100,df=9)
  1. > Fox <- rexp(100,rate=0.4)   #rate表示发生率
  1. > Fox <- rbeta(100,shape1=2,shape2=10)
  1. > Fox <- rf(100,df1=4,df2=8)
  1. > Fox <- rgamma(100,3,7)
  1. > Fox <- rlnorm(100,meanlog=0,sdlog=2)
  1. > Fox <-  rt(100,df=3)
  1. > Fox <- rweibull(100,shape=3,scale=4)

2、生成可再生的随机数
好处：当我每次运行某一脚本或程序时，都是用同一个序列。对比一下：

  1. > runif(5)
  2. [1] 0.02626516 0.11781895 0.55625933 0.76238102 0.34511263
  3. > runif(5)
  4. [1] 0.9751362 0.6268066 0.8686730 0.4429201 0.5387621
  6. > set.seed(178)
  7. > runif(5)
  8. [1] 0.97328778 0.30880372 0.78763989 0.03226390 0.09998702
  9. > set.seed(178)
  10. > runif(5)
  11. [1] 0.97328778 0.30880372 0.78763989 0.03226390 0.09998702

3、sample()

  1. #随机地从一个向量中选定n项，不放回抽样
  2. > sample(c(1,2,3,4,5,6,7,8,9,10),5)
  3. [1] 5 2 1 3 6
  1. #随机地从一个向量中选定n项，放回抽样
  2. > sample(c(1,2,3,4,5,6,7,8,9,10),5,replace=TRUE)
  3. [1] 10  4  2  8  4
  1. #以随机顺序选定vec的所有元素，并且每个元素只使用一次。
  2. #相当于随机地对向量vec重新排序。
  3. > sample(1:9)
  4. [1] 8 3 1 5 2 9 4 6 7

Part 2:假设检验

1、基础知识点

  1. 当P<α(显著水平)时，拒绝H0，结论表述为：很可能不是...
  2. 当P>α(显著水平)时，接受H0，结论表述为：很可能是...
  1. 如果在一次大选中，某人的支持率为55%，置信水平0.95上的置信区间是（50%,60%），那么他的真实支持率有95%的概率落在50%和60%之间。因此他的真实支持率不足50%的可能性小于2.5%（假设分布是对称的）。

2、计算相对频数
用一个逻辑表达式，来识别感兴趣的观测值。例如

  1. > lens <- rt(200,df=3)
  2. > mean(lens>0)
  3. [1] 0.525

3、数据集的分位数及分位数的逆

  1. #寻求观测值x，使得小于x的观测值比例为f
  2. > lens <- rt(200,df=3)
  3. > quantile(lens,0.07)
  4.        7% 
  5. -1.725398 
  7. #此函数更常用的用法：
  8. #1、找出数据集的四分位数：
  9. > lens <- rt(200,df=3)
  10. > quantile(x)
  11.          0%         25%         50%         75%        100% 
  12. -4.76590614 -0.67000543  0.05765649  0.87783149  8.66798922 
  14. #2、识别数据集中间某一比例（例如90%）的观测值
  15. > lens <- rt(200,df=3)
  16. > quantile(x,c(.05,.95))
  17.        5%       95% 
  18. -1.998840  2.188455 
  1. #给定观测值x,想知道小于x的数据的比例。
  2. > lens <- rt(200,df=3)
  3. > mean(lens < -1.725398)
  4. [1] 0.07

4、t检验

  1. > novo <- rgamma(50,2,5)
  2. > t.test(novo,conf.level=0.99)
  4.     One Sample t-test
  6. data:  novo
  7. t = 9.5248, df = 49, p-value = 9.879e-13
  8. alternative hypothesis: true mean is not equal to 0
  9. 99 percent confidence interval:
  10.  0.2558597 0.4562121
  11. sample estimates:
  12. mean of x 
  13. 0.3560359 
  14. # 输出结果表明，置信水平为99%时，总体均值的置信区间为(0.2558597,0.4562121)
  1. # 现有两个来自不同总体的样本，利用t检验，分析这两个总体是否有相同的均值。
  2. > lens1 <- rweibull(100,shape=3,scale=4)
  3. > lens2 <- rt(200,df=1)
  4. > t.test(lens1,lens2)
  6.     Welch Two Sample t-test
  8. data:  lens1 and lens2
  9. t = -0.0592, df = 199.472, p-value = 0.9528
  10. alternative hypothesis: true difference in means is not equal to 0
  11. 95 percent confidence interval:
  12.  -7.170954  6.752886
  13. sample estimates:
  14. mean of x mean of y 
  15.  3.394187  3.603221 
  16. # 若是配对数据，可以设置参数paired=TRUE。若两个总体有相同的方差，可以设置参数var.equal=TRUE。
  17. # 若其中一个样本量少于20个数据点，则其必须服从正态分布。
  18. # 若两个总体有相同的方差，可以设置参数var.equal=TRUE使检验更有效。
  1. # 现有一个来自总体的样本，利用t检验分析总体的均值是否等于mu
  2. > lens1 <- rweibull(100,shape=3,scale=4)
  3. > t.test(lens1,mu=3.3)
  5.     One Sample t-test
  7. data:  lens1
  8. t = 0.7751, df = 99, p-value = 0.4402
  9. alternative hypothesis: true mean is not equal to 3.3
  10. 95 percent confidence interval:
  11.  3.153061 3.635313
  12. sample estimates:
  13. mean of x 
  14.  3.394187 
  15. # 由于p=>0.05，所以有证据（样本数据）支持原假设：总体均值为3.3

5、中位数的置信区间

  1. > lens1 <- rweibull(100,shape=3,scale=4)
  2. > wilcox.test(lens1,conf.int=TRUE)
  4.     Wilcoxon signed rank test with continuity correction
  6. data:  lens1
  7. V = 5050, p-value < 2.2e-16
  8. alternative hypothesis: true location is not equal to 0
  9. 95 percent confidence interval:
  10.  3.118144 3.632256
  11. sample estimates:
  12. (pseudo)median     #伪中位数，它与真正的中位数是有区别的。
  13.       3.388438 
  14. > median(lens1)
  15. [1] 3.459298

6、比例的置信区间

  1. # 有一个来自总体的样本数据，它由成功或失败两种时间构成，现需分析总体成功比例的置信区间
  2. > prop.test(6,9,alternative="greater")
  4.     1-sample proportions test with continuity correction
  6. data:  6 out of 9, null probability 0.5
  7. X-squared = 0.4444, df = 1, p-value = 0.2525
  8. alternative hypothesis: true p is greater than 0.5
  9. 95 percent confidence interval:
  10.  0.3496555 1.0000000
  11. sample estimates:
  12.         p 
  13. 0.6666667 
  14. # 设置参数alternative="greater",使该检验变为单尾检验。由于p>0.05,故接受H0。

7、检验正态性

  1. > library(MASS)
  2. > attach(Cars93)
  3. > shapiro.test(Horsepower)      #W统计量检验
  5.     Shapiro-Wilk normality test
  7. data:  Horsepower
  8. W = 0.9358, p-value = 0.0001916
  10. #加载nortest包，运用更高端的统计检验来判定样本数据的正态性。
  11. > library(nortest)
  12. > ad.test(Horsepower)           #Anderson-Darling检验
  14.     Anderson-Darling normality test
  16. data:  Horsepower
  17. A = 1.2873, p-value = 0.002276
  19. > cvm.test(Horsepower)          #Carmer-Von Mises检验
  21.     Cramer-von Mises normality test
  23. data:  Horsepower
  24. W = 0.1605, p-value = 0.01689
  26. > lillie.test(Horsepower)       #Lillierfors检验
  28.     Lilliefors (Kolmogorov-Smirnov) normality test
  30. data:  Horsepower
  31. D = 0.1018, p-value = 0.01876
  33. > pearson.test(Horsepower)      #Pearson卡方检验正态性复合假设
  35.     Pearson chi-square normality test
  37. data:  Horsepower
  38. P = 28.4731, p-value = 0.001516
  40. > sf.test(Horsepower)           #Shapiro-Francia检验
  42.     Shapiro-Francia normality test
  44. data:  Horsepower
  45. W = 0.9372, p-value = 0.0004632

8、检验两个未知分布的总体，是否有类似的形状

  1. #这部分的检验，可用作Kruskal-Wallis检验的准备工作。
  2. > lens1 <- rweibull(100,shape=3,scale=4)
  3. > lens2 <- rt(200,df=1)
  4. > wilcox.test(lens1,lens2)
  6.     Wilcoxon rank sum test with continuity correction
  8. data:  lens1 and lens2
  9. W = 17512, p-value < 2.2e-16
  10. alternative hypothesis: true location shift is not equal to 0
  11. # p<0.05，表明两个总体不具有类似的形状。

9、检验相关系数的显著性
当我们把数据输入到计算机，计算出相关系数。此时应该保持清醒的头脑，问问自己：数据是否足够？相关系数足够大吗？此时应该用R运行cor.test()

  1. #xk和yk都取自正态总体。
  2. > cor(xk,yk)
  3. [1] 0.8352458
  4. > cor.test(xk,yk,method="Pearson")  
  5. #method参数可以选择Pearson、Spearman、Kendall三种方法。
  6.          Pearson's product-moment correlation
  8. data:  xk and yk
  9. t = 2.1481, df = 2, p-value = 0.1648
  10. alternative hypothesis: true correlation is not equal to 0
  11. 95 percent confidence interval:
  12.  -0.6491634  0.9857821
  13. sample estimates: 
  14.       cor 
  15. 0.8352458
  17. 如果在调用cor(xk,yk)后，就盲目接受结果，是愚蠢的。从cor.test()的结果可以看出：
  18. 不仅P>0.05（应接受H0），而且置信区间中也包含了0,所以不能确信0.8352458这个相关系数是显著的。

来自psych包的更高大上的相关性检验方法

  1. > library(MASS)
  2. > library(psych)
  3. > attach(Cars93)
  4. > car <- Cars93[,c("Weight","Horsepower","RPM","EngineSize")]
  5. > corr.test(car,method="kendall")
  6. Call:corr.test(x = car, method = "kendall")
  7. Correlation matrix     #相关系数矩阵
  8.            Weight     Horsepower     RPM     EngineSize
  9. Weight       1.00       0.63        -0.31       0.74
  10. Horsepower   0.63       1.00        -0.05       0.65
  11. RPM         -0.31      -0.05         1.00      -0.40
  12. EngineSize   0.74       0.65        -0.40       1.00
  13. Sample Size 
  14. [1] 93
  15. Probability values (Entries above the diagonal are adjusted for multiple tests.) 
  16.            Weight     Horsepower     RPM     EngineSize
  17. Weight          0       0.00         0.01       0
  18. Horsepower      0       0.00         0.63       0
  19. RPM             0       0.63         0.00       0
  20. EngineSize      0       0.00         0.00       0
  22. #从结果看出，RPM与Horsepower的相关系数为-0.31并不显著地不为0。

10、检验两个样本是否来自同一个分布
使用Kolmogorov-Smirnov方法进行检验，该方法适用于所有分布。

  1. > lens1 <- rweibull(500,shape=3,scale=4)
  2. > lens2 <- rt(300,df=1)
  3. > lens3 <- rweibull(10,shape=3,scale=4)
  4. > ks.test(lens1,lens2)
  6.     Two-sample Kolmogorov-Smirnov test
  8. data:  lens1 and lens2
  9. D = 0.794, p-value < 2.2e-16
  10. alternative hypothesis: two-sided
  11. # P<0.05,表明来自不同的分布。
  13. > ks.test(lens1,lens3)
  15.     Two-sample Kolmogorov-Smirnov test
  17. data:  lens1 and lens3
  18. D = 0.336, p-value = 0.1747
  19. alternative hypothesis: two-sided
  20. # P>0.05,表明来自同一分布。
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注