# 3.2.1項
set.seed(60)
(test <- floor(rnorm(40, 50, 10)))
##  [1] 57 54 44 48 32 32 60 53 40 54 47 49 63 44 46 65 48 55 46 54 47 45 56 56 53
## [26] 54 47 32 71 66 47 66 48 44 53 35 42 40 53 48
max(test); min(test)
## [1] 71
## [1] 32
range(test)
## [1] 32 71
summary(test)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   32.00   44.75   48.00   49.85   54.25   71.00
hist(test)

stem(test)
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   3 | 222
##   3 | 5
##   4 | 002444
##   4 | 566777788889
##   5 | 33334444
##   5 | 5667
##   6 | 03
##   6 | 566
##   7 | 1
qqnorm(test)
qqline(test)

sum(test)
## [1] 1994
length(test)
## [1] 40
sum(test)/length(test)
## [1] 49.85
mean(test)
## [1] 49.85
median(test)
## [1] 48
var(test)
## [1] 86.43846
sum((test-mean(test))^2)/length(test); sum((test-mean(test))^2)/(length(test)-1)
## [1] 84.2775
## [1] 86.43846
# 3.2.2項
# ?cars
# cars
plot(cars) # データフレームから散布図を描く

plot(cars[,1],cars[,2]) # 同じ長さの2つのベクトルでも散布図が描ける
cov(cars[,1],cars[,2]); cor(cars[,1],cars[,2]) # 共分散と相関係数
## [1] 109.9469
## [1] 0.8068949
kaiki <- lm(cars[,2]~cars[,1]) # 単回帰直線を求める
summary(kaiki) # 回帰式の情報を得る
## 
## Call:
## lm(formula = cars[, 2] ~ cars[, 1])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## cars[, 1]     3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
abline(kaiki) # 散布図中に回帰直線を引く

# 3.2.3項
# ?iris
# iris
boxplot(iris[1:50,1], iris[51:200,1]) # 箱ひげ図を描く

t.test(iris[1:50,1], iris[51:200,1])
## 
##  Welch Two Sample t-test
## 
## data:  iris[1:50, 1] and iris[51:200, 1]
## t = -15.144, df = 147.39, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.419898 -1.092102
## sample estimates:
## mean of x mean of y 
##     5.006     6.262
# 3.3節
p <- 1/6; q <- 1-p
choose(10,3)*p^3*q^7; dbinom(3,10,p) # 10回中3回1が出る確率
## [1] 0.1550454
## [1] 0.1550454
sum(choose(10,0:3)*p^(0:3)*q^(10-0:3)); pbinom(3,10,p) # 10回中1が出るのが3回以下の確率
## [1] 0.9302722
## [1] 0.9302722
qbinom(0.5, 10, p) # F(x)が0.5以上になる最小のx
## [1] 2
rbinom(5,10,p) # 10回振ることを5セット繰り返したときの1の出る回数のセットを表す乱数
## [1] 4 4 3 0 2
dnorm(175, 171, 5) # x=175での密度関数f(x)の値
## [1] 0.05793831
pnorm(175, 171, 5) # x=175での分布関数F(x)の値
## [1] 0.7881446
qnorm(0.7, 171, 5) # 分位数すなわちF(x)=0.7となるxの値
## [1] 173.622
rnorm(6, 171, 5) # 乱数すなわち確率変数の実現値を6個
## [1] 182.8953 174.9794 173.4342 161.9171 172.1278 170.9544
dnorm(mean=171, 175, sd=5) # 引数の名前を指定することで順番を変える
## [1] 0.05793831
dnorm(2) # 母数を指定しないときはデフォルトが適用され、N(0,1)の密度関数f(x)の値
## [1] 0.05399097
# ?distributions
rnorm(5)
## [1] -0.72403577 -0.39867237 -0.51408822  0.03022384  1.00878882
rnorm(5)
## [1]  1.61075571 -0.58752882  0.15324015 -0.02206241 -0.86460052
set.seed(10); rnorm(5)
## [1]  0.01874617 -0.18425254 -1.37133055 -0.59916772  0.29454513
set.seed(10); rnorm(5)
## [1]  0.01874617 -0.18425254 -1.37133055 -0.59916772  0.29454513
# 3.4.1項
curve(sin(x))

curve(sin(x), xlim=c(-3*pi, 3*pi))
curve(cos(x), add=T)

x <- 1:10; y <- sin(x) # x=(1,2,3,...,10), y=(sin(1),sin(2),...,sin(10))を代入
plot(x,y)

barplot(x)

hist(y)

pie(x)

plot(x)

plot(y)

plot(x,y,type="b", xlab="yoko", main="sin x", xlim=c(2,7), log="x", col="blue", pch="猫", lty=3)

plot(x,y)

oldpar <- par() # 既存のグラフィックスパラメータの保存
par(pch=5)      # グラフィックスパラメータを指定
plot(x,y)

