盆暗の学習記録

データサイエンスを中心として,日々学んだことの備忘録としていく予定です。初心者であり独学なので内容には誤りが含まれる可能性が大いにあります。

[R]ゼロから作る最小二乗法1:単回帰

「自分で数式をコードに落としていって動かす」という作業は非常に勉強になると思ったので,いろんなアルゴリズムをゼロから作っていきたいと思います。第一弾として単回帰からはじめていきます。

理論の要点を整理してから実装する構成で述べていきます。

最小二乗法だけでなく統計的推論についても触れていきます。

モデル

線形回帰(linear regression)は,予測の目的変数\(Y_i\)と説明変数\(X_{i1}, X_{i2}, ..., X_{ik}\)の間に次のような線形関係を仮定したモデルを置いて予測する手法です。

\[ Y_i = \beta_0 + \beta_1X_{i1} + \beta_2 X_{i2} + \cdots + \beta_k X_{ik} + u_i \]

今回は線形回帰のなかでも説明変数が1つの単回帰モデル(simple regression model) \[ Y_i = \beta_0 + \beta_1X_{i} +u_i \] での場合で考えます。

誤差関数

パラメータ\(\beta_0, \beta_1\)を推定するために多く使われる方法は,実測値\(Y_i\)と予測値\(\hat{Y}_i\)誤差2乗和(Sum of Squared Error)

\[ SSE = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 \]

あるいは平均2乗誤差(Mean of Squared Error: MSE)

\[ MSE = \frac{1}{n} \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 \]

を最小にするパラメータを採用するという最小二乗法(least squares method)です。

どちらの誤差関数も同様にU字型(下に凸)の関数になり,関数の接線の傾きが0になる地点(関数を微分して0になる点)で誤差が最小になるパラメータが得られます1

どちらの誤差関数でも,最小二乗法でのパラメータ推定式は同じものになります。

パラメータの推定

回帰係数のパラメータは「誤差関数(例えばSSE)を微分してゼロになる点」という次の条件 \[ \begin{cases}\frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) = 0\\\frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i = 0\end{cases} \] の解となるものであり, \[ \begin{align} \hat{\beta}_0 &= \bar{Y} - \beta_1 \bar{X}\\ \hat{\beta}_1 &= \frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y_i})} {\sum_{i=1}^n (X_i - \bar{X})^2} \end{align} \] になります(\(\bar{Y},\bar{X}\)はそれぞれ\(Y_i,X_i\)の平均値)。

Rで実装

式に入れるだけなのでシンプルですね

# 関数を定義
SimpleRegression_OLS <- function(x, y) {
  # β_1の推定
  beta_1 = sum((x - mean(x))*(y - mean(y))) / sum((x - mean(x))^2)
  # β_0の推定
  beta_0 = mean(y) - beta_1 * mean(x)
  
  # 結果を格納
  result = c("beta_0" = beta_0, "beta_1" = beta_1)
  return(result)
}

# womenデータセットで試す
SimpleRegression_OLS(x = women$height, y = women$weight)
##    beta_0    beta_1 
## -87.51667   3.45000
# lm()と比較
res <- lm(weight ~ height, women)
res
## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Coefficients:
## (Intercept)       height  
##      -87.52         3.45

同じ値が出ました。

統計的検定

:誤差項は相互に独立に平均0,分散正規分布に従う

という仮定を置いた単回帰モデルをsimple normal regression modelと呼び,そこでは次のような方法で統計的推論を行うことができます。

標準誤差

定量\(\hat{\beta}_0, \hat{\beta}_1\)の標準誤差は で,ここでは誤差項の分散の推定量(残差分散)であり,残差平方和を\(n-2\)で割ったもの です。

検定統計量

検定統計量は次のようになります2\[ t = \frac{\hat{\beta}_1 - \beta_1}{SE(\hat{\beta}_1)} \]

これは自由度\(n-2\)\(t\)分布に従うことが知られています。

仮説体系は多くの場合 \[ H_0: \beta_1 = 0, \hspace{2em} H_1: \beta_1 \neq 0 \] なので,実際に使う検定統計量は \[ t = \frac{\hat{\beta}_1}{SE(\hat{\beta}_1)} \] になります。

そして,有意水準\(\alpha\)のもとで次の判定を行います。

あるいは\(p\)値の算出をします。

Rで実装

まず標準誤差を出します

# 標準誤差
StandardError <- function(x, y){
  # 残差の取得
  coef = SimpleRegression_OLS(x,y)
  y_hat = coef[1] + coef[2] * x
  u = y - y_hat
  
  # 誤差項の推定
  n = length(x)
  sigma_hat = sum(u^2)/(n-2)
  
  # 標準誤差(β_0)の推定
  A = (1/n)+(mean(x)^2/sum((x - mean(x))^2))
  SE_0 = sqrt(A * sigma_hat)
  
  # 標準誤差(β_1)の推定
  SE_1 = sqrt(sigma_hat / sum((x - mean(x))^2))
  
  # 結果の格納
  result = c("SE_0" = SE_0, "SE_1" = SE_1)
  return(result)
}

