2020-12-15miki笔记

R数据科学###

20200919 8:30 raining####

shif+enter 换行

planes %>%
count(tailnum)%>%
filter(n>1)

weather%>%
count(year,month,day,hour,origin)%>%

后面的年月日等条件捆绑计数

filter(n>1)

flights2%>%
select(-origin,-dest)%>%
left_join(airlines,by="carrier")

left_join左连接

x%>%
right_join(y,by="key")

right_join右连接

x%>%
inner_join(y,by="key")

inner_join内连接取交集

x%>%
full_join(y,by="key")

full_join全链接取并集

semi_join(x,y,by="key")

semi_join保留x表中与Y表中的观测相匹配的所有x

anti_join(x,y,by="key")

anti_join去除x表中与Y表中的观测相匹配的所有x

top_dest <-flights%>%
count(dest,sort=TRUE)%>%

sort=TRUE降序

head(10)

head()前几位

top_dest

优化的方法optim()

best<-optim(c(0,0),measure_distance,data=sim1)
best$par

sim1_mod<-lm(y~x,data=sim1)

lm()线性模型

coef(sim1_mod)

coef()调出

线性回归

普通最小二乘

20200920 8:30 经管院b224 raining####

y是均值为0，标准差为2计算10000次的正态分布，x是男女随机取10000次，y=b+ax

y<-rnorm(10000, mean=0, sd=2)
x<-sample(
0:1,10000,replace=TRUE
)
mm<-lm(y~x)

mm<-lm(y~x-1) 减掉1表示去除截距

coef(mm)

textA=scan("D:/mikistudysoftware/RStudio/text A.txt",what="",sep="",na.strings="NA")

导入文本文件，注意文本路径""改成"/"

what:声明读入为字符类型数据，可能指定读入的精度/类型

sep：指定各个读入的数据之间的分隔符；默认情况下分隔符：空格、tab；如果不是其它分隔符，例如“：/”通过SEP来指定

可以通过list指定读入变量的变量名，同时生成的对象为列表，则可以同时读入字符与数字

Skip 从第几行开始读入数据

Nlines 指定最大读入行数

如果通过键盘输入的时候，不typeof(hw1_a$ID)希望出现下标提示，则可以使用：quiet=TRUE

encoding =””指定的编码格式，有时候读入的中文可能会出现乱码的时候，可能通过这个参数来指定：Latin-1 或者 UTF-8

juhaoA=sum(str_detect(textA,"\."))

统计文本A中的.的个数

wordA=length(textA)

统计文本A文字的个数

willA=sum(str_detect(textA,"will"))

统计文本A中“will”这个词语的个数

20200927

箱线图适用于连续变量异常值检测

箱线图

ggplot(hw1_fj_dna,aes(,Income))+

geom_boxplot()

去除极端值的箱线图

boxplot(hw1_fj_dna$Income, outline = F)

boxplot(hw1_fj_dna$Income)

极端值

boxplot.stats(hw1_fj_dna $Income)$ out

去除极端值

hw1_fj_dna_IncomeFilter<-hw1_fj_dna $Income[!hw1_fj_dna$ Income%in% boxplot.stats(hw1_fj_dna $Income)$ out]
hw1_fj_dna_IncomeFilter

不使用科学计数法

options(scipen=3)

20200925 14:01 RHW1###

1. 请将数据hw1_a和hw1_b分别读入R，查看数据并指出各个变量的形式，最小值，最大值，中值，均值，标准差。

读取hw1_a.csv文件

install.packages("openxlsx")
library(openxlsx)
hw1_a<-read.csv("D:/mikiStudyFiles/les/hw1_a.csv")
hw1_a

方法1：陈肯-查看hw1_a各变量类型，最小值，最大值，中值，均值，标准差

str(hw1_a)

直接查看hw1_a中所有字段的类型

str(hw1_b)
summary(hw1_a)
summary(hw1_b)

直接汇总求出hw1_a中所有字段的max min median mean 1st Qu. 3rd Qu.

方法2：miki

typeof(hw1_a $ID) typeof(hw1_a$ Age)
typeof(hw1_a $Years_at_Employer) typeof(hw1_a$ Years_at_Address)
typeof(hw1_a $Income) max(hw1_a$ ID,na.rm=TRUE)
max(hw1_a $Age,na.rm=TRUE) max(hw1_a$ Years_at_Employer,na.rm=TRUE)
max(hw1_a $Years_at_Address,na.rm=TRUE) max(hw1_a$ Income,na.rm=TRUE)
min(hw1_a $ID,na.rm=TRUE) min(hw1_a$ Age,na.rm=TRUE)
min(hw1_a $Years_at_Employer,na.rm=TRUE) min(hw1_a$ Income,na.rm=TRUE)
min(hw1_a $Years_at_Address,na.rm=TRUE) median(hw1_a$ ID,na.rm=TRUE)
median(hw1_a $Age,na.rm=TRUE) median(hw1_a$ Years_at_Employer,na.rm=TRUE)
median(hw1_a $Income,na.rm=TRUE) median(hw1_a$ Years_at_Address,na.rm=TRUE)
mean(hw1_a $ID,na.rm=TRUE) mean(hw1_a$ Age,na.rm=TRUE)
mean(hw1_a $Years_at_Employer,na.rm=TRUE) mean(hw1_a$ Years_at_Address,na.rm=TRUE)
mean(hw1_a $Income,na.rm=TRUE) sd(hw1_a$ ID,na.rm=TRUE)
sd(hw1_a $Age,na.rm=TRUE) sd(hw1_a$ Years_at_Employer,na.rm=TRUE)
sd(hw1_a $Years_at_Address,na.rm=TRUE) sd(hw1_a$ Income,na.rm=TRUE)