par(oldpar)     # 元のグラフィックスパラメータに戻す
plot(x,y)


# 3.4.2項
plot(x,y)
lines(x,y)
points(x, cos(x))
# text(locator(1),"y=sin x")
# identify(x,y)
# identify(x,y,plot="FALSE")
# i <- identify(x,y)
# i


# 3.5.1項
library() # インストールされているパッケージの一覧
search() # ロードされているパッケージの一覧
## [1] ".GlobalEnv"        "package:stats"     "package:graphics" 
## [4] "package:grDevices" "package:utils"     "package:datasets" 
## [7] "package:methods"   "Autoloads"         "package:base"
# install.packages("evd") # evdパッケージのインストール
library(evd) # evdパッケージのロード
# install.packages("ismev") # evdパッケージのインストール
library(ismev)
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

# install.packages("extRemes") # extRemesパッケージのインストール
library(extRemes)
## Loading required package: Lmoments
## Loading required package: distillery
## 
## Attaching package: 'extRemes'
## The following objects are masked from 'package:evd':
## 
##     fbvpot, mrlplot
## The following objects are masked from 'package:stats':
## 
##     qqnorm, qqplot
library()
search()
##  [1] ".GlobalEnv"         "package:extRemes"   "package:distillery"
##  [4] "package:Lmoments"   "package:ismev"      "package:mgcv"      
##  [7] "package:nlme"       "package:evd"        "package:stats"     
## [10] "package:graphics"   "package:grDevices"  "package:utils"     
## [13] "package:datasets"   "package:methods"    "Autoloads"         
## [16] "package:base"
# 3.6.1項
test <- floor(rnorm(40,50,10))
s <- sqrt(sum((test-mean(test))^2)/length(test))
(47-mean(test))/s*10+50
## [1] 51.39866
z <- (test-mean(test))/s*10+50; z; mean(z)
##  [1] 57.98057 40.42880 50.30167 36.04085 51.39866 66.75646 62.36852 51.39866
##  [9] 64.56249 62.36852 54.68962 43.71975 52.49564 64.56249 59.07756 48.10770
## [17] 30.55592 47.01071 30.55592 40.42880 50.30167 47.01071 44.81674 52.49564
## [25] 51.39866 33.84688 53.59263 64.56249 55.78660 39.33181 38.23482 57.98057
## [33] 34.94387 50.30167 47.01071 65.65948 45.91373 44.81674 63.46550 43.71975
## [1] 50
m <- mean(test)
hensachi <- function(x){(x-m)/s*10+50}
hensachi(47)
## [1] 51.39866
hensachi
## function(x){(x-m)/s*10+50}
hensachi2 <- function(x, group){
 # xは素点 groupは集団全員の得点ベクトル
 m1 <- mean(group) # 集団の平均
 s1 <- sqrt(sum((group-m1)^2)/length(group)) # 集団の標準偏差
 (x-m1)/s1*10+50}
hensachi2(47,test)
## [1] 51.39866
# m1; s1


# 3.6.2項
hensachi <- function(x) {(x-m)/s*10+50}
hensachi
## function(x) {(x-m)/s*10+50}
var
## function (x, y = NULL, na.rm = FALSE, use) 
## {
##     if (missing(use)) 
##         use <- if (na.rm) 
##             "na.or.complete"
##         else "everything"
##     na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs", 
##         "everything", "na.or.complete"))
##     if (is.na(na.method)) 
##         stop("invalid 'use' argument")
##     if (is.data.frame(x)) 
##         x <- as.matrix(x)
##     else stopifnot(is.atomic(x))
##     if (is.data.frame(y)) 
##         y <- as.matrix(y)
##     else stopifnot(is.atomic(y))
##     .Call(C_cov, x, y, na.method, FALSE)
## }
## <bytecode: 0x00000000256bddf8>
## <environment: namespace:stats>
# ?mean
methods(mean)
##  [1] mean,ANY-method          mean,Matrix-method       mean,sparseMatrix-method
##  [4] mean,sparseVector-method mean.Date                mean.default            
##  [7] mean.difftime            mean.POSIXct             mean.POSIXlt            
## [10] mean.quosure*           
## see '?methods' for accessing help and source code
mean.default
## function (x, trim = 0, na.rm = FALSE, ...) 
## {
##     if (!is.numeric(x) && !is.complex(x) && !is.logical(x)) {
##         warning("argument is not numeric or logical: returning NA")
##         return(NA_real_)
##     }
##     if (na.rm) 
##         x <- x[!is.na(x)]
##     if (!is.numeric(trim) || length(trim) != 1L) 
##         stop("'trim' must be numeric of length one")
##     n <- length(x)
##     if (trim > 0 && n) {
##         if (is.complex(x)) 
##             stop("trimmed means are not defined for complex data")
##         if (anyNA(x)) 
##             return(NA_real_)
##         if (trim >= 0.5) 
##             return(stats::median(x, na.rm = FALSE))
##         lo <- floor(n * trim) + 1
##         hi <- n + 1 - lo
##         x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
##     }
##     .Internal(mean(x))
## }
## <bytecode: 0x000000002494b320>
## <environment: namespace:base>
sin
## function (x)  .Primitive("sin")
# 3.6.3項
factsum <- function(n){
  p <- 1; s <- 0 
  for (i in 1:n){ p <- p*i; s <- s+i}
  c(p,s)}