# womenデータセットで試す
StandardError(x = women$height, y = women$weight)

次に検定を処理します

tTest <- function(x, y){
  # 係数を算出
  beta = SimpleRegression_OLS(x, y)
  # 標準誤差を算出
  SE = StandardError(x, y)
  # 検定統計量を算出
  t_0 = beta[1] / SE[1]
  t_1 = beta[2] / SE[2]
  # p値
  n = length(x)
  p_0 = pt(q = -abs(t_0), df = (n-2))*2 # 両側検定なので2倍する
  p_1 = pt(q = -abs(t_1), df = (n-2))*2 # 両側検定なので2倍する
  # 結果
  result = matrix(c(t_0, t_1, p_0, p_1), ncol = 2,
                  dimnames = list(c("beta_0","beta_1"),c("t_value","p_value")))
  return(result)
}

# womenデータセットで試す
tTest(x = women$height, y = women$weight)
##          t_value      p_value
## beta_0 -14.74103 1.711082e-09
## beta_1  37.85531 1.090973e-14

lmのほうと比較してみます

summary(res)
## 
## Call:
## lm(formula = weight ~ height, data = women)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7333 -1.1333 -0.3833  0.7417  3.1167 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
## height        3.45000    0.09114   37.85 1.09e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared:  0.991,  Adjusted R-squared:  0.9903 
## F-statistic:  1433 on 1 and 13 DF,  p-value: 1.091e-14

ちゃんと同じ値になっています。

おまけ

OLS推定量の導出

SSEを誤差関数にした最小二乗推定量の導出についてメモ。

まず, \[ SSE = \sum_{i=1}^n (Y_i - \hat{Y}_i)^2 = \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)^2 \]微分は,合成微分 (f \circ g)' = f'(g(x)) g'(x)から, \[ \frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) \\ \frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i \] となる。

誤差を最小化するパラメータは次の条件の解である \[ \begin{cases} \frac{\partial SSE}{\partial \beta_0} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i) = 0 \\ \frac{\partial SSE}{\partial \beta_1} = -2 \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 X_i)X_i = 0 \end{cases} \] これを整理すると正規方程式(normal equations) \[ \begin{cases} \begin{align} n \beta_0 + \sum_{i=1}^n \beta_1 X_i &= \sum_{i=1}^n Y_i \\ \sum_{i=1}^n \beta_0 X_i + \sum_{i=1}^n \beta_1 X_i^2 &= \sum_{i=1}^n X_iY_i \end{align} \end{cases} \] が得られ,このうち上の式に\((1/n)\sum_{i=1}^n X_i\)を乗じて \[ \begin{cases} \begin{align} \sum_{i=1}^n \beta_0 X_i + \frac{1}{n}\beta_1 (\sum_{i=1}^n X_i)^2 &= \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \\ \sum_{i=1}^n \beta_0 X_i + \sum_{i=1}^n \beta_1 X_i^2 &= \sum_{i=1}^n X_iY_i \end{align} \end{cases} \] 下の式から引けば \[ \begin{align} \sum_{i=1}^n \beta_1 X_i^2 - \frac{1}{n}\beta_1 (\sum_{i=1}^n X_i)^2 &= \sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \\ (\sum_{i=1}^n X_i^2 - \frac{1}{n}(\sum_{i=1}^n X_i)^2)\beta_1 &= \sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i) \end{align} \] となるから,\(\beta_1\)\[ \begin{align} \beta_1 &= \frac{\sum_{i=1}^n X_iY_i - \frac{1}{n}(\sum_{i=1}^n X_i)(\sum_{i=1}^n Y_i)} {\sum_{i=1}^n X_i^2 - \frac{1}{n}(\sum_{i=1}^n X_i)^2} \\ &= \frac{\sum_{i=1}^n(X_i-\frac{1}{n}\sum_{i=1}^n X_i)(Y_i-\frac{1}{n}\sum_{i=1}^n Y_i)} {\sum_{i=1}^n (X_i - \frac{1}{n}\sum_{i=1}^n X_i)^2} \\ &= \frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y_i})} {\sum_{i=1}^n (X_i - \bar{X})^2} \end{align} \] となる。

一方,\(\beta_0\)の値は正規方程式の上の式を\(n\)で割れば \[ \beta_0 + \frac{1}{n}\sum_{i=1}^n \beta_1 X_i = \frac{1}{n}\sum_{i=1}^n Y_i \] から \[ \begin{align} \beta_0 &= \frac{1}{n}\sum_{i=1}^n Y_i - \frac{1}{n}\sum_{i=1}^n \beta_1 X_i \\ &= \bar{Y} - \beta_1 \bar{X} \end{align} \] として得られる。