读取hw1_b.csv文件

hw1_b<-read.csv("D:/mikiStudyFiles/les/hw1_b.csv")
hw1_b

查看hw1_b各变量类型，最小值，最大值，中值，均值，标准差

typeof(hw1_b $ID) typeof(hw1_b$ Credit_Card_Debt)
typeof(hw1_b $Automobile_Debt) typeof(hw1_b$ Is_Default)
max(hw1_b $ID) max(hw1_b$ Credit_Card_Debt)
max(hw1_b $Automobile_Debt) max(hw1_b$ Is_Default)
min(hw1_b $ID) min(hw1_b$ Credit_Card_Debt)
min(hw1_b $Automobile_Debt) min(hw1_b$ Is_Default)
median(hw1_b $ID) median(hw1_b$ Credit_Card_Debt)
median(hw1_b $Automobile_Debt) median(hw1_b$ Is_Default)
mean(hw1_b $ID) mean(hw1_b$ Credit_Card_Debt)
mean(hw1_b $Automobile_Debt) mean(hw1_b$ Is_Default)
sd(hw1_b $ID) sd(hw1_b$ Credit_Card_Debt)
sd(hw1_b $Automobile_Debt) sd(hw1_b$ Is_Default)

2. 结合上课我们所学的几种数据join 的形式，将两个数据集进行合并。对于每种数据合并的方式，请说明key, 并且报告合并后的数据样本总行数。

left_join(data,by="连接字段")，nrow() 表示行数

left_join

hw1_lf<-hw1_a%>%
left_join(hw1_b,by="ID")
hw1_lf

right_join

hw1_rj<-hw1_a%>%
right_join(hw1_b,by="ID")
hw1_rj
nrow(hw1_lf)
nrow(hw1_rj)

inner_join

hw1_ij<-hw1_a%>%
inner_join(hw1_b,by="ID")
hw1_ij
nrow(hw1_ij)

full_join

hw1_fj<-hw1_a%>%
full_join(hw1_b,by="ID")
hw1_fj
nrow(hw1_fj)

semi_join

hw1_sj<-hw1_a%>%
semi_join(hw1_b,by="ID")
hw1_sj
nrow(hw1_sj)

anti_join

hw1_aj<-hw1_a%>%
anti_join(hw1_b,by="ID")
hw1_aj
nrow(hw1_aj)

3. 请筛选出hw1_a 中收入大于4000的样本，并将此样本和hw1_b 中Is_Default=1的样本合并，你可以使用inner join的方式。这一问中你可以用pipe的书写形式。

比较的时候一定要用==，=表示赋值

hw1_filter<-hw1_a%>%
filter(hw1_a $Income>40000)%>% inner_join( hw1_b%>% filter(hw1_b$ Is_Default==1),by="ID"
)
hw1_filter

4. 在第2问的基础上, 请给出Income对Years_at_Employer的散点图，你发现了哪些趋势和现象?

ctr+shif+c多行注释

# #画图 ggplot(data,aes(x,y))+

# # geom_point()

na.rm=TRUE 排除空值

ggplot(hw1_fj,aes(Years_at_Employer,Income))+
geom_point(na.rm = TRUE)+
scale_y_continuous(labels = scales::comma)

5. 在第4问的基础上按照Is_Default 增加一个维度，请展示两变量在不同违约状态的散点图。请使用明暗程度作为区分方式

alpha表示明暗程度 as.factor()表示转换成向量

ggplot(hw1_fj,aes(Years_at_Employer,Income,alpha=as.factor(Is_Default)))+
geom_point(na.rm = TRUE)+
scale_y_continuous(labels = scales::comma)

6. 对于第5问，请使用形状作为另外一种区分方式。

shape表示形状

ggplot(hw1_fj,aes(Years_at_Employer,Income,shape=as.factor(Is_Default)))+
geom_point(na.rm = TRUE)+
scale_y_continuous(labels = scales::comma)

7. 请找出各个列的缺失值，并删除相应的行。请报告每一变量的缺失值个数，以及所有缺失值总数。

sum()表示计数，is.na()判断是否为空值

方法1：陈肯-筛选出hw1_fj第2列不为空的所有记录

hw1_fj=filter(hw1_fj,!is.na(hw1_fj[2]))

方法2：miki

计算缺失值总数

sum(is.na(hw1_fj))

计算每个变量的缺失值数

