R语言怎么检验分布是不是T分布

Python011

R语言怎么检验分布是不是T分布,第1张

ks.test()实现了KS检验,可以检验任意样本是不是来自给定的连续分布。

你这里的用法就是:

ks.test(data,pt,df=df) #data是样本的数据,df是要检验的t分布的自由度

我们可以用很多方法分析一个单变量数据集的分布。最简单的办法就是直接看数

字。利用函数summary 和fivenum 会得到两个稍稍有点差异的汇总信息。此外,stem

(\茎叶"图)也会反映整个数据集的数字信息。

>attach(faithful)

>summary(eruptions)

Min. 1st Qu. Median Mean 3rd Qu. Max.

1.600 2.163 4.000 3.488 4.454 5.100

>fivenum(eruptions)

[1] 1.6000 2.1585 4.0000 4.4585 5.1000

>stem(eruptions)

The decimal point is 1 digit(s) to the left of the |

16 | 070355555588

18 | 000022233333335577777777888822335777888

20 | 00002223378800035778

22 | 0002335578023578

24 | 00228

26 | 23

28 | 080

30 | 7

32 | 2337

34 | 250077

36 | 0000823577

38 | 2333335582225577

40 | 0000003357788888002233555577778

42 | 03335555778800233333555577778

44 | 02222335557780000000023333357778888

46 | 0000233357700000023578

48 | 00000022335800333

50 | 0370

茎叶图和柱状图相似,R 用函数hist 绘制柱状图。

>hist(eruptions)

>## 让箱距缩小,绘制密度图

>hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)

>lines(density(eruptions, bw=0.1))

>rug(eruptions) # 显示实际的数据点

更为精致的密度图是用函数density 绘制的。在这个例子中,我们加了一条

由density 产生的曲线。你可以用试错法(trial-and-error)选择带宽bw(bandwidth)

因为默认的带宽值让密度曲线过于平滑(这样做常常会让你得到非常有\意思"的密度

分布)。(现在已经有一些自动的带宽挑选方法2,在这个例子中bw = "SJ"给出的结

果不错。)

我们可以用函数ecdf 绘制一个数据集的经验累积分布(empirical cumulative

distribution)函数。

>plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)

显然,这个分布和其他标准分布差异很大。那么右边的情况怎么样呢,就是火山

爆发3分钟后的状况?我们可以拟合一个正态分布,并且重叠前面得到的经验累积密

度分布。

>long <- eruptions[eruptions >3]

>plot(ecdf(long), do.points=FALSE, verticals=TRUE)

>x <- seq(3, 5.4, 0.01)

>lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)

分位比较图(Quantile-quantile (Q-Q) plot)便于我们更细致地研究二者的吻合

程度。

par(pty="s") # 设置一个方形的图形区域

qqnorm(long)qqline(long)

上述命令得到的QQ图表明二者还是比较吻合的,但右侧尾部偏离期望的正态分布。

我们可以用t 分布获得一些模拟数据以重复上面的过程

x <- rt(250, df = 5)

qqnorm(x)qqline(x)

这里得到的QQ图常常会出现偏离正态期望的长尾区域(如果是随机样本)。我们可以用

下面的命令针对特定的分布绘制Q-Q图

qqplot(qt(ppoints(250), df = 5), x, xlab = "Q-Q plot for t dsn")

qqline(x)

最后,我们可能需要一个比较正规的正态性检验方法。R提供了Shapiro-Wilk 检

>shapiro.test(long)

Shapiro-Wilk normality test

data: long

W = 0.9793, p-value = 0.01052

和Kolmogorov-Smirnov 检验

>ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))

One-sample Kolmogorov-Smirnov test

data: long

D = 0.0661, p-value = 0.4284

alternative hypothesis: two.sided

(注意一般的统计分布理论(distribution theory)在这里可能无效,因为我们用同样

的样本对正态分布的参数进行估计的。)

转载于:

http://www.biostatistic.net/thread-2413-1-1.html

R语言中一组数据服从威布尔分布。

首先可以利用fitdistr函数求得weibull分布的形状参数和尺寸参数,假设数据为x:library(MASS) #fitdistr需要利用MASS包fitdistr(x, densfun = "weibull",lower=0)得到形状参数shape与尺度参数scale

然后利用ks.test进行检验:ks.test(jitter(x),"pweibull",shape,scale)

上边的jitter用来做小扰动,因为如果x中有重复数据的话ks.test会报错,如果x中没有重复数据则不需要jitter。shape是得到的形状参数,scale是得到的尺度参数。

ks.test得到两个结果,一个是D,越小越好,一个是p-value,这个值要大于0.05

1. 判断存在:一个元素是不是在向量中用 a%in%b

>a="TT"

>b=c("AA","AT","TT")

>a %in% b

[1] TRUE

2. 判断某一元素这向量中的索引(第几个位置): index.TT=which(b==”TT”)

>index.TT=which(b=="TT")#index.TT是想知道的索引号,which是判断函数,b是想知道的元素所在的向量

>index.TT

[1] 3

3. 相当于 python 中的字典, names 函数

>b

[1] "AA" "AT" "TT"

>names(b)=c("geno1","geno2","geno3")#geno mean genotype

>names(b)

[1] "geno1" "geno2" "geno3"

>names(b)[1]

[1] "geno1"

>names(b)[1]="test"

>names(b)

[1] "test""geno2" "geno3"

>names(b)=NULL

>b

[1] "AA" "AT"

>b["geno2"]

"AT"

pop_name=c(“CEU”,"YRI")

