钱老师在书中的这一章介绍了R语言的数据类型以及语法结构特别推荐了R Commander给初学者使用。这些内容本笔记略过。主要讲述本章中提到的一个案例：【如果水体检测值有10%超过标准就被列为“受损的”水体】这个10%是否合理。
要检验这个标准是否合理，最简单的办法就是分批次地在已达标额水体中不断重复采样并测定TP，然后根据这个规定来确定我们把水体判断为“受损”水体的几率，或者反过来采受损水体。显然自然水体的变化多样，这样的操作实际上并不可行。而，如果我们知道水体的TP浓度变量的分布，我们可以用计算机模拟实际的采样过程。采集水样和测定浓度可以从随机分布中抽取一个随机数来模拟。模拟是评估随机变量行为的一种统计学工具。
由于大多数水质浓度变量都可以看做服从对数正态分布，我们可以用浓度的对数值并假设满足正态分布。
在我们的案例中假设水质标准的对数值是3 ，并且污染水体水质浓度对数的分布为。 水质不达标的时间超过10%，该水体被认为受损水体。我们可以估算出水体在环保局10%的规则下被误认为是达标水体的概率。假设我们有10个样本，未受损意味着只有一个或没有观测值是超过3的。对每次测量浓度小于3是达标的，对某一水体10个样本中总超标的样本个数小于2（即，0和1）是达标的。
Norm1 <- as.data.frame(matrix(rnorm(10000*10, mean=2, sd=1), ncol=10))
rownames(Norm1) <- paste("sample", 1:10000, sep="")
colnames(Norm1) <- paste("obs", 1:10, sep="")
Norm1$violations <- with(Norm1,
(obs1>3)+(obs2>3)+(obs3>3)+(obs4>3)+(obs5>3)+
(obs6>3)+(obs7>3)+(obs8>3)+(obs9>3)+(obs10>3))
Norm1$impaired <- with(Norm1, as.numeric(violations<2))
numSummary(Norm1[,"impaired", drop=FALSE], statistics=c("mean", "sd", "IQR",
"quantiles"))
mean sd IQR 0% 25% 50% 75% 100% n
0.5136 0.49984 1 0 0 1 1 1 10000
可以看出，我们用10000和水体，每个水体测10个样本，把污染水体判断为达标的概率是0.51。也就是将近一半的都判断错了。
通常，犯错的概率高要归结于观测值的数量少。我们可以用同样的模拟方法来看如果每个水体样本量增大到100的话，犯错的概率有没有降低。那么如果超过3的样本数小于10的话就认为是达标的。
Norm2 <- as.data.frame(matrix(rnorm(10000*100, mean=2, sd=1), ncol=100))
rownames(Norm2) <- paste("sample", 1:10000, sep="")
colnames(Norm2) <- paste("obs", 1:100, sep="")
Norm2$violations <- apply(X=Norm2,MARGIN = 1,FUN = function(x){
return(sum(x>3))
})
Norm2$impaired <- with(Norm2, as.numeric(violations<=10))
numSummary(Norm2[,"impaired", drop=FALSE], statistics=c("mean", "sd", "IQR",
"quantiles"))
mean sd IQR 0% 25% 50% 75% 100% n
0.064 0.2447652 0 0 0 0 0 1 10000
当样本量增加时把污染样本判断为达标样本的概率减小了。但是我们看看把达标样本判断为污染样本的概率。假设达标水体水质浓度对数的分布为
Norm1 <- as.data.frame(matrix(rnorm(10000*10, mean=2, sd=0.75), ncol=10))
rownames(Norm1) <- paste("sample", 1:10000, sep="")
colnames(Norm1) <- paste("obs", 1:10, sep="")
Norm1$violations <- with(Norm1,
(obs1>3)+(obs2>3)+(obs3>3)+(obs4>3)+(obs5>3)+
(obs6>3)+(obs7>3)+(obs8>3)+(obs9>3)+(obs10>3))
Norm1$impaired <- with(Norm1, as.numeric(violations>2))
numSummary(Norm1[,"impaired", drop=FALSE], statistics=c("mean", "sd", "IQR",
"quantiles"))
mean sd IQR 0% 25% 50% 75% 100% n
0.0536 0.2252379 0 0 0 0 0 1 10000
Norm2 <- as.data.frame(matrix(rnorm(10000*100, mean=2, sd=0.75), ncol=100))
rownames(Norm2) <- paste("sample", 1:10000, sep="")
colnames(Norm2) <- paste("obs", 1:100, sep="")
Norm2$violations <- apply(X=Norm2,MARGIN = 1,FUN = function(x){
return(sum(x>3))
})
Norm2$impaired <- with(Norm2, as.numeric(violations>10))
numSummary(Norm2[,"impaired", drop=FALSE], statistics=c("mean", "sd", "IQR",
"quantiles"))
mean sd IQR 0% 25% 50% 75% 100% n
0.3061 0.4608948 1 0 0 0 1 1 10000
虽然当样本量增加时把污染样本判断为达标样本的概率减小了，但是把达标样本判断为污染的概率却增加了。
网友评论