sum(is.na(hw1_fj $ID)) sum(is.na(hw1_fj$ Age))
sum(is.na(hw1_fj $Years_at_Employer)) sum(is.na(hw1_fj$ Years_at_Address))
sum(is.na(hw1_fj $Income)) sum(is.na(hw1_fj$ Automobile_Debt))
sum(is.na(hw1_fj $Credit_Card_Debt)) sum(is.na(hw1_fj$ Is_Default))

查看缺失值的行

hw1_fj[!complete.cases(hw1_lf),]

删除缺失值的行

hw1_fj_dna<-hw1_lf[complete.cases(hw1_lf),]

删除第25列有空值的行

manghe[complete.cases(manghe[,25]),]

8. 找出Income中的极端值并滤掉对应行的数据

quantile()表示取百分位比的值

quantile(hw1_a$Income,c(0.025,0.975))
hw1_a2=filter(hw1_a,Income>14168.81&Income<173030.92)

先做直方图

ggplot(hw1_fj_dna,aes(Income,))+
geom_histogram(bins = 30)+
scale_x_continuous(labels = scales::comma)

将频数在1-5的数据放大

ggplot(hw1_fj_dna,aes(Income))+
geom_histogram(bins = 30)+
coord_cartesian(ylim = c(1,5))+
scale_x_continuous(labels = scales::comma)

通过图形筛选出Income中的极端值

(hw1_fj_dna%>%
filter(hw1_fj_dna $Income > 100000))$ Income

删除极端值行

hw1_fj_dna_IncomeFilter<-hw1_fj_dna[!hw1_fj_dna $Income%in% (hw1_fj_dna%>% filter(hw1_fj_dna$ Income > 100000))$Income,]
hw1_fj_dna_IncomeFilter
nrow(hw1_fj_dna_IncomeFilter)

箱线图连续型变量异常值判断##

ggplot(hw1_fj_dna,aes(,Income))+

geom_boxplot()

去除极端值的箱线图

boxplot(hw1_fj_dna$Income, outline = F)

boxplot(hw1_fj_dna$Income)

极端值

boxplot.stats(hw1_fj_dna $Income)$ out

去除极端值

hw1_fj_dna_IncomeFilter<-hw1_fj_dna[!hw1_fj_dna $Income%in% boxplot.stats(hw1_fj_dna$ Income)$out]

hw1_fj_dna_IncomeFilter

9. 将Income对数化，并画出直方图和density curve.

方法1：陈肯-密度直方图简便画法

inc<-hw1_a$Income
lninc<-log(inc)
hist(lninc,prob=T)
lines(density(lninc),col="blue")

方法2：miki

对数化Income直方图

ggplot(hw1_fj_dna_IncomeFilter,aes(Income,))+
geom_histogram(na.rm = TRUE,bins=30,color="gray")+
scale_x_log10(labels = scales::comma)

对数化Income密度曲线

ggplot(hw1_fj_dna_IncomeFilter,aes(Income,))+
geom_density(na.rm = TRUE,color="pink",size=1)+
scale_x_log10(labels = scales::comma)

二个图合在一个图里很奇怪密度曲线会很平缓不明显，这部分还需要琢磨琢磨

10. 以Income作为因变量，Years at Employer作为自变量，进行OLS回归，写出回归的方程，并指出自变量系数是否在某一显著性水平上显著。同时，解释你的结果(这一问你自己发挥可以找code解决)。

OLS回归方程lm(y~x,data) summary()计算统计量

y<-hw1_fj_dna_IncomeFilter $Income x<-hw1_fj_dna_IncomeFilter$ Years_at_Employer
miki<-lm(y~x,hw1_fj_dna_IncomeFilter)
miki

计算描述性统计量

summary(miki)

###################################20201007 RHW2####################################

1.编写函数get.root(a,b,c),求解一元二次方程ax^2+bx+c=0的实根;

创建函数function(),if(){}else if(){}else{}

rt<-function(a,b,c){
if(a==0&b==0&c==0){

a,b,c均为0,等式恒成立，有无穷解,无意义

root<-NA
}
else if(a==0&b==0){

a,b均为0，c不为0,无解

root<-NULL
}
else if(a==0){

a为0,b,c不为0，有一个根

root<- -c/b
} else{

a,b,c均不为0，判断delta

d<-b^2-4ac
if(d>=0){
root<-(-b+c(1,-1)sqrt(d))/(2a)
} else
(root<-NULL)
}
return(root)}

2.已知一元二次方程ax^2+bx+c=0的三个系数都是随机变量，其中a服从[1,5]上的均匀分布，b服从正态分布N(3,10),c服从均值为1的指数分布，请编写函数get.prop(),计算该方程有实根的概率；

rnorm(n, mean=0, sd=1) 高斯（正态）分布

rexp(n, rate=1) 指数分布

rgamma(n, shape, scale=1) γ分布

rpois(n, lambda) Poisson分布

rweibull(n, shape, scale=1) Weibull分布

rcauchy(n, location=0, scale=1) Cauchy分布

rbeta(n, shape1, shape2) β分布

rt(n, df) t分布

rf(n, df1, df2) F分布