参考

経済分析のための統計的方法

経済分析のための統計的方法


  1. こういう形の誤差関数は少しずつ誤差を下げていく勾配降下法のような数値解でも,いきなり誤差の最小値を算出して使う最小二乗法のような解析解でも解くことができます。

  2. 切片の係数\(\hat{\beta}_0\)のほうも同じ形の式で検定統計量が算出されます

愛媛県に行ってきた

f:id:nigimitama:20190209235824j:plain

愛媛県に行って松山市をうろうろしたので,せっかくなので印象に残ったこととかをメモしておきます。

(記憶力が極めて貧弱なので,せっかくどこか行って写真撮っても言葉とともに整理して記録しないと忘れちゃうんですよね…)

松山城

県庁のすぐ後ろに松山城があったので登ってみました。

道中の坂や石を積んだ簡素な階段はコンクリートで舗装されて歩きやすいようになっていて,地元の人や中国人観光客がたくさんいました。

f:id:nigimitama:20190209234137j:plain

f:id:nigimitama:20190209234146j:plain

お城の方は外壁の色合いなどがすごく綺麗でした。

f:id:nigimitama:20190209234159j:plain

内装も綺麗に保存されていました。

f:id:nigimitama:20190209234217j:plain

皇居とかお城とかに行くと,少し先に現代的建築物が見えたりして時代が混じった景色が面白いですよね。

f:id:nigimitama:20190209234228j:plain

屋根にはシャチホコが載ってるタイプのお城でした

萬翠荘

萬翠荘(ばんすいそう)は久松定謨(ひさまつさだこと)伯爵が1922年に別邸として建てたフランス風の建物で,外装は残念ながら修繕中だったのですが内装は問題なく見ることができました

f:id:nigimitama:20190209234255j:plain

f:id:nigimitama:20190209234739j:plain

天皇陛下がお泊りになった際にここで朝食を召し上がられた」という貴賓室。この部屋はやけに狭かったですが,いい雰囲気でした。

シャワー付きの客室や水洗式トイレなど,当時の最先端の建築をしていたようです。(一方で工費は30万円=現在の価値で19億円だったとか)

Assassin's Creed Unityというゲームをプレイしたとき,画面上に出てくるフランス式の建物に感動していた私にとってはとてもいい場所でした(下の画像はゲーム画面)。

f:id:nigimitama:20190209235146j:plain

道後温泉

道後温泉本館にも行ってきました。明治時代の建物のようですが3階建ての立派な建物でした。

ただ,今は保存修理工事中とかで1階の一部の湯しか入れませんでしたが…。

f:id:nigimitama:20190209235502j:plain

ここはまぁ,正直なところ外装だけみて楽しむべきでした。内装は微妙です。

地元の人も多く利用されているようで,背中に大きな入れ墨を入れた人や,湯の吹出口に股間をすりつけて気持ちよさそうにしているお爺さんとかが居てカオスでした…

伊予鉄道

町中に路面電車が走っていて,道のほうも線路の面積を考慮しているためにすごく幅が広くて,電車用の信号機があって,電車用の電線がたくさん張ってあって…といった雰囲気がなかなか良かったです

f:id:nigimitama:20190209235516j:plain

一部の車両はレトロな形をしていて,実用性と観光資源とを兼ね備えていそうな感じがすごく良かったです。

f:id:nigimitama:20190209235534j:plain

この車両は夏目漱石坊っちゃん』にも登場したタイプの車両ゆえに坊っちゃん列車と呼ばれているのだとか。いい車両です。

f:id:nigimitama:20190209235548j:plain

街が綺麗

他に印象的だったのは,街が全体的に綺麗だったことです。

f:id:nigimitama:20190209235623j:plain

上の図の建物は愛媛県庁です。歴史のありそうな立派な建物でした。敷地内にはみかんの木もいくつか植えられていて,実もいっぱい生ってました。

f:id:nigimitama:20190209235639j:plain

また県庁周辺は道路が広くて,道路が直線的に格子状になっていて,「計画的に開発された都市」という感じがしました。

普段は東京のゴチャゴチャした道を見ているので,こういう整備された町並みを見るだけでも少しテンションが上がります。

f:id:nigimitama:20190209235702p:plain

私は貧乏学生で旅行はほとんどしてこなかったので,今回はすばらしい卒業旅行を体験できてよかったです。

[R]交差項や2乗項を作る

どうやるのかちょっと悩んだのでメモ。

交差項

交差項等を作りたいときはmodel.matrix()を使うといいようです。