factsum(5)
## [1] 120  15
# 3.6.4項
(d <- as.data.frame(matrix(1:9,3,3)))
##   V1 V2 V3
## 1  1  4  7
## 2  2  5  8
## 3  3  6  9
a <- NULL # aを空の変数とする
for (j in 1:3){a <- c(a, max(d[,j]))} # 既存のaにj列目の最大値を付け加える
a
## [1] 3 6 9
apply(d,2,max) # 列ごとの最大値、2は列を表す
## V1 V2 V3 
##  3  6  9
apply(d,1,max) # 行ごとの最大値、1は行を表す
## [1] 7 8 9
# 3.7.1項
# getwd() # 現在の作業ディレクトリを知る
# setwd("C:/Users/saigo/Documents/R") # 作業ディレクトリを変更する
# getwd()
# setwd("C:/Users/saigo/Documents/R2") # 存在しないフォルダに変更する


# 3.7.2項
rm(list=ls()) # 変数がすべて消去されるので注意
(x <- 3:10)
## [1]  3  4  5  6  7  8  9 10
# q()
x
## [1]  3  4  5  6  7  8  9 10
objects()
## [1] "x"
tenki <- read.csv("tenki.csv")
tenki
##      X 最高気温 最低気温 降水確率
## 1 東京       35       27       40
## 2 札幌       26       19       10
read.csv("tenki2.csv") # 表頭を列名とみなしてしまう
##   東京 X35 X27 X40
## 1 札幌  26  19  10
read.csv("tenki2.csv", header=F) # 表頭を列名とせず、新たな列名をつける
##     V1 V2 V3 V4
## 1 東京 35 27 40
## 2 札幌 26 19 10
read.table("tenki.csv") # スペースやタブで区切られていないので3行1列とみなす
##                            V1
## 1 ,最高気温,最低気温,降水確率
## 2               東京,35,27,40
## 3               札幌,26,19,10
read.table("tenki.csv", sep=",") # 表頭を列名とみなさずデータとみる
##     V1       V2       V3       V4
## 1      最高気温 最低気温 降水確率
## 2 東京       35       27       40
## 3 札幌       26       19       10
read.table("tenki.csv", sep=",", header=T) # 正しい読み込みとなった
##      X 最高気温 最低気温 降水確率
## 1 東京       35       27       40
## 2 札幌       26       19       10
# ?BOD
BOD
##   Time demand
## 1    1    8.3
## 2    2   10.3
## 3    3   19.0
## 4    4   16.0
## 5    5   15.6
## 6    7   19.8
write.table(BOD, "BOD.txt") # txtファイルだが、内容としてはssvとなっている。
write.table(BOD, "BOD.csv") # csvとしたがスペース区切りで正しくないファイルとなる
write.table(BOD, "BOD2.csv", sep=",") # 正しいcsvファイルができる
write.csv(BOD, "BOD3.csv") # 正しいcsvファイルができる
# ?datasets
# saikoro
source("sai.R") # プログラムの読み込み
saikoro
## function (n) 
## {
##     floor(runif(n) * 6 + 1)
## }
saikoro(4)
## [1] 3 5 4 4
x <- 1:10; y <- 1:10
save(x,y, file="suuretsu.Rdata") # xをファイルsuuretsu.Rdataに保存
# rm(x,y)
# x
# y
load("suuretsu.Rdata") # ファイルsuuretsu.Rdataを読み込み
x; y
##  [1]  1  2  3  4  5  6  7  8  9 10
##  [1]  1  2  3  4  5  6  7  8  9 10
x <- 1:10
dput(x, "suuretsu.txt")
dget("suuretsu.txt")
##  [1]  1  2  3  4  5  6  7  8  9 10
# 4章



# 4.2節
library(evd)
curve(drweibull(x, loc=3, scale=2, shape=1.5), xlim=c(-3,3)) # ワイブルの密度関数のグラフ

pfrechet(5, loc=3, scale=2, shape=1.5) # フレシェの分布関数のx=5における値
## [1] 0.3678794
qgumbel(0.2,loc=3,scale=2)
## [1] 2.04823
rgev(5, shape=-1) # 一般極値分布に従う乱数を5個
## [1] -0.2196579  0.3304433 -0.4070970 -2.0283460  0.6692334