rchisq(n, df) χ 2 分布

rbinom(n, size, prob)二项分布

rgeom(n, prob)几何分布

rhyper(nn, m, n, k) 超几何分布

rlogis(n, location=0, scale=1) logistic分布

rlnorm(n, meanlog=0, sdlog=1)对数正态

rnbinom(n, size, prob)负二项分布

runif(n, min=0, max=1)均匀分布

rwilcox(nn, m, n), rsignrank(nn, n) Wilcoxon分布

方法1：陈肯-使用循环计算

get.prob<-function(n){
a=runif(n,min=1,max=5)
b=rnorm(n,mean=3,sd=sqrt(10))
c=rexp(n,rate=1)
k=0
for (i in 1:n) {
if(sign(b[i]b[i]-4a[i]*c[i])==1|0)
{k=k+1}
}
return(k/n)
}

get.prob(100000)

方法2

创建函数function(),sum()计数注意：N(3,10) 表示期望是3，方差是10，标准差是根号10，所以表示为rnorm(n,3,sqrt(10))

p<-function(n){
a<-runif(n,1,5)
b<-rnorm(n,3,sqrt(10))
c<-rexp(n,1)
d<-b^2-4ac
sum(d>=0)/n
}

问题2：给定数据，请完成以下任务，请给出code 和输出结果。

1.请读入数据，使用软件分别给出 price, marketshare，和brand的缺失值数量。请按照每一个brand, 将数据按照先marjetshare 后price 进行从高到低排序；

读数read.csv(),sum(is.na())计数缺失值

导入csv数据，并统计缺失值数量

library(openxlsx)
library(tidyverse)
hw2<-read.csv("D:/mikiStudyFiles/les/hw2.csv")
hw2
sum(is.na(hw2 $price)) sum(is.na(hw2$ marketshare))
sum(is.na(hw2$brand))

group_by()分组,summarise()分组摘要,arrange(desc())降序

按照brand分组，并按marketshare,price降序排序

hw2_arr<-hw2[complete.cases(hw2),]%>%
group_by(brand)%>%
arrange(desc(marketshare,price))%>%
summarise(brand,marketshare,price)
hw2_arr

2.请按照brand 的种类，对price和marketshare 求均值。

hw2_brand_avg<-hw2[complete.cases(hw2),]%>%
group_by(brand)%>%
summarise(avg_price=mean(price,na.rm=TRUE),avg_marketshare=mean(marketshare,na.rm=TRUE))
hw2_brand_avg

3.请按照brand 的种类，对price和marketshare 画散点图。

方法1：陈肯-分几个小块

ggplot(data=hw2[complete.cases(hw2),])+
geom_point(mapping = aes(x=price,y=marketshare))+
facet_wrap(~brand,nrow=2)

方法2：miki -按brand分颜色

hw2_point<-ggplot(hw2[complete.cases(hw2),],aes(price,marketshare,color=as.factor(brand)))+
geom_point(na.rm=TRUE)
hw2_point

4.请按照价格的均值，产生新的变量price_new, 低于均值为“低价格”，高于均值为“高价格”。同样对市场份额也是，产生变量marketshare_new, 数值为“低市场份额”和“高市场份额”

新增变量mutate(),将某个变量分组ifelse(条件,true,false)

price_avg<-mean(hw2[complete.cases(hw2),]$price,na.rm=TRUE)
hw2_mutate<- hw2[complete.cases(hw2),] %>%
mutate(
price_new = ifelse(
price >price_avg,'高价格','低价格'
),
marketshare_new=ifelse(
price >price_avg,'高市场份额','低市场份额'
)
)
hw2_mutate

case_when的用法

hw2_mutate <- hw2 %>%

mutate(

price_new = case_when(

price >price_avg ~ '低',

TRUE ~ '高'

))

5.请估计模型，marketshare为Y，price为X.

一般线性模型lm()

y<-hw2 $marketshare x<-hw2$ price
miki<-lm(y~x,hw2)

coef()表示调用函数

coef(miki)

summary()表示计算描述性统计量

summary(miki)

6. 请画出（5）的拟合直线

方法1：陈肯-geom_abline()

ggplot(data=hw2)+
geom_point(aes(x=price,y=marketshare))+
geom_abline(data= m1,col= "blue")

方法2： miki- modelr::add_predictions()计算预测值

library(modelr)
hw2_mod<-hw2%>%
add_predictions(miki)
ggplot(hw2,aes(x),na.rm=TRUE)+
geom_point(aes(y=y),na.rm=TRUE)+
geom_line(aes(y=pred),na.rm=TRUE,
hw2_mod,
color="pink",
size=1)

7.请随机产生若干直线，验证（5）的结果是最优的

方法1：陈肯-使用最小二乘法验证拟合，拟合准则是使与的距离的平方和最小

b0=runif(20000,-5,5)
b1=runif(20000,-5,5)

d<-NA
sum<-NA
n<-1

while(n<=20000){
for(i in 1:24){
d[i]<-(marketshare[i]-b0[n]-b1[n]*price[i])^2}
sum[n]<-sum(d)
n<-n+1
}