object引数にformulaを指定して,変数名を*でかけ合わせて交差項を作ります。

library(tidyverse)
model.matrix(object = ~ Sepal.Length * Sepal.Width, data = iris) %>% head()
##   (Intercept) Sepal.Length Sepal.Width Sepal.Length:Sepal.Width
## 1           1          5.1         3.5                    17.85
## 2           1          4.9         3.0                    14.70
## 3           1          4.7         3.2                    15.04
## 4           1          4.6         3.1                    14.26
## 5           1          5.0         3.6                    18.00
## 6           1          5.4         3.9                    21.06

切片を除く

-1を入れれば切片は作られません

model.matrix(object = ~  -1 + Sepal.Length * Sepal.Width, data = iris) %>% head()
##   Sepal.Length Sepal.Width Sepal.Length:Sepal.Width
## 1          5.1         3.5                    17.85
## 2          4.9         3.0                    14.70
## 3          4.7         3.2                    15.04
## 4          4.6         3.1                    14.26
## 5          5.0         3.6                    18.00
## 6          5.4         3.9                    21.06

もっと多数の場合

「手持ちのデータセットのうち,2つの変数の交差項だけを,すべての変数について組み合わせる」というようなことをやりたい場合は次のようにやればよさそうです。

# 2つの変数名の全ての組み合わせ
colnames_df <- expand.grid(colnames(iris), colnames(iris))

# formulaにする
fmla <- as.formula(str_c("~ -1 + ", str_c(colnames_df[,1], " * ", colnames_df[,2], collapse = " + ")))
> fmla
## ~-1 + Sepal.Length * Sepal.Length + Sepal.Width * Sepal.Length + 
##     Petal.Length * Sepal.Length + Petal.Width * Sepal.Length + 
##     Species * Sepal.Length + Sepal.Length * Sepal.Width + Sepal.Width * 
##     Sepal.Width + Petal.Length * Sepal.Width + Petal.Width * 
##     Sepal.Width + Species * Sepal.Width + Sepal.Length * Petal.Length + 
##     Sepal.Width * Petal.Length + Petal.Length * Petal.Length + 
##     Petal.Width * Petal.Length + Species * Petal.Length + Sepal.Length * 
##     Petal.Width + Sepal.Width * Petal.Width + Petal.Length * 
##     Petal.Width + Petal.Width * Petal.Width + Species * Petal.Width + 
##     Sepal.Length * Species + Sepal.Width * Species + Petal.Length * 
##     Species + Petal.Width * Species + Species * Species

このformulaをmodel.matrix()に入れます

# これをmodel.matrix()に入れると全組み合わせができる
origin_dumy_interact <- model.matrix(object = fmla, data = iris)

# (1)元の変数,(2)カテゴリカル変数をダミー変数にしたもの,(3)交差項の3つが入っている
origin_dumy_interact %>% as_data_frame() %>% glimpse()
## Observations: 150
## Variables: 21
## $ Sepal.Length                     <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4,...
## $ Sepal.Width                      <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9,...
## $ Petal.Length                     <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7,...
## $ Petal.Width                      <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4,...
## $ Speciessetosa                    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ Speciesversicolor                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ Speciesvirginica                 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Sepal.Length:Sepal.Width`       <dbl> 17.85, 14.70, 15.04, 14.26, 1...
## $ `Sepal.Length:Petal.Length`      <dbl> 7.14, 6.86, 6.11, 6.90, 7.00,...
## $ `Sepal.Length:Petal.Width`       <dbl> 1.02, 0.98, 0.94, 0.92, 1.00,...
## $ `Sepal.Length:Speciesversicolor` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Sepal.Length:Speciesvirginica`  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Sepal.Width:Petal.Length`       <dbl> 4.90, 4.20, 4.16, 4.65, 5.04,...
## $ `Sepal.Width:Petal.Width`        <dbl> 0.70, 0.60, 0.64, 0.62, 0.72,...
## $ `Sepal.Width:Speciesversicolor`  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Sepal.Width:Speciesvirginica`   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Petal.Length:Petal.Width`       <dbl> 0.28, 0.28, 0.26, 0.30, 0.28,...
## $ `Petal.Length:Speciesversicolor` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Petal.Length:Speciesvirginica`  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Petal.Width:Speciesversicolor`  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ `Petal.Width:Speciesvirginica`   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...

できました。交差項だけでなく,交差項に使った元の変数も(factor,characterはダミー変数化されて)入っていますね。

交差項だけ取り出したい場合は":"でselectすればいいし,このまま使っても良い‥という感じですね。

留意点

ただ,留意点はmodel.matrix()がna.omit()のようなNAのリストワイズ除去も行ってしまう点ですね…。

