# 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