交通数据分析基础 实验指导3 随机数

随机数实际分析中常用到概率函数，R语言中概率函数的形式为：[dpqr]函数缩写()。其中[dpqr]表示从中取一个字母。d：密度函数p：分布函数q：分位数函数r：生成随机数函数分布名称函数分布名称函数Beta分布betaLogistic分布logis二项分布binom多项分布multinom柯西分布cauchy负二项分布nbinom卡方分布chisq正态分布norm指数分布exp泊松分布poisF分布fWilcoxon符号秩分布signranGamma分布gammat分布t几何分布geom均匀分布unif超几何分布hyperWeibull分布weibull对数正态分布lnormWilcoxon秩和分布wilcox#生成10个服从标准正态分布随机数

rnorm(10,mean=0,sd=1)##[1]0.6490.612-0.452-0.2590.8610.6330.0310.7581.825-0.225#生成10个服从均匀分布的随机数

runif(10)##[1]0.8340.3490.5310.8830.5730.5540.8330.7640.9620.164#设置随机数种子

#R生成的是伪随机数，因此设置随机数种子后，每次生成的随机数都一样

set.seed(123)

runif(10)##[1]0.28760.78830.40900.88300.94050.04560.52810.89240.55140.4566#分位数

qnorm(c(0.025,0.5,0.975))##[1]-1.960.001.96#分布函数

pnorm(c(-1.96,0,1.96))##[1]0.0250.5000.975数据分布的QQ图检验生成数据框set.seed(1)

data<-data.frame(

value1=rnorm(300,mean=1,sd=2),#正态分布

value2=runif(300,min=0,max=10)#均匀分布

)#检验value1是否服从正态分布

library(ggplot2)

ggplot(data,aes(sample=value1))+

geom_qq()+

geom_qq_line()#检验value2是否服从正态分布

#从图中可以看出不符合正态分布

library(ggplot2)

ggplot(data,aes(sample=value2))+

geom_qq()+

geom_qq_line()#观察value2的密度曲线图，可以看出可能服从均匀分布

ggplot(data,aes(x=value2))+

geom_density()#观察value2的密度曲线图，可以看出可能服从均匀分布

value2_min<-min(data$value2)

value2_max<-max(data$value2)

ggplot(data,aes(sample=value2))+

geom_qq(distribution=stats::qunif,dparams=list(min=value2_min,max=value2_max))+

geom_qq_line(distribution=stats::qunif,dparams=list(min=value2_min,max=value2_max))区间估计#生成30个均值为50、标准差为5的数据

set.seed(123)

data<-rnorm(n=30,mean=50,sd=5)#均值μ的置信区间估计

#方差σ已知

#自定义函数计算Z分布置信区间

z_test_ci<-function(data,sigma,conf.level=0.95){

n<-length(data)

x_bar<-mean(data)

alpha<-1-conf.level

z<-qnorm(1-alpha/2)

margin<-z*sigma/sqrt(n)

ci<-c(x_bar-margin,x_bar+margin)

return(ci)

}

#调用函数

sigma_known<-5

ci_mu_known<-z_test_ci(data,sigma=sigma_known,conf.level=0.95)

print(paste("均值置信区间（方差已知）:","(",ci_mu_known[1],",",ci_mu_known[2],")"))##[1]"均值置信区间（方差已知）:(47.9752870761238,51.5536753635581)"#均值μ的置信区间估计

#方差σ未知

#直接使用内置函数t.test()：

t_result<-t.test(data,conf.level=0.95)

ci_mu_unknown<-t_result$

print(paste("均值置信区间（方差未知）:","(",ci_mu_unknown[1],",",ci_mu_unknown[2],")"))##[1]"均值置信区间（方差未知）:(47.9328667909283,51.5960956487537)"#方差σ²的置信区间估计

var_ci<-function(data,conf.level=0.95){

n<-length(data)

s2<-var(data)

alpha<-1-conf.level

chi_low<-qchisq(1-alpha/2,df=n-1)

chi_high<-qchisq(alpha/2,df=n-1)

ci<-c((n-1)*s2/chi_low,(n-1)*s2/chi_high)

return(ci)

}

#调用函数

ci_var<-var_ci(data,conf.level=0.95)

print(paste("方差置信区间:","(",ci_var[1],",",ci_var[2],")"))##[1]"方差置信区间:(15.2607287770605,43.4817900717857)"假设检验单正态总体set.seed(1)

sample_data<-rnorm(100,mean=5,sd=3)sample_mean<-mean(sample_data)

sample_sd<-sd(sample_data)

cat("样本均值:",sample_mean,"\n样本标准差:",sample_sd)##样本均值:5.33

##样本标准差:2.69#方差已知（σ=3）

#Z检验，检验均值是否为5

z_score<-(sample_mean-5)/(3/sqrt(100))