resi=m1$residuals
resi2=sum(resi^2)

check=sum(as.numeric(sum<resi2))

方法2：mikii-使用原理取N条线，求最优解，计算点到直线的距离最小的线

添加斜线函数geom_abline(),随机取N个斜率和截距

models<-tibble(
k=runif(1000,-1,1),
b=runif(1000,-1,1)
)
ggplot(hw2,aes(x,y))+
geom_abline(aes(intercept=b,slope=k),models,alpha=1/4)+
geom_point(na.rm=TRUE)

optim()求最优解

a1<-runif(1000,-1,1)
a2<-runif(1000,-1,1)
model1 <- function(a, hw2) {
a[1] + x * a[2]
}
model1(c(7, 1.5), hw2)
measure_distance <- function(mod, hw2) {
diff <- y[!is.na(y)] - model1(mod,hw2)[!is.na(model1(mod,hw2))]
sqrt(mean(diff ^ 2))
}
hw2_dist <- function(a1, a2) {
measure_distance(c(a1, a2), hw2)
}
models <- models %>%
mutate(dist = purrr::map2_dbl(a1, a2, hw2_dist))
models
best<-optim(c(0, 0),measure_distance,hw2)
best$par

计算每一个点到直线的垂直距离

计算预测值

model1<-function(k,b,hw2){
b+x*k
}

计算均方根误差

measure_distance <- function(k1,b1,hw2) {
diff <- y[!is.na(y)] - model1(k1,b1,hw2)[!is.na(model1(k1,b1,hw2))]
sqrt(mean(diff ^ 2))
}

计算点到模型间的距离

library(purrr)
hw2_dist<-function(k,b){
measure_distance(k,b,hw2)
}
models<-models%>%
mutate(dist=map2_dbl(k,b,hw2_dist))
models

画图，取距离最小的10条线,颜色最亮的为最优解

ggplot(hw2,aes(x,y))+
geom_point(size=2,na.rm=TRUE)+
geom_abline(aes(intercept=b,slope=k,color = -dist),
data=filter(models,rank(dist)<=10),na.rm=TRUE)

8.请估计模型，marketshare为Y，price和brand 为X.

多变量线性模型lm(y~x1+x2+...,data)

y<-hw2 $marketshare x1<-hw2$ price
x2<-hw2$brand
miki2<-lm(y~x1+x2,hw2)
coef(miki2)
summary(miki2)

#####################统计学基于R 20201017 8:30################

绘制核密度图###

load("D:/mikiStudyFiles/les/example/ch2/example2_2.Rdata")
par(mfrow=c(1,2),cex=0.8,mai=c(0.7,0.7,0.1,0.1))
d<-example2_2$销售额
plot(density(d),col=1,xlab="销售额",ylab="密度",main="")
rug(d,col="pink")
plot(density(d),col="gold",border="black")
polygon(density(d),col="gold",border="black")
rug(d,col="brown")

用lattice包绘制核密度图###

load("D:/mikiStudyFiles/les/example/ch2/example2_3_1.Rdata")
library(lattice)
dp1<-densityplot(~射击环数|运动员,data=example2_3_1,col="blue",cex=0.5,par.strip.text=list(cex=0.6),sub="(a)栅格图")
dp2<-densityplot(~射击环数,group=运动员,data=example2_3_1,auto.key=list(columns=1,x=0.01,y=0.95,cex=0.6),cex=0.5,sub="(b)比较图")
plot(dp1,split=c(1,1,2,1))
plot(dp2,split=c(2,1,2,1),newpage=F)

lattice包绘制的6名运动员射击成绩的和密度图###

不同收入等级的脸谱图###

注意安装包括号内需要加双引号install.packages(""),加载包不需要library()

install.packages("aplpack")
library(aplpack)
load("D:/mikiStudyFiles/les/example/ch2/matrix2_5.Rdata")
faces(matrix2_5,nrow.plot=4,ncol.plot=2,face.type=0)

计算绘制洛伦茨曲线所需的各百分比数值###

install.packages("DescTools")
library(DescTools)
Lc(example2_10 $组中值,example2_10$ 人数)
Lc(example2_10 $组中值,example2_10$ 人数)
Lc(example2_10 $组中值,example2_10$ 人数)

绘制洛伦茨曲线###

plot(Lc(example2_10 $组中值,example2_10$ 人数),xlab="人数比例",ylab="收入比例",col=4,panel.first=grid(10,10,col="gray70"))

ggplot绘图###

三个不同专业的50名学生和统计学的考试分数数据###

set.seed(12345)
n=100
x<-round(rnorm(n,80,6))
y<-round(2+0.4*x+rnorm(n,50,3))
a<-sample(c("管理学","会计学","金融学"),20,replace=T)
b<-sample(c("男","女"),20,replace=T)
d<-data.frame(专业=a,性别=b,数学=x,统计学=y)
head(d)

20201102RHW3第三次作业################

问题 1：

A 和 B 约定在某篮球场见面。他俩都不太守时，出现时间服从均匀分布。他俩也都没有耐心，每个人都会只等对方十分钟就会离开。已知 A 到篮球场的时间为下午 4 点到 5点之间。