一応,出力されるmatrixのrownamesを使えば何行目が除去されたかがわかるので,この情報を使えば他のデータと行を合わせることは容易にできます。

2乗項

ダミー変数や数値でない変数を避けつつ2乗項をつくります。

model.matrix()のように一発でやってくれる関数を知らないため,自作関数をselect_if()mutate_all()rename_all()に通していくことで処理してみます。

# 関数を定義 ----------------------------------------
# 値が2水準より多いかどうかを判定する関数(select_if()用)
is_more_than_2_level <- function(x){
  is_more_than_2_level = unique(x) %>% na.omit() %>% length() > 2
  return(is_more_than_2_level)
}
# 2乗する関数(mutate_all()用)
as_quadratic_term <- function(x){
  x^2
}
# 2乗項であることを示す名前にする(rename_all()用)
name_add_sqr <- function(name) {
  newname <- str_c("sqr_",name)
  return(newname)
}

# 処理の実行 -----------------------------------------
quadratic_terms <- iris %>% 
  # ダミー変数でないものを取り出し,数値のものを取り出す
  select_if(is_more_than_2_level) %>% select_if(is.numeric) %>% 
  # 2乗する
  mutate_all(as_quadratic_term) %>% 
  # 2乗項であることを示す名前にする
  rename_all(name_add_sqr)
> head(quadratic_terms)
##   sqr_Sepal.Length sqr_Sepal.Width sqr_Petal.Length sqr_Petal.Width
## 1            26.01           12.25             1.96            0.04
## 2            24.01            9.00             1.96            0.04
## 3            22.09           10.24             1.69            0.04
## 4            21.16            9.61             2.25            0.04
## 5            25.00           12.96             1.96            0.04
## 6            29.16           15.21             2.89            0.16

できました。

使用を避けるべきグラフ

最近グラフつくってて思ったことをメモ

円グラフを使わない

「円グラフわかりにくいから,棒グラフを使え」という話です(人間は面積より長さのほうが大小関係を把握できるため)。

例えば以下の円グラフではgear = 3と gear = 4の差がパッと見にはわからないと思いますが,棒グラフではひと目で大小関係がわかります。

f:id:nigimitama:20190208224540p:plain

この主張は最近では色んな所で目にするようになり,データ分析に関わる人の間では大部分の人の合意がとれていそうに思います*1

帯グラフもあまり良くないかも

それに加えて最近私が思うのが,「帯グラフも棒グラフでいいんじゃね?」説です。

帯グラフのカテゴリが2つなら非常に有効だと思うのですが,カテゴリが3以上のものになると始点が異なり始めるので,それだったら棒グラフにして始点を統一したほうがいいのではないかな,と。

f:id:nigimitama:20190208224551p:plain

例えば,上の帯グラフでの「直列エンジンのgear = 4」と「V型エンジンのgear = 3」では,どちらの割合が大きいかがパッと見はわからないと思います。

でも,棒グラフにすれば始点が揃うので差がわかるようになります。

棒グラフのほうがちょっと面積はとってしまいますが,誠実なグラフ(トリックをしにくいグラフ)は棒グラフかなと思います。

参考:グラフづくり全般で参考になる記事

とくに後者の記事の方は本当に細かいポイントまでわかりやすくまとめられていていいです。

(参考)記事中のグラフのコード

せっかくコード書いたので載せておきます。

円グラフと棒グラフ

library(tidyverse)
library(RColorBrewer)
library(gridExtra)

pie <- ggplot(mtcars, aes(x = "", fill = factor(gear))) +
  geom_bar(width = 1)+
  coord_polar(theta = "y", direction = -1)+ # bar to pie
  scale_fill_brewer(palette = "Blues")+ # color
  theme(axis.text.x=element_blank())+ # 文字を消す
  labs(x = "", y = "", fill = "gear",
       title = "ギアの数")

bar <- ggplot(mtcars, aes(x = gear)) +
  geom_bar(fill = brewer.pal(3, "Blues")[3]) + 
  labs(title = "ギアの数")

grid.arrange(pie, bar, ncol = 2)

帯グラフと棒グラフ

# 集計
df <- mtcars %>% mutate(vs = case_when(vs == 0 ~ "V型エンジン", vs == 1 ~ "直列エンジン"),
                        gear = as.factor(gear))
df_m <- df %>% group_by(vs) %>% count() %>% rename(m = n)
table_df <- df %>% group_by(gear, vs) %>% count() %>% full_join(df_m) %>% mutate(p = n / m)

# plot
obi <- ggplot(table_df, aes(x = vs, y = p, fill = gear)) + 
  geom_bar(stat = "identity", position = "fill", color = "White") + 
  scale_fill_brewer(palette = "Blues") + 
  coord_flip() +
  labs(x = "", y = "")

