# encording with UTF-8
#
# Chapter 5, 推定の考え方
#

# C:\EconDat\5_est
setwd("C:/EconDat/5_est")

getwd()


#
# 例5.8　医薬品製造企業の売上高の標本統計量
#

## 医薬品製造企業売上高の母集団
df_medsale <- read.csv("med_sales.csv", header=TRUE)
## 単位: 億円, 2016
## df_medsaleを母集団とみなす

summary(df_medsale$sales)

## 観測値の数
n_pop <- length(df_medsale$sales)
n_pop

dev.off() # 今あるグラフを消す
par(ps = 20)          #  文字を大きく
par(mex=1.5)   # 1行の高さを1.5倍とする

## 図5.6

hist(df_medsale[,2], breaks=seq(0,5000,200),freq=F, main="",xlab="売上高(億円)", ylab="密度", col="red")

## 母平均
mean_medsale <- mean(df_medsale$sales)
mean_medsale

## 母分散
var_medsale <- sum((df_medsale$sales-mean_medsale)^2)/n_pop
var_medsale


## 標本

## 表5.2 医薬品製造企業の2016年売上高から抽出した10社の標本
df_med_smpl <- read.csv("med_smpl.csv", header=TRUE)

## 標本の大きさ
n_smpl <- length(df_med_smpl$med_sales)
n_smpl

## 平均
mean_med_smpl <- sum(df_med_smpl$med_sales)/n_smpl
mean_med_smpl

## 標本分散
var_med_smpl <- sum((df_med_smpl$med_sales-mean_med_smpl)^2)/(n_smpl-1)
var_med_smpl 

## 標本標準偏差
se_med_smpl <-sqrt(var_med_smpl)
se_med_smpl

# 図5.8 平均値の標本分布

## 標本抽出実験
s<- 400 # 実験回数
rec<-numeric(s)
set.seed(1234) # 乱数のseedの設定

for(i in 1:s){
  df_medsale_sample <- sample(df_medsale[,2], 10, replace= FALSE) # step1
  # replace = FALSE  sampling without replacement 非復元抽出
  # replace= TRUE  sampling with replacement 復元抽出

  rec[i] <- mean(df_medsale_sample) #step2
}
summary(rec)
mean_mean <- mean(rec)
mean_mean
sd_mean <- sd(rec)
sd_mean

hist(rec, breaks=seq(0,2000,200), freq=F, main="",xlab="平均売上高(億円)", ylab="密度", col="blue")


#
# 例5.9: 都道府県データによる標本平均の分布
#

## 都道府県データの読み込み
df_pref <- read.csv("../data/prefdat/pref_dat_2019.csv", header=TRUE)

## 図5.9 
## 可処分所得（二人以上の世帯のうち勤労者世帯）（家計調査結果）（円）(2017)

summary(df_pref$income)
## 2×1のグラフに表示するためのコマンド
## par(mfrow=c(2,1))
hist(df_pref$income, breaks=seq(300000, 600000, 20000), freq=F, main="", xlab="可処分所得（円）", ylab="密度", col="red")

set.seed(12345) #乱数の種の設定
s<- 400 # 実験回数
rec <- NULL # 空のベクトル
for(i in 1:s){
  income_sample <- sample(df_pref$income, 33, replace= FALSE) 
  # replace = FALSE  sampling without replacement 非復元抽出
  # replace= TRUE  sampling with replacement 復元抽出
  
  rec[i] <- mean(income_sample)
}
summary(rec)
sd(rec)

hist(rec, breaks=seq(420000, 460000,2000), xlab="平均可処分所得(円)", ylab="密度", main="",freq=F, col="blue")



## 図5.10
## 日照時間（年間）（時間）(2016)

summary(df_pref$sunlight)
hist(df_pref$sunlight, breaks=seq(1500,2500,50), xlab="日照時間(年間)(時間)", ylab="密度", main="",freq=F, col="red")

set.seed(12345) #乱数の種の設定
s<- 400 # 実験回数
rec <- NULL # 空のベクトル
for(i in 1:s){
  sun_sample <- sample(df_pref$sunlight, 20, replace= FALSE) 
  # replace = FALSE  sampling without replacement 非復元抽出
  # replace= TRUE  sampling with replacement 復元抽出
  
  rec[i] <- mean(sun_sample)
}
summary(rec)
sd(rec)

hist(rec, breaks=seq(1850,2150,20), xlab="平均日照時間(年間)(時間)", ylab="密度", main="",freq=F, col="blue")

## 図5.11
## 降水量（年間）（mm）(2017)

summary(df_pref$rain)
hist(df_pref$rain, breaks=seq(1000,2800,100), xlab="降水量（年間）（mm）(2017)", ylab="密度", main="",freq=F, col="red")