(1) 如果 B 到达篮球场的时间也为下午 4 点到 5 点之间，模拟运行 50000 次，看看他们成功相遇的概率。

方法1：陈肯-是否见面建一个表 tibble()表示简单数据框,abs()表示绝对值

n_Sim <- 50000
sim_meet <- tibble(
A = runif(n_Sim, min = 0, max = 60),
B = runif(n_Sim, min = 0, max = 60)
) %>%
mutate(result = ifelse(abs(A - B) <= 10,
"They meet", "They do not"))
p_meet <- sim_meet %>% count(result) %>%
arrange(n) %>%
mutate(percent = n / n_Sim)
p_meet

方法2：miki

library(tidyverse)
a<-runif(50000,0,1)
b<-runif(50000,0,1)
c<-a-b
p<-sum(c>=-1/6&c<=1/6)
p/50000

(2) 对上一问的 50000 次模拟，用不同颜色在一张图中展示成功相遇与否。

d<-c>=-1/6&c<=1/6
hw3_1<-tibble(
A到达时间=a,
B到达时间=b,
AB到达时间差=c,
AB是否相遇=d
)
hw3_point<-ggplot(hw3_1,aes(A到达时间,B到达时间,color=AB是否相遇))+
geom_point()
hw3_point

(3) B 应该如何选择 4 点到 5 点之间的哪个时间段，来提升他们成功相遇的概率？用模拟展示你的理由

分成6段，看哪段的概率高

a1<-runif(50000,0,1)
b1<-runif(50000,0,1/6)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000
a1<-runif(50000,0,1)
b1<-runif(50000,1/6,1/3)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000
a1<-runif(50000,0,1)
b1<-runif(50000,1/3,1/2)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000
a1<-runif(50000,0,1)
b1<-runif(50000,1/2,2/3)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000
a1<-runif(50000,0,1)
b1<-runif(50000,2/3,5/6)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000
a1<-runif(50000,0,1)
b1<-runif(50000,5/6,1)
c1<-a1-b1
p1<-sum(c1>=-1/6&c1<=1/6)
p1/50000

问题 2：请使用 nycflights13 和 pipe 语法

（1）从 flights 数据表中挑选出以下变量：(year, month, day, hour, origin, dep_delay, distance, carrier)，将生产的新表保存为 flight1。

select(data,字段名1,字段名2...)

library(nycflights13)
flight1<-select(flights,year, month, day, hour, origin, dep_delay, distance, carrier)

（2）从 weather 数据表中挑选出以下变量： (year, month, day, hour, origin, humid, wind_speed)，将生产的新表保存为 weather1。

weather1<-select(weather,year, month, day, hour, origin, humid, wind_speed)

（3）将 flight1 表和 weather1 表根据共同变量进行内连接，随机抽取 100000 行数据，将生产的结果保存为 flight_weather。 (提示：sample_n()函数，不用重复抽取)

无放回随机抽样 sample_n(数据表,抽取行数,replace=Flase)

flight_weather<-sample_n(flight1%>%
inner_join(weather1,by=c("year","month","day","hour","origin")),100000, replace = FALSE)
flight_weather

（4）从 flight_weather 表中对三个出发机场按照平均出发延误时间排降序，并将结果保留在 longest_delay 表中。把结果展示出来。

longest_delay<-flight_weather%>%
group_by(origin)%>%
summarise(delay_avg=mean(dep_delay,na.rm=TRUE))%>%
arrange(desc(delay_avg))
longest_delay

（5）根据出发地（origin）在同一个图中画出风速 wind_speed（x 轴）和出发延误时间 dep_delay（y 轴）的平滑曲线图。

geom_smooth()表示平滑曲线

hw3_smooth<-ggplot(flight_weather,aes(wind_speed,dep_delay,color=origin))+
geom_smooth(se = FALSE,na.rm=TRUE)

se=False 去除阴影

hw3_smooth

（6）根据不同出发地（origin）在平行的 3 个图中画出风速 wind_speed（x 轴）和出发延误时间 dep_delay（y 轴）的散点图。

将多个模型的可视化结果放在一张图中 facet_wrap(~变量)

hw3_2_point<-ggplot(flight_weather,aes(wind_speed,dep_delay,color=origin))+
geom_point(na.rm=TRUE)+
facet_wrap(~origin)
hw3_2_point

（7）根据 flight_weather 表，画出每个月航班数的直方分布图，x 轴为月份，y 轴是每个月份航班数所占的比例。

geom_bar()表示条形图,..prop..表示分类变量中每个类别占总量的比，group=1就是将这些类别当作一组的这样一个整体去分别计算各个类别的占比

hw3_bar<-ggplot(flight_weather)+
geom_bar(aes(x=month,y= ..prop.., group = 1))
hw3_bar

（8）根据 flight_weather 表，画出每个月航班距离的 boxplot 图，x轴为月份，y 轴为航行距离, 根据的航行距离的中位数从低到高对 x 轴的月份进行重新排序。

箱线图