bar <- ggplot(table_df, aes(x = vs, y = p, fill = gear)) + 
  geom_bar(stat = "identity", position = "dodge", color = "White") + 
  scale_fill_brewer(palette = "Blues") + 
  coord_flip() +
  labs(x = "", y = "")

grid.arrange(obi, bar)

[R]BBCニュースで使われるグラフを描く{bbplot}パッケージ

BBCニュースで使われるグラフのスタイルを再現する{bbplot}というパッケージの存在を知りました

github.com

ちょっと使ってみたいと思います。

インストールと実行

ggplotのテーマを一括でセットするbbc_style()という関数がメインの機能のようです。

使い方はggplotの関数に+でつなげるだけ。

# install.packages('devtools')
devtools::install_github('bbc/bbplot')

library(bbplot)
library(tidyverse)

ggplot(longley, aes(x = Year, y = GNP)) +
  geom_line(colour = "#007f7f", size = 1) +
  geom_hline(yintercept = 0, size = 1, colour="#333333") +
  labs(title = "GNP",
       subtitle = "GNP in US, 1947-1962") +
  bbc_style()

f:id:nigimitama:20190208191025p:plain

Windowsだと警告がでる

ただし,Windowsの場合,次のような警告がでます

Warning message:
In grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
   Windows のフォントデータベースにフォントファミリが見付かりません 

bbplot()の中身を見てみると,Macの英字フォントのHelveticaが指定されています。