set.seed(12345) #乱数の種の設定
s<- 400 # 実験回数
rec <- NULL # 空のベクトル
for(i in 1:s){
  rain_sample <- sample(df_pref$rain, 20, replace= FALSE) 
  # replace = FALSE  sampling without replacement 非復元抽出
  # replace= TRUE  sampling with replacement 復元抽出
  
  rec[i] <- mean(rain_sample)
}
summary(rec)
sd(rec)

hist(rec, breaks=seq(1400,1900,50), xlab="平均降水量（年間）（mm）", ylab="密度", main="",freq=F, col="blue")

#
# 例5.10: 医薬品製造企業の売上高の標本平均の標準誤差 
#

## 医薬品製造企業の2016年売上高から抽出した10社の標本
df_med_smpl <- read.csv("med_smpl.csv", header=TRUE)

## 標本の大きさ
n_smpl <- length(df_med_smpl[,2])
n_smpl

## 平均
mean_med_smpl <- sum(df_med_smpl[,2])/n_smpl
mean_med_smpl

## 標本分散
var_med_smpl <- sum((df_med_smpl[,2]-mean_med_smpl)^2)/(n_smpl-1)
var_med_smpl 

## 標本標準偏差
sd_med_smpl <- sqrt(var_med_smpl)
sd_med_smpl

## 平均の標準誤差

var_mean <- var_med_smpl/n_smpl
var_mean

## 平均の標準誤差 
se_mean <- sqrt(var_mean)
se_mean


#
# 例5.12: 平均アルバイト時給の標本分布
# 

## 学生生活アンケート2018
setdw("C:/EconDat/5_est")
df_survey <- read.csv("../data/student_survey2018/student_survey2018_data.csv", header=TRUE)
head(df_survey)
tail(df_survey)

## データの構造を知る
str(df_survey)

## Q45	アルバイトの時給はいくらですか？(単位：円)
wage <- df_survey$Q45

## 欠損値の削除
wage <- na.omit(wage)

## 標本サイズ
n<- length(wage)
n

## 平均
wage_mean <-mean(wage)
wage_mean

## 標準偏差
wage_sd <-sd(wage)
wage_sd

## 平均の標準誤差
wage_mean_se <- wage_sd/sqrt(n)
wage_mean_se

#
# 図5.13: 平均アルバイト時給の標本分布
#

## 平均の標本分布

##  dev.off() # 今あるグラフを消す
##  par(ps = 20)          #  文字を大きく
##  par(mex=1.5)   # 1行の高さを1.5倍とする

curve(dnorm(x, mean=wage_mean, sd=wage_mean_se), 1000, 1150, type="l",  col="black",lwd=3, xlab="アルバイト時給（円)", ylab="確率密度")
abline(v=mean(wage_mean), col="red",lty="dashed", lwd=3)
# text(wage_mean, -0.001, labels="1074.2")
mtext("平均 1074.2", side=1, line=1) # side= 1(下), 2(左), 3(上), 4(右)。line= はグラフから離す行数


#
# 図5.14:正規分布のかたち
#

## dev.off() # 今あるグラフを消す
## par(ps = 20)          #  文字を大きく
## par(mex=1.5)   # 1行の高さを1.5倍とする

## 標準正規分布
curve(dnorm(x, mean=0, sd=1), -7, 7, lty="solid", col="black", lwd=2,  xlab="x", ylab="確率密度")
abline(v=0, col="black",lty="dashed", lwd=2)
text(-3, 0.3,labels = "標準正規分布 N(0,1)",ps=40)

## 平均を2としてみる
curve(dnorm(x, mean=2, sd=1), -7, 7, type="o", lty="dashed", add=TRUE, col="blue", lwd=2)
text(4.5, 0.3,labels = "正規分布 N(2,1)")

##　分散を2にしてみる
curve(dnorm(x, mean=0, sd=2), -7, 7, lty="dotted", add=TRUE, col="red", lwd=3)
text(-4, 0.15,labels = "正規分布 N(0,4)")

## 標準正規分布からの確率
pnorm(0, mean=0, sd=1)

pnorm(1, mean=0, sd=1)

pnorm(1.96, mean=0, sd=1)

#
# 例5.14: 平均アルバイト時給の95%信頼区間
#

## 分布左側の臨界値
z_lower <- qnorm(0.025)
z_lower

## 分布右側の臨界値
z_upper <- qnorm(0.975)
z_upper

## 信頼区間の計算
wage_mean_lower <- wage_mean+z_lower*wage_mean_se
wage_mean_upper <- wage_mean+z_upper*wage_mean_se

wage_mean_confidence_95 <- c(wage_mean_lower, wage_mean_upper)
wage_mean_confidence_95