hw3_box<-ggplot(flight_weather,aes(x=month,y=distance,group = month))+
geom_boxplot()
hw3_box

根据中位数排序

reorder()重新排序

hw3_box_reorder<-ggplot(flight_weather)+
geom_boxplot(
aes(x=reorder(month,distance,FUN = median),y=distance
))
hw3_box_reorder

问题三：

（1)定义每个区域变种的熵为一定模式，熵是衡量事物分布均匀性的指标，写出一个名为H的函数表示它

先建立矩阵

library(purrr)
a1<-c(0.2,0.2,0.2,0.2,0.2)
a2<-c(0.8,0.1,0.05,0.025,0.025)
a3<-c(0.05,0.15,0.7,0.05,0.05)
hw3_3<-rbind(a1,a2,a3)
hw3_3

给矩阵的字段命名

colnames(hw3_3) <- c("v1","v2","v3","v4","v5")
hw3_3

建立函数function(x,y,..,data)

H<-function(p,hw3_3){
-sum(p*log(p))
}
H

（2）定义区域与区间的距离为一个函数，写一个名为KLD的函数表示它

KLD<-function(p,q,hw3_3){
sum(p*(log(p)-log(q)))
}
KLD
KLD(a1,a2)
KLD(a1,a3)
KLD(a2,a3)

（3）应用purr：：map计算每个区域内变种分布的熵，观察并解释各区该熵的大小顺序

map(x,函数), list()函数合并多个列表成一个列表

hw3_map<-map(list(a1,a2,a3),H)

（4）计算每个区域与其他区域的KL距离，将结果保存为矩阵形式。该矩阵的可行可以代表每个区域到其他区域的变种分布的距离。

方法1：陈肯先建立一个空矩阵，循环计算，将计算得出的KLD值分别放入每个矩阵对应位置

Dm <- matrix( NA , nrow=3 , ncol=3 )
for ( i in 1:3 ) {
for ( j in 1:3 ) {
Dm[i,j] <- KLD( a[[j]] , a[[i]] )
}
}
Dm

方法2：miki-rbind(a,b)按行合并矩阵，cbind(a,b)按列合并矩阵

k1<-KLD(a1,a2)
k2<-KLD(a1,a3)
k3<-KLD(a2,a3)
hw3_KL<-rbind(k1,k2,k3)
hw3_KL

20201213 15:00统计学小组作业-盲盒######

词云

library(wordcloud)
wordcloud<-read.csv("D:/mikiStudyFiles/les_business analyst/wordcloud.csv")
colors=c('red','blue','green','orange','MediumPurple','pink','DarkGray','blue','Gold')
wordcloud(wordcloud $word,wordcloud$ amount,scale=c(5,0.3),min.freq=-Inf,max.words=100,colors=colors,random.order=F,random.color=F,ordered.colors=F)

是否购买男女占比

library(tidyverse)
manghe<-read.csv("D:/mikiStudyFiles/les_business analyst/manghe.csv")
manghe<- manghe %>%
mutate(
是否购买 = ifelse(
是否购买过 == 0,'购买','未购买'
)
)

manghe[complete.cases(manghe[,25]),]
ggplot(manghe[complete.cases(manghe[,25]),]) +
geom_bar(aes(x = 性别, fill = 是否购买),na.rm=TRUE,position="dodge")

mh_new %>%
group_by(年龄,是否购买) %>% # calculate the counts
summarize(counts = n()) %>%
arrange(-counts) # sort by counts

是否购买年龄占比

ggplot(manghe[complete.cases(manghe[,25]),]) +
geom_bar(aes(x = 年龄, fill = 是否购买),na.rm=TRUE,position="dodge")

是否购买职业占比

ggplot(manghe[complete.cases(manghe[,25]),]) +
geom_bar(aes(x = 职业, fill = 是否购买),na.rm=TRUE,position="dodge")

是否购买职业地域占比

ggplot(manghe[complete.cases(manghe[,25]),]) +
geom_bar(aes(x = 省, fill = 是否购买),na.rm=TRUE,position="dodge")

购买过盲用户评价月消费

ggplot(manghe%>%filter(manghe $平均月消费!=""&manghe$ 是否购买过==1)) +
geom_bar(aes(x = 平均月消费),na.rm=TRUE)

抽到隐藏款概率

manghe<- manghe %>%
mutate(
是否抽到隐藏款 = ifelse(
是否抽到过隐藏款== 0,'未抽到过隐藏款','抽到过隐藏款'
)
)

ggplot(manghe[complete.cases(manghe[,26]),]) +
geom_bar(aes(x = 是否抽到隐藏款,y=..prop..,group=1),na.rm=TRUE)

上届大神R：https://www.jianshu.com/p/468862281e19

最后编辑于：2020.12.18 23:33:03

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 194,390评论 5赞 459
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 81,821评论 2赞 371
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 141,632评论 0赞 319
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,170评论 1赞 263
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,033评论 4赞 355
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,098评论 1赞 272
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,511评论 3赞 381
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,204评论 0赞 253
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,479评论 1赞 290
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,572评论 2赞 309
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,341评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,213评论 3赞 312
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,576评论 3赞 298
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,893评论 0赞 17
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,171评论 1赞 250
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,486评论 2赞 341
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,676评论 2赞 335

