---
title: "contrw04"
author: "Saburo Higuchi"
date: "2019/6/7"
output: html_document
---

```{r setup, include=FALSE}
Sys.setenv(LANG="C")
Sys.setlocale('LC_ALL','C')
knitr::opts_chunk$set(echo = TRUE)
```

## データの読み込み

ここは, Cのほうで出力の仕様にしたがっておけば変更しなくて済むはず.
```{r}
d<-read.csv("contrw03.csv",comment.char='#') # 先頭の行を ヘッダ(列名)と見なす. 先頭英字でなかったら先頭にXを補う
d<-d[1:(ncol(d)-1)] # 1列目から(列総数-1) までを残す. コンマで区切られた最後の列は空だから
n<-nrow(d)  # 行数=サンプルサイズ
tseq<-seq(0,ncol(d)-1,1) # 時刻のベクトル c(0,1,2,...,20) を作る
#ヘッダのないとき. read.csv で header=FALSEとして
#xlabel<-paste("X",tseq,sep="") # 列名のベクトル X0,X1,X2,...,X20 を作る
#names(d)<-xlabel # 列名を使う. これしないと列名は V1,...,V21
```

## グラフ

時刻t=1,4,9,16におけるヒストグラム
```{r}
hist(d$X1) # hist(d[,2])でも同じ. 2列目.
hist(d$X4)
hist(d$X9)
hist(d$X16)
hist(d$X9,breaks=seq(-30,20,2),ylim=c(0,100))
hist(d$X16,breaks=seq(-30,20,2),ylim=c(0,100))
```
異なる時刻を比較するには, 最後の2行のように, 共通の階級と範囲を指定するとよい.

breaks=seq(-30,20,2)は, breaks=c(-30,-28,-26,...,20) と同じで階級を定めている.

ylim=c(0,100)は, y軸方向の範囲を定めている.


1,2,3回目のウォーク. 横軸t, 縦軸x
```{r}
plot(tseq,d[1,],ylim=c(-10,10),xlab="t",ylab="x")
plot(tseq,d[2,],ylim=c(-10,10),xlab="t",ylab="x")
plot(tseq,d[3,],ylim=c(-10,10),xlab="t",ylab="x")
```

## 統計量
X(1)=R(1)の標本平均値, 不偏標本分散, 不偏標本標準偏差
```{r}
mean(d$X1) # mean(d[,2])
var(d$X1)
sd(d$X1)
```

すべての時刻(すべての列)の mean, var, sd をまとめて計算する(文法的には理解しなくてよい)
```{r}
xav<-lapply(d,mean)
xvar<-lapply(d,var)
xsd<-lapply(d,sd)
```

標本平均値, 不偏標本分散, 不偏標本標準偏差の時間依存性のグラフ.
```{r}
plot(tseq,xav,xlab="t",ylab="mean x")
plot(tseq,xvar,xlab="t",ylab="var x")
plot(tseq,xsd,xlab="t",ylab="sd x")
boxplot(d,ylab="x")
```

傾きを最小二乗法で求める
```{r}
lm(unlist(xav)~0+tseq) # 最小二乗法で比例定数=係数Coefficientを求める 0+ は切片Interceptを0とおけ
lm(unlist(xvar)~0+tseq) # 最小二乗法で比例定数=係数Coefficientを求める 
```

### 問題1
母ナントカの相当する量については, 上で最小二乗法で求めた直線の傾きはいくつになるか.



### 問題2
X(12)>-1 となる標本比率を求めよう. 必要なら近似を使って, 相当する母比率を求めよう.