p_value_z<-2*pnorm(-abs(z_score))#双侧检验

cat("Z值:",z_score,"\np值:",p_value_z)##Z值:1.09

##p值:0.276#方差已知（σ=3）

#Z检验，检验均值是否为7

z_score<-(sample_mean-7)/(3/sqrt(100))

p_value_z<-2*pnorm(-abs(z_score))#双侧检验

cat("Z值:",z_score,"\np值:",p_value_z)##Z值:-5.58

##p值:2.44e-08#方差未知

#T检验，检验均值是否为5

t_score<-(sample_mean-5)/(sample_sd/sqrt(30))

p_value_t<-2*pt(-abs(t_score),df=29)#双侧检验

cat("T值:",t_score,"\np值:",p_value_t)##T值:0.664

##p值:0.512双正态总体假设检验#生成数据

set.seed(1)

n1<-30

n2<-40

mu1<-5

mu2<-7

sd1<-3

sd2<-4

sample1<-rnorm(n1,mean=mu1,sd=sd1)

sample2<-rnorm(n2,mean=mu2,sd=sd2)#假设方差未知

#假设两个样本的总体方差相等

t_test_equal<-t.test(sample1,sample2,var.equal=TRUE)

cat("方差相等假设下的T检验结果:\n")##方差相等假设下的T检验结果:print(t_test_equal)##

##TwoSamplet-test

##

##data:sample1andsample2

##t=-3,df=68,p-value=0.002

##alternativehypothesis:truedifferenceinmeansisnotequalto0

##95percentconfidenceinterval:

##-4.222-0.983

##sampleestimates:

##meanofxmeanofy

##5.257.85#假设方差未知

#假设两个样本的总体方差不相等

t_test_unequal<-t.test(sample1,sample2,var.equal=FALSE)

cat("方差不相等假设下的T检验结果:\n")##方差不相等假设下的T检验结果:print(t_test_unequal)##

##WelchTwoSamplet-test

##

##data:sample1andsample2

##t=-3,df=68,p-value=0.001

##alternativehypothesis:truedifferenceinmeansisnotequalto0

##95percentconfidenceinterval:

##-4.15-1.05

##sampleestimates:

##meanofxmeanofy

##5.257.85#假设方差已知

mean1<-mean(sample1)

mean2<-mean(sample2)

z_value<-(mean1-mean2)/sqrt(sd1^2/n1+sd2^2/n2)

p_value<-2*pnorm(-abs(z_value))

cat("Z检验统计量:",z_value,"\n")##Z检验统计量:-3.11cat("Z检验的p值:",p_value,"\n")##Z检验的p值:0.00187方差分析单因素方差分析set.seed(1)

#创建数据框

group<-factor(rep(c("Group1","Group2","Group3"),each=30))

value<-c(rnorm(30,mean=5,sd=2),

rnorm(30,mean=7,sd=2),

rnorm(30,mean=6,sd=2))

data<-data.frame(group,value)

#查看数据结构

str(data)##'data.frame':90obs.of2variables:

##$group:Factorw/3levels"Group1","Group2",..:1111111111...

##$value:num3.755.373.338.195.66...#使用aov()函数拟合ANOVA模型

model_anova<-aov(value~group,data=data)

#查看模型摘要

summary(model_anova)##DfSumSqMeanSqFvaluePr(>F)

##group266.233.110.39.6e-05***

##Residuals87279.43.2

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1多因素方差分析#生成示例数据

set.seed(1)

A<-factor(rep(c("A1","A2"),each=30))#因素A（2水平）

B<-factor(rep(c("B1","B2","B3"),each=10))#因素B（3水平）

y<-c(rnorm(10,mean=5),rnorm(10,mean=6),rnorm(10,mean=7),

rnorm(10,mean=5.5),rnorm(10,mean=6.5),rnorm(10,mean=7.5))

#创建数据框

data<-data.frame(A,B,y)

#查看数据结构

str(data)##'data.frame':60obs.of3variables:

##$A:Factorw/2levels"A1","A2":1111111111...

##$B:Factorw/3levels"B1","B2","B3":1111111111...

##$y:num4.375.184.166.65.33...#执行无交互作用的双因素ANOVA

model_manova1<-aov(y~A+B,data=data)

#查看结果

summary(model_manova1)##DfSumSqMeanSqFvaluePr(>F)

##A14.54.545.950.018*

##B235.517.7523.254.5e-08***

##Residuals5642.80.76

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1#执行有交互作用的双因素ANOVA

model_manova2<-aov(y~A*B,data=data)

#查看结果

summary(model_manova2)##DfSumSqMeanSqFvaluePr(>F)

##A14.54.545.790.02*

##B235.517.7522.647.2e-08***

##A:B20.40.210.260.77

##Residuals5442.30.78

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1线性回归set.seed(1)