2020-12-15miki笔记

R数据科学###

20200919 8:30 raining####

shif+enter 换行

后面的年月日等条件捆绑计数

left_join左连接

right_join右连接

inner_join内连接 取交集

full_join全链接 取并集

semi_join保留x表中与Y表中的观测相匹配的所有x

anti_join去除x表中与Y表中的观测相匹配的所有x

sort=TRUE降序

head()前几位

优化的方法optim()

lm()线性模型

coef()调出

线性回归

普通最小二乘

20200920 8:30 经管院b224 raining####

y是均值为0，标准差为2计算10000次的正态分布，x是男女 随机取10000次，y=b+ax

mm<-lm(y~x-1) 减掉1表示去除截距

导入文本文件，注意文本路径""改成"/"

what:声明读入为字符类型数据，可能指定读入的精度/类型

sep：指定各个读入的数据之间的分隔符；默认情况下分隔符：空格、tab；如果不是其它分隔符，例如“：/”通过SEP来指定

可以通过list指定读入变量的变量名，同时生成的对象为列表，则可以同时读入字符与数字

Skip 从第几行开始读入数据

Nlines 指定最大读入行数

如果通过键盘输入的时候，不typeof(hw1_a$ID)希望出现下标提示，则可以使用：quiet=TRUE

encoding =””指定的编码格式，有时候读入的中文可能会出现乱码的时候，可能通过这个参数来指定：Latin-1 或者 UTF-8

统计文本A中的.的个数

统计文本A文字的个数

统计文本A中“will”这个词语的个数

20200927

箱线图适用于连续变量异常值检测

箱线图

ggplot(hw1_fj_dna,aes(,Income))+

geom_boxplot()

去除极端值的箱线图

boxplot(hw1_fj_dna$Income, outline = F)

boxplot(hw1_fj_dna$Income)

极端值

去除极端值

不使用科学计数法

20200925 14:01 RHW1###

1. 请将数据hw1_a和hw1_b分别读入R，查看数据并指出各个变量的形式，最小值，最大值，中值，均值，标准差。

读取hw1_a.csv文件

方法1：陈肯-查看hw1_a各变量类型，最小值，最大值，中值，均值，标准差

直接查看hw1_a中所有字段的类型

直接汇总求出hw1_a中所有字段的max min median mean 1st Qu. 3rd Qu.

方法2：miki

读取hw1_b.csv文件

查看hw1_b各变量类型，最小值，最大值，中值，均值，标准差

2. 结合上课我们所学的几种数据join 的形式，将两个数据集进行合并。对于每种数据合并的方式，请说明key, 并且报告合并后的数据样本总行数。

left_join(data,by="连接字段")，nrow() 表示行数

left_join

right_join

inner_join

full_join

semi_join

anti_join

3. 请筛选出hw1_a 中收入大于4000的样本，并将此样本和hw1_b 中Is_Default=1的样本合并，你可以使用inner join的方式。这一问中你可以用pipe的书写形式。

比较的时候一定要用==，=表示赋值

4. 在第2问的基础上, 请给出Income对Years_at_Employer的散点图，你发现了哪些趋势和现象?

ctr+shif+c多行注释

# #画图 ggplot(data,aes(x,y))+

# # geom_point()

na.rm=TRUE 排除空值

5. 在第4问的基础上 按照Is_Default 增加一个维度，请展示两变量在不同违约状态的散点图。请使用明暗程度作为区分方式

alpha表示明暗程度 as.factor()表示转换成向量

6. 对于第5问，请使用形状作为另外一种区分方式。

shape表示形状

7. 请找出各个列的缺失值，并删除相应的行。请报告每一变量的缺失值个数，以及所有缺失值总数。

sum()表示计数，is.na()判断是否为空值

方法1：陈肯-筛选出hw1_fj第2列不为空的所有记录

方法2：miki

计算缺失值总数

计算每个变量的缺失值数

查看缺失值的行

删除缺失值的行

删除第25列有空值的行

inner_join内连接取交集

full_join全链接取并集

y是均值为0，标准差为2计算10000次的正态分布，x是男女随机取10000次，y=b+ax

5. 在第4问的基础上按照Is_Default 增加一个维度，请展示两变量在不同违约状态的散点图。请使用明暗程度作为区分方式

boxplot.stats(hw1_fj_dna $Income)$ out

hw1_fj_dna_IncomeFilter<-hw1_fj_dna[!hw1_fj_dna $Income%in% boxplot.stats(hw1_fj_dna$ Income)$out]

创建函数function(),sum()计数注意：N(3,10) 表示期望是3，方差是10，标准差是根号10，所以表示为rnorm(n,3,sqrt(10))

4.请按照价格的均值，产生新的变量price_new, 低于均值为“低价格”，高于均值为“高价格”。同样对市场份额也是，产生变量marketshare_new, 数值为“低市场份额”和“高市场份额”