names(pop_name)=c(1,2)

names(pop_name[1])=1

4. 去除某一元素: b[-index.nu]

#想去除元素”TT”,如果你不知道是第几个索引,可以先判断索引,再删除。

>b=c("AA","AT","TT")

>names(b)=c("geno1","geno2","geno3")

>index.TT=which(b=="TT")

>b=b[-index.TT]

>b

geno1 geno2

"AA""AT"

5. 相当于 Python 中的 set() 函数 和 count() 函数: unique() , table()

>b=c("TT","AT","AT","TT","AA")

>unique(b)#即相当于去除所有的重复,只保留一个

[1] "TT" "AT" "AA"

>table(b)#以元素为name,统计各元素的个数

b

AA AT TT

122

6. 字符串的分割: strsplit()

>test="AA"

>strsplit(test)

错误于strsplit(test) :缺少参数"split",也没有缺省值

>strsplit(test,split='')

[[1]]

[1] "A" "A"

>test=strsplit(test,split='')[[1]]

>test

[1] "A" "A"

7. 文本文档的写入: write.table()

write.table( res.matrix,file=new.file,sep='\t',quote=F,row.names=F,col.names=F,append=T)#quote=F去掉引号后写入,row.names=F去掉行的名字写入,否则会把名字写进去

##写入数据时候最好把数据存储成一个matrix然后直接写。要是每行每行写的话要注意数据的格式了。先建立一个空的matrix,见8,然后通过rbind或者cbind叠加上去。

方法一:

a=c()

b=c(“AA”,”TT”,”CC”)

for (i in 1:3){

a=c(a,b)

}

write.table(a,file=”test.txt”)#你会发现结果是

AA

TT

CC

….

##而且还有行和列的名字,因为没有设置参数。因为对于c向量来说,写的话默认是竖着写的,每个元素占一行。所以比较方便的就是rbind

方法二:

a=c()

b=c(“AA”,”TT”,”CC”)

for (i in 1:3){

a=rbind(a,b)

}

write.table(a,file=”test.txt”,quote=F,row.names=F,col.names=F)#你会发现结果是

AA TT CC

AA TT CC

AA TT CC

##原因是rbind把最总结果当做矩阵了。对于R数据的写入最好能生成最后的矩阵再写入。但是西面的梅一行写一次和方法二的效果是想通的,但是要用到append参数。

a=c()

b=c(“AA”,”TT”,”CC”)

for (i in 1:3){

a=rbind(a,b)

write.table(a,file=”test.txt”,quote=F,row.names=F,col.names=F,append=T)

}

8. 建立一个空的 matrix :

res.matrix <- matrix( ,nrow=0,ncol=6 )##这样就建立了一个0行6列的空matrix了。

9. 如何将 R 运行结果输出到文件

>x=read.table("F:/my/work/chengxu/PValue/pc2jieguo/pc2302.txt")

>z=t(x)

>ks.test(y,z)

Two-sample Kolmogorov-Smirnov test

data:y and z

D = 0.207, p-value <2.2e-16

alternative hypothesis: two-sided

如上面运行结果,我想将p-value <2.2e-16自动保存到一个文件中,如何用R程序实现,谢谢!

sink("output.txt")

print(ks.test(y,z)$p.value)

sink()

http://cos.name/cn/topic/16300

10 降序排列:

>a=c(1,1.2,0.1,4,5,-0.1)

>a=sort(a,decreasing=T)

>a

[1]5.04.01.21.00.1 -0.1

11. 取前1%的数

>a=c(1:10,4:20,1:100,1:1000)

>a=sort(a,decreasing=T)#先降序

>sig=a[round(length(a)*0.01)]

>sig

[1] 990

12.在shell中直接执行R脚本

R CMD BATCH --argstest.R

13. R中高级作图的方法

http://qizhi502.blog.163.com/blog/static/11497002520120611451736/

14:设置字体类型:

par(family='Times New Roman')

15:控制图形四周的空白大小

par(mfrow=c(3,1),mar=c(0,0,0,0))

其中mar是四周的间距,分别为x,y上下的距离

16控制作图区域的大小layout

layout(c(1,2,3),height=c(1,1,0.5))

分成竖着三份, 其中三份比列依次为(高度依次为2:2:1)

17保留两位小数

round(0.123,digits=2)

18 在原有图的基础上画图:

par(fig=c(0.1,0.5,0.43,0.65), new=TRUE)

19 只显示y轴

plot(1:10,1:10,axes=F)

axis(2,at.....)

20 调节刻度方向 las

plot(1:10,1:10,las=1)

21 屏幕分割

layout(matrix(1:16,4,4))###竖着plot

par(mfrow=c(4,4))##横着plot

22.逻辑表示或者

xor为异或,两值不等为真,两值相等为假。例:xor(0, 1)

23. 从向量中随机取几个数sample

sample(rep(1:1000),10)

23 字符串转换成小数浮点型

as.numeric("0.123")

24. 读取不规范的文本

f=readLines(afile,n=1)#n表示读几行

f=strsplit(f,'\t')##分割

f[1][[1]]##第一行

f[1][[1]][1]##第一行 第一个字符串

25. write 写入文件

write(afile, "a\tb\t",append=T) #沿着每行一次 写入

26. 不需要循环,这直接对matrix没行或者每列进行筛选操作apply()

apply(data,col2 or row1, max>0)

27.保留2位小数

a=2.300

a=as.numeric(sprintf(“%.3f”,a))

28。调出假设检验的p value

t.test(data1,data2)$p.value