function () 
{
  font <- "Helvetica"
  ggplot2::theme(plot.title = ggplot2::element_text(family = font, 
(中略)
}

上のグラフでは警告を出しつつもHelveticaの代わりにArialフォントを使って実行していたようです。

fontを游ゴシックに書き換えます

# 1. フォントファミリーを定義
windowsFonts(YuGothic = windowsFont("游ゴシック"))

# 2. bbc_styleの書き換え
bbc_style2 <- function () 
{
  font <- "YuGothic"
  ggplot2::theme(plot.title = ggplot2::element_text(family = font, 
    size = 28, face = "bold", color = "#222222"), plot.subtitle = ggplot2::element_text(family = font, 
    size = 22, margin = ggplot2::margin(9, 0, 9, 0)), plot.caption = ggplot2::element_blank(), 
    legend.position = "top", legend.text.align = 0, legend.background = ggplot2::element_blank(), 
    legend.title = ggplot2::element_blank(), legend.key = ggplot2::element_blank(), 
    legend.text = ggplot2::element_text(family = font, size = 18, 
      color = "#222222"), axis.title = ggplot2::element_blank(), 
    axis.text = ggplot2::element_text(family = font, size = 18, 
      color = "#222222"), axis.text.x = ggplot2::element_text(margin = ggplot2::margin(5, 
      b = 10)), axis.ticks = ggplot2::element_blank(), 
    axis.line = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(), 
    panel.grid.major.y = ggplot2::element_line(color = "#cbcbcb"), 
    panel.grid.major.x = ggplot2::element_blank(), panel.background = ggplot2::element_blank(), 
    strip.background = ggplot2::element_rect(fill = "white"), 
    strip.text = ggplot2::element_text(size = 22, hjust = 0))
}

折れ線グラフ

ggplot(longley, aes(x = Year, y = GNP)) +
  geom_line(colour = "#007f7f", size = 1) +
  geom_hline(yintercept = 0, size = 1, colour="#333333") +
  labs(title = "アメリカのGNPの推移(1947-1962)") +
  bbc_style2()

f:id:nigimitama:20190208191040p:plain

英語でなく日本語になるとやっぱり印象が変わりますが,結構綺麗…なのかな?

軸の線を消して横軸だけつくるというのがBBCスタイルの模様

棒グラフ

# 前処理
temp <- infert %>% mutate(case = ifelse(case == 1, "症例群", "対照群"))
# plot
ggplot(temp, aes(x = education, fill = case)) +
  geom_bar(position = "dodge", size = 1) +
  geom_hline(yintercept = 0, size = 1, colour="#333333") +
  labs(title = "教育年数",
       subtitle = "群別の被験者の教育年数") +
  scale_fill_manual(values = c("#1380A1", "#FAAB18")) +
  bbc_style2()

f:id:nigimitama:20190208191135p:plain

BBCのグラフはsubtitleも使われるので,そこも真似るとより綺麗になるかも。

ただ,このようにBBC的な色合いを持ってくるのはbbc_style()の仕事ではないので,公式のサンプルコードを見ながら色合いを自分で指定するしかないというのがちょっと面倒ですね

まとめ

  • 文字サイズとかのスタイルをBBC風に一括指定できる
    • BBC風の色使いやグラフの装飾は手動
  • でもこれ,Helveticaフォントが格好いいだけでは‥?
    • Windowsだと再現率が低い
    • 日本語だと再現率はさらに下がる
  • 「美しいグラフとは何か」について考える上では大いに参考になる
    • BBCの実際の記事ではplotly的な動的なグラフになっているのでこれが再現できるとさらに綺麗なグラフに近づきそう

参考

[Python]パソコンのバッテリがおかしいのでログを取る

動機

昨年の11月にVAIO S13を買ったのですが,購入当初から「バッテリが少なくなっています」という警告が出ずにいきなりバッテリ切れでスリープに入ることがありました。

指紋リーダーもフリーズすることがよくあったので,「ソフトウェア(Windows?)のインストールに失敗してるのかな」と考えてとりあえず放置していました

Windowsのアップデートをいくつか経て,今はバッテリ残量低下の警告がでるようになってきたのですが,電池残量が30~20%のところから一気に0%に減ってPCがダウンすることがたまにあります。

ただ,毎回そうなるわけではなく,いまいち再現性がない不具合なのでログをとってみたいと思うようになりました。

pythonによるバッテリ残量の取得

バッテリに関するセンサーが取得した情報はpsutilパッケージで得られるようです。

psutil documentation — psutil 5.5.0 documentation

で,日時とバッテリ残量(%),充電器接続の有無といった情報を一定時間ごとに取得してcsvに保存するコードを書いてみました(私はエンジニアリング技術が皆無なのでたぶんスマートなコードじゃないと思いますが…)

import psutil
import datetime
import pandas as pd
import os
import time

i = 0
while True:    
    # バッテリ情報
    btr = psutil.sensors_battery()
    # 現在時刻情報
    dt_now = datetime.datetime.now() # date and time
    # make DF
    log_i = pd.DataFrame(data = {"date": [dt_now.strftime('%Y-%m-%d')], "time": [dt_now.strftime('%H:%M:%S')], 
                                "battery_left_percent": [btr.percent], "battery_power_plugged": [btr.power_plugged]})
    
    # print
    print(log_i)

    # DF binding
    if i == 0:
        log = log_i
    else:
        log = pd.concat([log, log_i])
    
    # save
    file_name = "log_"+dt_now.strftime('%Y-%m-%d')+".csv"
    log.to_csv(file_name, index = None)

    # sleep
    time.sleep(60)
    i = i + 1

データの様子

取得したcsvデータを使って差分をとってバッテリ残量の変化量を算出し,plotするとこんな感じ。

f:id:nigimitama:20190206185445p:plain

-15前後の2つの点は,1分間にバッテリが約15%減少したことを示します。

やっぱおかしいですね。サポートに連絡します…

パソコン蛾物故割れたのでMBを替えた

パソコン(のBIOS)が起動しなくなって,いろいろと試行錯誤した結果,マザーボード(MB)の載せ替えで直った,という経験をしたのでメモしておきます。

問題発生時の状況

グラボの載せ替え

昨年末にGTX1060を買いまして,載せ替えました。

ところが,載せ替えたらパソコンが起動しなくなりました。

このときは1060を外して元の構成にしたりグラボを外して最小構成にしたりすればWindowsが起動できたので,

「1060に初期不良があったのだろう」と考えてサポートに連絡し,新品と交換してもらいました。

Windowsの破損

その後,新品の1060が来るまでの間にGTX760(元のグラボ)に戻して起動しようとしたところ, 「Windowsの起動に使う一部のファイルが読み込めない」という表示がでて Windowsが起動できなくなりました

復元ポイントがあれば復元できたみたいですがそういうのは作ってませんでした…

元々購入から4年経つ寿命間近なPC(パソコンの法定耐用年数は4年)で,いろいろとガタが来ているPCだったので,いくつかパーツが壊れてるのかなと判断し,問題の特定とパーツの交換を行っていこうと考えました。

問題発生時のPCの状態

なお,このときの「壊れている可能性があり,かつ,起動に関係すると思われるパーツ」は以下の通り。ボロボロです。

  • マザーボード(MB)
    • 購入からの経過年数:4年
    • ホコリ掃除のたびにBIOSの時刻がリセットされたり起動が不安定になる状態にあった(電池切れ?)
  • HDD1(OSが入っている方)
    • 購入からの経過年数:4年
    • これまで問題なく動いてきたし,CrystalDiskInfo上のエラーチェックでも「正常」と診断されてきた
  • HDD2
    • 購入からの経過年数:8年(前のPCから引き継いだ)
    • CrystalDiskInfo上で「注意」と表示されてて,そろそろヤバいなと思っていた

試した解決策1:Windowsの修復

Windows10のインストールプログラムWindowsの修復ができる)をUSBメモリに入れてパソコンを起動しようとしたところ,

BIOSすら起動しなくなりました

試した解決策2:HDDの取り替え

HDDのSATAケーブルも全部はずして起動するとBIOSを起動することができたので,「繋いだMBに悪影響を与えるくらいにHDDが壊れている」という仮説を立ててHDDを買い換えることにしました。

今の時期は4TBが一番コスパが良さそうだったのでこれにしました。

結果:HDDを取り替えてもBIOSは起動せず

そして新しいHDDを取り付け,そのHDDのSATAケーブルのみをMBに接続して起動しようとしましたが, 相変わらずSATAケーブルをつないでいるとBIOSも起動しない状態でした

試した解決策3:MBの電池交換

マザーボードにはBIOSの設定などの情報を保持させておくのに必要な電力を供給するボタン電池がついています。

ネットで試した解決策を探していた時,「マザーボードボタン電池が尽きるとパソコンが起動しなくなることがある」という旨のことが書いてあったのでそれを参考に電池交換も行うことに決めました。

結果:電池を替えても状況は改善せず

近所のスーパーでCR2032Pを買って交換しましたが…ダメでした

そもそもBIOSが起動しないときにMBの電池交換は不要だった

後で別のサイトをみてわかったのですが,「マザーボードボタン電池が尽きるとパソコンが起動しなくなることがある」という記述の「パソコンが起動しない」はWindowsが起動しないことを指していて*1BIOSが起動しないことを意味しているわけではありませんでした…早とちりしました。

試した解決策4:CMOSクリア

BIOSが起動しない問題を直す方法のひとつとしてネットで見つけた方法が「マザーボード上にあるBIOSの設定などの情報(CMOS情報)をリセットする」というものでした。

結果:効果なし

電池を抜いてしばらく(5分ほど)放置することでもCMOSクリアできる,と書いてあるサイトもあったのでその方法でクリアを試みたものの,相変わらずBIOSは起動しませんでした

試した解決策5:MBの取り替え

ここまでくるとMBが壊れている可能性が高いな,と考えてMBを買い替えました。

元のマザーボード

f:id:nigimitama:20190131201854j:plain

もともとG-Tuneで買ったPCで,MBはH81チップセットGIGABYTE GA-H81M-D3V-JP というものがついていました。

買い替えたマザーボード

f:id:nigimitama:20190131201830j:plain

元のマザーボードはH81で,メモリスロットが2つ(8GBのメモリで考えるなら16GBが上限)のものでした。

ただ,普段の作業でメモリ16GBでは不足を感じていたので,元のマザーボードよりちょっと上のグレードのB85チップセット(H81に比べてメモリスロットが2つ増える)のこれに替えました。約1万円。

マザーボード載せ替えの様子

配線がどうなっているのかなどを写真で記録に残しながらコードを抜いていきます…

f:id:nigimitama:20190131202215j:plain

こちらが取り外した古いマザーボード

f:id:nigimitama:20190131202226j:plain

CPUは問題が起こってないと思われるので新しいMBに載せ替えます。

実はCPUに触れるのは数年ぶり(載せ替えるのは初めて)で,CPUの裏面はピンだらけの剣山みたいなやつかと思ってたんですが,今の時代のCPU(写真のはi7-4790k)は剣山タイプじゃないんですね。

f:id:nigimitama:20190131202251j:plain

新しいマザーボードをネジ止めしてCPU,CPUクーラー,メモリ等を載せて配線すれば終了です。(グラボ等はまた後日にして今は最小構成での復旧に努めます) f:id:nigimitama:20190131202302j:plain

結果:BIOSが起動できるようになった

BIOSは正常に起動しました

f:id:nigimitama:20190131202319j:plain

マザーボードのメーカーも世代も変わったのでBIOSというよりUEFIがデフォルトになってますね。

マウスが使えるし日本語も対応しているし非常に使いやすい。

f:id:nigimitama:20190131202331j:plain

Windowsの修復

BIOSが起動できるようになったので,Windows10のインストールプログラムWindowsの修復ができる)を入れたUSBメモリを使ってWindowsの修復を行いました。

それにより,OSが入った旧いHDDのほうでWindowsが起動できるようになりました。

まとめ

BIOSが起動しなくて,CMOSクリアしても起動しないなら,マザーボードを買い換えるのがよさそうです。

マザーボードに1万円くらいかかるし,マザーボードの載せ替えには時間と気力を使うけど ,新品のパソコンを買い直すのに比べれば安く済むはずです。

このパソコンはちょこちょこパーツを替えたりして15万円くらいは投資してるので,それを買い換えることを考えると1~2万円程度で回復できたのはホッとしました。このパソコンにはCPUが死ぬまでは頑張ってもらいたいと考えています。


*1:「MBの電池が尽きる→ BIOSの設定が失われる→ BIOSの初期設定ではOSが入っていないドライブを読み込む→ Windowsが起動しない」というシナリオを想定している