# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

# ## 习题 1.3
#
# 数据：

data <- read.table("./ex_1_3.txt", header=TRUE)
data

summary(data)

attach(data)

# 方便直接使用 Year, Nationwide, Rural, Urban 这几个变量。

# ### (1) 
# 求均值、方差、标准差、变异系数、偏度、峰度
#
# **均值**：

c(mean(Nationwide), mean(Rural), mean(Urban))

# 或者，R 提供有一个更方便的函数 `colMeans`，可以直接计算 data frame 各列均值 (`data[-1]` 排除了第一列 Year)：

colMeans(data[-1])

# colMeans 函数等价于如下 apply:

apply(data[-1], MARGIN=2, FUN=mean)

# 注：`MARGIN = 2` 就是取列。

# **方差**：

apply(data[-1], 2, var)

# **标准差**：

apply(data[-1], 2, sd)

# **变异系数**：

# +
cv <- function(x) sd(x)/mean(x)

apply(data[-1], 2, cv)
# -

# **偏度**（skewness）：

# 可以调用一个库完成计算：

# install.packages("psych")
library(psych)

# +
# g1 <- function(x) skew(x, type=2)  # g1、g2 是用 type=2: seehelp(skew)

# or 按照书上（P6）的手写这个函数
g1 <- function(x) {
    n <- length(x)
    
    A <- n / ((n-1) * (n-2))
    B <- 1 / sd(x)^3
    S <- sum((x - mean(x))^3)
    
    A * B * S
}

apply(data[-1], 2, g1)
# -

# **峰度**：

# +
# g2 <- function(x) kurtosi(x, type=2)  # help(kurtosi)

g2 <- function(x) {
    n <- length(x)
    
    A <- (n * (n+1)) / ((n-1) * (n-2) * (n-3))
    B <- 1 / sd(x)^4
    S <- sum((x - mean(x))^4)
    C <- (3 * (n-1)^2) / ((n-2) * (n-3))
    
    A * B * S - C
}

apply(data[-1], 2, g2)
# -

# 注：实际上，psych 包提供了一个 `describe` 函数，可以一次性得到各种常用值：

describe(data[-1], type=2)

# ### (2)
#
# **中位数**：

apply(data[-1], 2, median)

# **四分位距**：

apply(data[-1], 2, quantile)

# 五数：

fn <- apply(data[-1], 2, fivenum)
fn

# 四分位极差：

# +
R1 <- function(Q3, Q1) Q3 - Q1

R1(Q3=fn[4,], Q1=fn[2,])
# -

# 三均值：

# +
M3 <- function(Q1, M, Q3) Q1/4 + M/2 + Q3/4

M3(Q1=fn[2,], M=fn[3,], Q3=fn[4,])
# -

# ### (3) 直方图

histogram <- function(x, xname="x") {
    hist(x, prob=TRUE, main=paste("Histogram of" , xname))
    lines(density(x))
    rug(x) # show the actual data points
}

# layout(matrix(c(1,2,3), nr=1, byrow=T))
histogram(Nationwide, xname="Nationwide")
histogram(Rural, xname="Rural")
histogram(Urban, xname="Urban")

# ### (4) 茎叶图

stem(Nationwide)

stem(Rural)

stem(Urban)

# ### (5) 异常值

abnormal <- function(x) {
    fn <- fivenum(x);
    Q1 <- fn[2];  Q3 <- fn[4];
    
    R1 <- Q3 - Q1
    
    QD <- Q1 - 1.5 * R1
    QU <- Q3 + 1.5 * R1
    
    x[(x < QD) | (x > QU)]
}

apply(data[-1], 2, abnormal)

# 结果为空：没有异常值。

detach(data)

StarcoderdataJupyterScriptsDedupFiltered