# 必要なライブラリをロード
library(car)

##  要求されたパッケージ carData をロード中です

library(DT)

# カラーパレット（透過色）
COLS <- c(rgb(1, 0, 1, .25), # ピンク
          rgb(0, 0, 1, .25), # ラベンダー
          rgb(0, 1, 0, .25)) # グリーン
GRAY <- rgb(.5, .5, .5, .25) # グレー

# 乱数のシードを設定
set.seed(10)

# 標本サイズ
n <- 100

# 正規乱数とt乱数を生成
xn <- rnorm(n)       # 正規乱数
xt <- rt(n, df = 10) # 自由度10のt乱数

# データの概要を表示
summary(data.frame(xn, xt))

##        xn                xt         
##  Min.   :-2.1853   Min.   :-2.0669  
##  1st Qu.:-0.8291   1st Qu.:-0.5995  
##  Median :-0.1933   Median : 0.0516  
##  Mean   :-0.1365   Mean   : 0.1114  
##  3rd Qu.: 0.5934   3rd Qu.: 0.8164  
##  Max.   : 2.2205   Max.   : 3.7912

# ヒストグラムの階級区切りを定義
bins <- seq(-3, 4, 1) # 階級区切り

# ヒストグラムをプロット
hist(xn, breaks = bins, col = COLS[1], 
     main = 'ヒストグラム 1', 
     xlab = '観測値', 
     ylab = '度数')

hist(xt, breaks = bins, col = COLS[2],
     main = 'ヒストグラム 2', 
     xlab = '観測値', 
     ylab = '度数')

# 四分位数を計算
q1 <- quantile(xn, probs = c(1/4, 2/4, 3/4))
q2 <- quantile(xt, probs = c(1/4, 2/4, 3/4))

# 四分位数の散布図をプロット
matplot(x = q1, y = q2, pch = 16, ylim = c(-10,10),
        main = '散布図（四分位数同士の比較）',
        xlab = '標本1分位数',
        ylab = '標本2分位数')
grid()
abline(a = 0, b = 10, col = 'red')

# Q-Qプロットの準備
p <- function(k) (k - 0.5) / n
k <- 1:n
q <- qnorm(p = p(k)) # 正規分位数

# 標本を順序付け
xn.o <- xn[order(xn)]
xt.o <- xt[order(xt)]

# Q-Qプロットをプロット
matplot(x = q, y = xn.o, type = 'n',
        xlim = c(-4, 4),
        ylim = c(-4, 4),
        main = '正規Q-Qプロット',
        xlab = '正規分位数',
        ylab = '標本分位数')
matpoints(x = q, y = xn.o, pch = 16, col = COLS[1])
matpoints(x = q, y = xt.o, pch = 16, col = COLS[2])
abline(a = 0, b = 1, col = 'red')
legend('topleft', col = COLS[1:2], pch = 16,
       legend = c('標本1（正規乱数)', '標本2（t乱数）'))

# carパッケージを用いたQ-Qプロット
qqPlot(xn,
       main = 'N(0, 1)',
       xlab = '正規分位数',
       ylab = '標本分位数')

## [1] 88 54

qqPlot(xt, 
       main = 't(d.f.=10)',
       xlab = '正規分位数',
       ylab = '標本分位数')

## [1] 52 49

# Shapiro-Wilk正規性検定を実施
shapiro.test(xn)

## 
##  Shapiro-Wilk normality test
## 
## data:  xn
## W = 0.98908, p-value = 0.5911

shapiro.test(xt)

## 
##  Shapiro-Wilk normality test
## 
## data:  xt
## W = 0.96269, p-value = 0.006286