n<-100#样本量

#生成连续变量

x1<-rnorm(n,mean=50,sd=10)

x2<-runif(n,min=0,max=10)

#生成分类变量（3类）

category<-factor(sample(c("A","B","C"),size=n,replace=TRUE,

prob=c(0.4,0.35,0.25)))

#生成响应变量（包含交互效应和异方差）

y<-2+1.5*x1-0.8*x2+

as.numeric(category=="B")*3+

as.numeric(category=="C")*5+

rnorm(n,sd=x2)#异方差设置

#创建数据框

data<-data.frame(y,x1,x2,category)

head(data)##yx1x2category

##169.643.72.68B

##281.751.82.19A

##371.541.65.17B

##4100.966.02.69B

##576.453.31.81A

##678.541.85.19C#不包含虚拟变量的回归

fit<-lm(y~x1+x2,data=data)

summary(fit)##

##Call:

##lm(formula=y~x1+x2,data=data)

##

##Residuals:

##Min1QMedian3QMax

##-16.078-2.5060.1173.24715.239

##

##Coefficients:

##EstimateStd.ErrortvaluePr(>|t|)

##(Intercept)2.56133.06920.830.4061

##x11.52230.058426.09<2e-16***

##x2-0.56500.1880-3.000.0034**

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1

##

##Residualstandarderror:5.2on97degreesoffreedom

##MultipleR-squared:0.875,AdjustedR-squared:0.873

##F-statistic:341on2and97DF,p-value:<2e-16#包含交互项的回归

fit<-lm(y~x1*x2,data=data)

summary(fit)##

##Call:

##lm(formula=y~x1*x2,data=data)

##

##Residuals:

##Min1QMedian3QMax

##-16.906-2.7690.2813.30215.457

##

##Coefficients:

##EstimateStd.ErrortvaluePr(>|t|)

##(Intercept)7.38455.36101.380.17

##x11.42800.103913.75<2e-16***

##x2-1.67951.0334-1.630.11

##x1:x20.02160.01971.100.28

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1

##

##Residualstandarderror:5.19on96degreesoffreedom

##MultipleR-squared:0.877,AdjustedR-squared:0.873

##F-statistic:228on3and96DF,p-value:<2e-16#包含虚拟变量的回归

fit<-lm(y~x1+x2+category,data=data)

summary(fit)##

##Call:

##lm(formula=y~x1+x2+category,data=data)

##

##Residuals:

##Min1QMedian3QMax

##-14.020-2.6300.1282.09512.554

##

##Coefficients:

##EstimateStd.ErrortvaluePr(>|t|)

##(Intercept)0.23162.80650.080.93442

##x11.52440.052629.00<2e-16***

##x2-0.66290.1715-3.860.00020***

##categoryB3.82161.12263.400.00097***

##categoryC5.43801.15724.708.8e-06***

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1

##

##Residualstandarderror:4.68on95degreesoffreedom

##MultipleR-squared:0.901,AdjustedR-squared:0.897

##F-statistic:216on4and95DF,p-value:<2e-16#包含虚拟变量、虚拟变量与一般变量的交互项的回归

fit<-lm(y~x1+x2+category+x1:category,data=data)

summary(fit)##

##Call:

##lm(formula=y~x1+x2+category+x1:category,data=data)

##

##Residuals:

##Min1QMedian3QMax

##-14.213-2.5950.0232.13412.144

##

##Coefficients:

##EstimateStd.ErrortvaluePr(>|t|)

##(Intercept)-1.38014.2348-0.330.74523

##x11.55580.081819.03<2e-16***

##x2-0.66230.1744-3.800.00026***

##categoryB5.85096.38300.920.36170

##categoryC9.25796.98771.320.18846

##x1:categoryB-0.03970.1229-0.320.74706

##x1:categoryC-0.07460.1345-0.550.58064

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1

##

##Residualstandarderror:4.72on93degreesoffreedom

##MultipleR-squared:0.901,AdjustedR-squared:0.895

##F-statistic:142on6and93DF,p-value:<2e-16#异方差检验

fit<-lm(y~x1+x2,data=data)

library(lmtest)

bp_test<-bptest(fit)

print(bp_test)#p<0.05表明是否存在异方差##

##studentizedBreusch-Pagantest

##

##data:fit

##BP=17,df=2,p-value=2e-04#异方差稳健标准误

library(sandwich)

coeftest(fit,vcov=vcovHC(fit,type="HC3"))#HC3类型稳健标准误##

##ttestofcoefficients:

##

##EstimateStd.ErrortvaluePr(>|t|)

##(Intercept)2.56133.47440.740.46

##x11.52230.067722.50<2e-16***

##x2-0.56500.2382-2.370.02*

##---

##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1#误差项独立性的Durbin-Watson检验

f