2019-03-25

ABC 122 反省会

C++ Python AtCoder

今回出遅れもあってBまでしか解けなかった。今回からC++でやっていきたいという気持ちがある。

A: Double Helix

最初の回答。

#include <bits/stdc++.h>
using namespace std;
 
int main() {
 char b;
 cin >> b;
 char p;
 if (b == 'A') {
   p = 'T';
 } else if (b == 'T') {
   p = 'A';
 } else if (b == 'C') {
   p = 'G';
 } else {
   p = 'C';
 }
 cout << p << endl;
}

辞書とかきっとあるよなと思って調べたらstd::mapがあったので次からスッと使えるようにしておきたい。

#include <bits/stdc++.h>
using namespace std;

int main() {
  map<char, char> mp{
    {'A', 'T'},
    {'T', 'A'},
    {'C', 'G'},
    {'G', 'C'}
  };

  char b;
  cin >> b;
  cout << mp[b] << endl;
}

B: AT Coder

最初の回答。

#include <bits/stdc++.h>
using namespace std;

int main() {
  string acgt = "ACGT";
  string s;
  cin >> s;
  int longest = 0;
  for (int i = 0; i < s.size(); i++) {
    for (int j = i; j < s.size(); j++) {
      int start = i;
      int end = j;
      bool is_acgt = true;
      for (int k = start; k <= end; k++) {
        bool in_acgt = false;
        for (int l = 0; l < acgt.size(); l++) {
          if (s.at(k) == acgt.at(l)) in_acgt = true;
        }
        if (!in_acgt) is_acgt = false;
      }
      if (is_acgt && longest < (end - start + 1)) {
        longest = end - start + 1;
      }
    }
  }
  cout << longest << endl;
}

多少整理した。文字列の検索を覚えておきたい。

#include <bits/stdc++.h>
using namespace std;

int main() {
  string acgt = "ACGT";
  string s;
  cin >> s;
  int longest = 0;
  for (int i = 0; i < s.size(); i++) {
    for (int j = i; j < s.size(); j++) {
      bool is_acgt = true;
      for (int k = i; k <= j; k++) {
        if (acgt.find(s.at(k)) == string::npos) {
          // 見つからないとstring::nposが返る
          is_acgt = false;
        }
      }
      if (is_acgt) {
        longest = max(longest, j - i + 1);
      }
    }
  }
  cout << longest << endl;
}

C: GeT AC

TLE連発の上ここで時間切れになった。そのまま考え続けて30分後くらいに解けた。出遅れが無かったらあるいは…

まだC++で試行錯誤しながら書くのに慣れていないのでこれはPythonで書いた。

# coding: utf-8
import sys
input = sys.stdin.readline

N, Q = [int(i) for i in input().split()]
S = input()
is_AC = [False] * N
AC_csum = [0] * N
for i in range(N):
  is_AC[i] = S[i:i+2] == "AC"
  AC_csum[i] = AC_csum[i-1] + is_AC[i]

for i in range(Q):
  l, r = [int(i) for i in input().split()]
  ans = AC_csum[r-1] - AC_csum[l-1] + is_AC[l-1] - is_AC[r-1]
  print(ans)

で、その後解説を見ながらC++で書いてみた。考え方は合っていたけど冗長なことをしていたっぽい。

#include <bits/stdc++.h>
using namespace std;

int main() {
  int N, Q;
  string S;
  cin >> N >> Q >> S;
  
  vector<int> at_cnt(N + 1, 0);
  
  for (int i = 0; i < N; i++) {
    at_cnt.at(i + 1) = at_cnt.at(i) + (S.substr(i, 2) == "AC");
  }

  for (int i = 0; i < Q; i++) {
    int l, r;
    cin >> l >> r;
    cout << at_cnt.at(r - 1) - at_cnt.at(l - 1) << endl;
  }

D: We Like AGC

~~今回Dまでたどり着けてないのでこれはまた後で。~~

解説を見ながらC++で書いてみた。

#include <bits/stdc++.h>
using namespace std;

const int MOD = 1e9 + 7;
vector<vector<int>> memo(200, vector<int>(200));
map<string, int> mp;
int mpidx = 0;

bool ok(string last4) {
  for (int i = 0; i < 4; i++) {
    string t = last4;
    if (i >= 1) {
      swap(t.at(i-1), t.at(i));
    }
    if (t.find("AGC") != string::npos) {
      return false;
    }
  }
  return true;
}

int dfs(int cur, string last3, int N) {
  if (mp.count(last3) == 0) {
    mp[last3] = mpidx;
    mpidx++;
  }
  if (memo.at(cur).at(mp[last3]) != 0) {
    return memo.at(cur).at(mp[last3]);
  }
  if (cur == N) {
    return 1;
  }
  int cnt = 0;
  for (int i = 0; i < 4; i++) {
    string acgt = "ACGT";
    char c = acgt.at(i);
    if (ok(last3 + c)) {
      cnt += dfs(cur + 1, last3.substr(1,2) + c, N);
      cnt %= MOD;
    }
  }
  memo.at(cur).at(mp[last3]) = cnt;
  return cnt;
}

int main() {
  int N;
  cin >> N;
  cout << dfs(0, "___", N) << endl;;
}

2019-03-11

RでFizzBuzz

PythonでFizzBuzz

PythonのFizzBuzz ワンライナーってどんなだったかなと思ってちょっと調べた。

for i in range(100):print(i%3//2*"Fizz"+i%5//4*"Buzz"or-~i)

最後の-~って何だと思ったら~はビット反転だった。2の補数、はるか昔に基本情報でやった記憶がうっすらと…

（参考: FizzBuzzを1byteで実装する - Qiita）

AtCoder ABC 121の反省

Python

AtCoder Beginner Contest 121 - AtCoderやったのでまたメモっておく。

A - White Cells

特に悩むことはない。白色のマス目を長方形に並べてカウントすれば良い。

# coding: utf-8

H, W = [int(i) for i in input().split()]
h, w = [int(i) for i in input().split()]

print((H - h) * (W - w))

B - Can you solve this?

特に考える部分は無かった。まずは最初の回答。

# coding: utf-8

N, M, C = [int(i) for i in input().split()]
B = [int(i) for i in input().split()]

correct = 0

for i in range(N):
    A = [int(i) for i in input().split()]
    if sum(i * j for i, j in zip(A, B)) + C > 0:
        collect += 1

print(correct)

ちょっと整理して内包表記で書き直してみたけど別に早くならず可読性だけ落ちてうーん…という感じ。

# coding: utf-8

def read_ints():
    return [int(i) for i in input().split()]

N, M, C = read_ints()
B = read_ints()

correct = sum(
    sum(i * j for i, j in zip(read_ints(), B)) + C > 0 
    for k in range(N)
    )

print(correct)

C - Energy Drink Collector

最初小細工しようとして時間超過してしまったので結局bruteforceっぽい感じで解いてしまった、ソートの速度次第では…と思ってたけど解説読んでも同じこと書いてあったのでこれで良しとする。

件数が多かったのでデータの読み込み方を少し調べた。sys.stdin.readlineを使うと早くなるらしい（cf. Pythonの知っておくと良い細かい処理速度の違い8個 - kumilog.net）。

# coding: utf-8
import sys
input = sys.stdin.readline
 
N, M = [int(i) for i in input().split()]
 
plices = [[int(i) for i in input().split()] for _ in range(N)]
 
plices.sort(key = lambda x:x[0])
 
m = 0
cost = 0
 
for i in range(N):
    if m >= M:
        break
    buy = min(M - m, plices[i][1])
    cost += plices[i][0] * buy
    m += buy
 
print(cost)

D - XOR World

締切に間に合わなかったが、筋は悪くないと思って考え続けた結果終了40分後に解けた。

0〜その数字までを対象に、ある桁に1が何回出現するのか、という点を考えて解いた。スピードは足りてたけど正直汚い。

# coding: utf-8

A, B = [int(i) for i in input().split()]

def count_1(bi, idx):
    bi = bi[::-1]
    cnt = 0
    if idx > 0:
        for i in range(0, idx):
            cnt += int(bi[idx]) * int(bi[i]) * int(2**i)
    if idx < (len(bi) - 1):
        for i in range(idx + 1, len(bi)):
            cnt += int(bi[i]) * int(2**i) / 2
    cnt += int(bi[idx])        
    return(int(cnt))

a = str(bin(max(A-1, 0)))[2:]
b = str(bin(B))[2:]
a = (len(b) - len(a)) * "0" + a # 長さ調整

def count_1s(bi):
    return([count_1(bi, i) for i in range(len(bi))[::-1]])

a_cnt = count_1s(a)
b_cnt = count_1s(b)


print(int("".join([str((i - j) % 2) for i, j in zip(b_cnt, a_cnt)]), 2))

で、解説を見て書き直した。

# coding: utf-8

def F_0n(n): # F(0, n)を求める
    if n % 2 == 0:
        return (n // 2 % 2) ^ n
    else:
        return (n + 1) // 2 % 2

def F(a, b):
    return F_0n(a - 1) ^ F_0n(b)

A, B = [int(i) for i in input().split()]

print(F(A, B))

きちんと排他的論理和の性質から攻めていくべきだった。半端に書き下して半端に法則を見つけてしまい、それをそのままコードにしようとしてしまったのが敗因な気がする。

2019-03-04

AtCoder ABC 120の反省

Python

AtCoder Beginner Contest 120 - AtCoder

やったので回答とメモを残しておく。

A - Faborite Sound

特に注意すべき点はない。

# coding: utf-8
 
a, b, c = [int(i) for i in input().split()]
 
print(min(b // a, c))

B - K-th Common Divisor

$B$の範囲がそれほど広くなかったので、bruteforceで解いてしまった。最初「大きい方から」というのを忘れていて一回ミスってしまった。

提出時は内包表記使ってなかった。慣れておきたい。

# coding: utf-8
 
a, b, k = [int(i) for i in input().split()]
 
divs = [i for i in range(1, a+1) if a % i ==0 and b % i == 0]
 
print(divs[-k])

C - Unification

列の中に異なる色がある限り必ず消滅するので、最後には0の数と1の数の差だけブロックが残る、という点に気づけば簡単。

# coding: utf-8
 
S = list(input())
 
n1 = sum(s == "1" for s in S)
n0 = sum(s == "0" for s in S)
 
print(len(S) - abs(n1 - n0))

D - Decayed Bridges

（ひっくり返せば行けそう…）とまでは気付いたものの時間内には手も足も出なかった。

Union-Find木というものを使う基本的な問題らしく、アリ本（プログラミングコンテストチャレンジブック [第2版]　?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?）ではp81に説明がある（もう少しちゃんと読んでおけばよかった）。

Union-Findの実装はPython:Union-Find木について - あっとのTECH LOGを参考に、グループ数をメモるように少し手を加えた。

# coding: utf-8
 
class UnionFind:
    def __init__(self, n):
        self.par = [i for i in range(n+1)]
        self.rank = [0] * (n+1) # 木の高さ
        self.n = [1] * (n+1)    # グループメンバー数
        
    def find(self, x):
        if self.par[x] == x:
            return x
        else:
            self.par[x] = self.find(self.par[x])
            return self.par[x]
    
    def find_n(self, x):
        return self.n[self.find(x)]
    
    def union(self, x, y):
        x = self.find(x)
        y = self.find(y)
        if self.rank[x] < self.rank[y]:
            self.par[x] = y
            self.n[y] = self.n[x] + self.n[y]
        else:
            self.par[y] = x
            self.n[x] = self.n[x] + self.n[y]
            if self.rank[x] == self.rank[y]:
                self.rank[x] += 1
    
    def is_same(self, x, y):
        return self.find(x) == self.find(y)
    
# -------------------------------------------------------------------

N, M = [int(i) for i in input().split()]

bridges = [tuple(int(i) for i in input().split()) for j in range(M)]

bridges = bridges[::-1]

result = []

now = int(N * (N-1) / 2) #最終不便さ

uf = UnionFind(N)

for i in range(M):
    result.append(now)
    a = bridges[i][0]
    b = bridges[i][1]
    if uf.is_same(a, b): # 既に同じグループ
        pass
    else:
        now -= uf.find_n(a) * uf.find_n(b)
        uf.union(uf.find(a), uf.find(b))

result = result[::-1]

for i in result:
    print(i)

以下Union-Find木に関するメモ。

Union-Find木はグループを管理できるデータ構造で、その名が示すように併合（Union）と検索（Find）の2つのステップが重要となる。

初期化

まず、初期状態（それぞれのノードが1つのメンバーだけを含むグループである状態）を作る。このとき、インデックスとノード番号を一致させておく。

f:id:Rion778:20190304074301p:plain

Union

併合時には、片方のグループの根からもう片方のグループの根へ辺を張る。検索時の計算量を節約するためには、木の高さを低く抑える必要がある。そこで、それぞれの木について高さを保持しておき、併合時には高さを比較して低いものから高いものへ辺を張るようにする（もし高さが同じなら、併合後に高さを+1する）。

f:id:Rion778:20190304074648p:plain

Find

あるノードがどのグループに属するかは、ノードの番号のindexを辿るようにして根を確認すれば良い。ノード番号とindexが同じになったとき、その番号がグループ番号となる。

f:id:Rion778:20190304074756p:plain

また、あるノードのグループが分かったら、ノードが直接根につながるように辺を張り直すと良い。こうすると木の高さを低く保てる。

f:id:Rion778:20190304074905p:plain

2018-11-23

統計検定2級の過去問

間違ったやつと自信をもって解けなかったやつ。

2016年11月問8[1]

問題（要約）

互いに無相関な標準化された確率変数 $X_1$ , $X_2$ , $X_3$ およびそれらの平均 $Y=(X_1+X_2+X_3)/3$ を考える。

$X_1$ と $Y$ の相関係数を求めよ。

回答

相関係数 $r_{X_1Y}$ は次の式で求められる。

$\displaystyle{ r_{X_1Y} = \frac{Cov [X_1, Y] }{\sqrt{V [X_1] V [Y] }} }$

ここで未知なのは共分散 $Cov[X_1, Y$ ]と $V[Y$ ]なので、この2つを求める。

共分散 $Cov[X_1, Y$ ]を求める

共分散は次の式で求められる。

$\displaystyle{ Cov [ X_1, Y ] = E [ X_1Y ] - E [ X_1 ] E [ Y ] }$

$X_1$ は標準化されているため、 $E[X_1$ = 0]であり、つまり $E[X_1Y$ ]を求めれば良い。 $Y$ を展開して整理すると、次のようになる。

$\displaystyle{ E [ X_1Y ] = \frac{1}{3}(E [ X_1^2 ] + E [ X_1X_2 ] + E [ X_1X_3 ]) }$

ここで、次の関係を利用する。各 $X$ は標準化されているという点に再度注意する。

$\begin{aligned} V[X_1] &= E[X_1^2] - E[X_1]^2 \\ &= E[X_1^2] \end{aligned}$

$\begin{aligned} Cov[X_1, X_2] &= E[X_1X_2] - E[X_1]E[X_2] \\ &= E[X_1X_2] \\ &= 0 \end{aligned}$

無相関な2変数の共分散は0である点にも注意。上記の関係を利用すると、

$\begin{align} E[X_1Y] &= \frac{1}{3}(E[X_1^2] + E[X_1X_2] + E[X_1X_3]) \\ &= \frac{1}{3}(V[X_1]) \\ &= \frac{1}{3} \end{align}$

すなわち、共分散は

$\displaystyle{ Cov[X_1, Y] = \frac{1}{3} }$

$V[Y$ ]を求める

3つの変数は無相関なので、和の分散は分散の和になる。

$\begin{align} V[Y] &= V[\frac{X_1 + X_2 + X_3}{3}] \\ &= \frac{1}{9}(V[X_1] + V[X_2] + V[X_3]) \\ &= \frac{1}{3} \end{align}$

答

$\begin{align} r_{X_1Y} &= \frac{Cov[X_1, Y]}{\sqrt{V[X_1]V[Y]}} \\ &= \frac{1/3}{\sqrt{1/3}} \\ &\approx 0.5773503 \end{align}$

2018-10-14

odbcパッケージ経由だとvarcharが255文字に切り詰められる問題

R SQL

ハマるの2回目なのでメモっておく。

PostgreSQLで“test”という名前のデータベースが存在していて、そこに接続している状態を想定する。

con <- odbc::dbConnect(odbc::odbc(), "test")

PostgreSQLのcharacter varying（varchar）は、長さ指定なしで使うと文字列の長さの制限がなくなる（cf. 8.3. 文字型）。

しかし、これをodbcパッケージ経由で取得すると255文字に切り詰められる（https://github.com/r-dbi/odbc/issues/202 で報告されているのと同じ現象？）。

library(dplyr)
query <- "SELECT (repeat('a', 1000))::varchar as str;"
odbc::dbGetQuery(con, glue::glue("CREATE TABLE x AS {query}"))
odbc::dbGetQuery(con, "SELECT str FROM x;") %>% nchar()

## str 
## 255

これは、TEXT型にキャストするととりあえず回避できる。

odbc::dbGetQuery(con, "SELECT CAST(str AS TEXT) FROM x;") %>% nchar()

##  str 
## 1000

dplyrの場合は、as.character()を使う。

tbl(con, "x") %>% mutate(str = as.character(str)) %>% pull() %>% nchar()

## [1] 1000

as.character()はTEXT型へのキャストをするSQLを生成するので、やっていることは同じ。

dbplyr::translate_sql(as.character(x))

## <SQL> CAST("x" AS TEXT)

2018-10-03

「マンガでわかる統計学[因子分析編]」の主成分分析をRでやる

R Statistics

やったのでメモっておく。

「第4章主成分分析」より。

主成分と主成分得点を求める

Step0. データ入力

library(tibble)
ramen <- tribble(
  ~店舗,    ~麺, ~具, ~スープ,
  "二郎",       2, 4, 5,
  "夢田屋",     1, 5, 1,
  "地回",       5, 3, 4,
  "なのはな",   2, 2, 3,
  "花ぶし",     3, 5, 5,
  "昇辰軒",     4, 3, 2,
  "丸蔵らあめん", 4, 4, 3,
  "海楽亭",     1, 2, 1,
  "なるみ家",   3, 3, 2,
  "奏月",       5, 5, 3
)

Step1. 変数ごとに基準化する

library(dplyr)

ramen %>% 
  mutate_if(is.numeric, scale) ->
  ramen_scaled
ramen_scaled

## # A tibble: 10 x 4
##    店舗             麺     具  スープ
##    <chr>         <dbl>  <dbl>   <dbl>
##  1 二郎         -0.671  0.341  1.45  
##  2 夢田屋       -1.34   1.19  -1.31  
##  3 地回          1.34  -0.511  0.759 
##  4 なのはな     -0.671 -1.36   0.0690
##  5 花ぶし        0      1.19   1.45  
##  6 昇辰軒        0.671 -0.511 -0.621 
##  7 丸蔵らあめん  0.671  0.341  0.0690
##  8 海楽亭       -1.34  -1.36  -1.31  
##  9 なるみ家      0     -0.511 -0.621 
## 10 奏月          1.34   1.19   0.0690

Step2. 相関行列を求める

ramen_scaled %>% 
  select(-店舗) %>% 
  cor ->
  ramen_cormat
ramen_cormat

##               麺        具    スープ
## 麺     1.0000000 0.1905002 0.3600411
## 具     0.1905002 1.0000000 0.3004804
## スープ 0.3600411 0.3004804 1.0000000

Step3. 固有値と固有ベクトルを求める

ramen_decomp <- eigen(ramen_cormat)
ramen_decomp

## eigen() decomposition
## $values
## [1] 1.5728539 0.8140083 0.6131378
## 
## $vectors
##            [,1]       [,2]       [,3]
## [1,] -0.5715110  0.6044710  0.5549685
## [2,] -0.5221161 -0.7896069  0.3223595
## [3,] -0.6330639  0.1055260 -0.7668731

※固有ベクトルは定数倍の任意性があるので、テキストと一緒にならない可能性があり、実際なっていない。

Step4

Step3における
- 最大の固有値に対応する固有ベクトル
- 最大から2番目の固有値に対応する固有ベクトル
- をもとに散布図を書く

library(ggplot2)
library(ggrepel)
target <- 1:2
target_vec <- ramen_decomp$vectors[, target]
as.data.frame(target_vec) %>% 
  mutate(item = c("麺", "具", "スープ")) %>% 
  ggplot(aes(x = -V1, y = -V2)) + # テキストに合わせるために反転
  geom_point() +
  geom_text_repel(aes(label = item), family = "IPAexGothic") +
  lims(x = c(0, 1), y = c(-1, 1)) +
  geom_vline(xintercept = 0, col = "gray") +
  geom_hline(yintercept = 0, col = "gray") +
  theme_minimal()

f:id:Rion778:20181003211033p:plain

Step5. 第一主成分と第二主成分を理解する。

第一主成分: 最大の固有値に対応する固有ベクトル。
- 上記例では-0.571511, -0.5221161, -0.6330639
- 頭から順に麺、具、スープに対応している。
第二主成分: 二番目の固有値に対応する固有ベクトル。
- 上記例では0.604471, -0.7896069, 0.105526
- 同様に、n番目の固有値に対応する固有ベクトルが第n主成分となる。

Step6 各個体の第一主成分と第二主成分における座標、つまり各個体の第一主成分得点と第二主成分得点を求める。

ramen_scaled %>% 
  rowwise() %>% 
  mutate(
    score1 = -sum(c(麺, 具, スープ) * ramen_decomp$vectors[, 1]),
    score2 = -sum(c(麺, 具, スープ) * ramen_decomp$vectors[, 2])
  )

## Source: local data frame [10 x 6]
## Groups: <by row>
## 
## # A tibble: 10 x 6
##    店舗             麺     具  スープ score1 score2
##    <chr>         <dbl>  <dbl>   <dbl>  <dbl>  <dbl>
##  1 二郎         -0.671  0.341  1.45    0.712  0.522
##  2 夢田屋       -1.34   1.19  -1.31   -0.974  1.89 
##  3 地回          1.34  -0.511  0.759   0.980 -1.29 
##  4 なのはな     -0.671 -1.36   0.0690 -1.05  -0.678
##  5 花ぶし        0      1.19   1.45    1.54   0.789
##  6 昇辰軒        0.671 -0.511 -0.621  -0.277 -0.744
##  7 丸蔵らあめん  0.671  0.341  0.0690  0.605 -0.144
##  8 海楽亭       -1.34  -1.36  -1.31   -2.31  -0.127
##  9 なるみ家      0     -0.511 -0.621  -0.660 -0.338
## 10 奏月          1.34   1.19   0.0690  1.43   0.124

Step7 第一主成分と第二主成分に基づいてプロットする。

ramen_scaled %>% 
  rowwise() %>% 
  mutate( # テキストに合わせるために-1を掛けて反転させる
    score1 = -sum(c(麺, 具, スープ) * ramen_decomp$vectors[, 1]),
    score2 = -sum(c(麺, 具, スープ) * ramen_decomp$vectors[, 2])
  ) %>% 
  ggplot(aes(x = score1, y = score2)) +
  geom_point() +
  geom_text_repel(aes(label = 店舗), family = "IPAexGothic") +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  theme_minimal()

f:id:Rion778:20181003211155p:plain

分析結果の精度を確認する

寄与率を計算。

ramen_decomp$values / 3

## [1] 0.5242846 0.2713361 0.2043793

累積寄与率を計算

cumsum(ramen_decomp$values) / 3 * 100

## [1]  52.42846  79.56207 100.00000

同じことをRでやる

Rで、というかprcompで。

ramen %>%
  select(-店舗) %>% 
  prcomp(scale = TRUE) -> result
result

## Standard deviations (1, .., p=3):
## [1] 1.2541347 0.9022241 0.7830312
## 
## Rotation (n x k) = (3 x 3):
##              PC1        PC2        PC3
## 麺     0.5715110 -0.6044710  0.5549685
## 具     0.5221161  0.7896069  0.3223595
## スープ 0.6330639 -0.1055260 -0.7668731

summary(result)

## Importance of components:
##                           PC1    PC2    PC3
## Standard deviation     1.2541 0.9022 0.7830
## Proportion of Variance 0.5243 0.2713 0.2044
## Cumulative Proportion  0.5243 0.7956 1.0000

先程の計算と一致していることを確認（※主成分に定数倍の差はある）。

標準偏差は固有値の平方根に一致する。

sqrt(ramen_decomp$values)

## [1] 1.2541347 0.9022241 0.7830312

ramen_decomp

## eigen() decomposition
## $values
## [1] 1.5728539 0.8140083 0.6131378
## 
## $vectors
##            [,1]       [,2]       [,3]
## [1,] -0.5715110  0.6044710  0.5549685
## [2,] -0.5221161 -0.7896069  0.3223595
## [3,] -0.6330639  0.1055260 -0.7668731

結果はbiplotで可視化できる。

biplot(result, family = "IPAexGothic", xlabs = ramen$店舗)

f:id:Rion778:20181003211949p:plain

ところで、ggplotを使いたいと思ってggfortify、ggbiplot、factoextraの3つのパッケージを試してみたのだが、現状日本語が問題なく使えるのはggfortifyだけのようだった。

ggfortifyはパッケージロード後にautoplotを通じてプロットをするが、引数を確認するには?ggfortify::ggbiplotを見る必要があった。

library(ggfortify)
autoplot(
  result2, 
  loadings = TRUE, 
  loadings.label = TRUE, 
  loadings.label.family = "IPAexGothic",
  loadings.label.repel = TRUE,
  label = TRUE,
  label.family = "IPAexGothic",
  label.repel = TRUE,
  scale = 0
) +
  theme_bw()

f:id:Rion778:20181003212928p:plain

2018-09-02

農薬の名前はポアソン分布に従っている

気がするのでちょっと見て欲しい。

まず、rvestを使ってFAMICから農薬登録情報のcsvをダウンロードして保存する。URLがちょいちょい変わるのでxpathとかで抽出したほうが良い。

library(rvest)
library(dplyr)
url <- "http://www.acis.famic.go.jp/ddata/index2.htm"
read_html(url) %>%
  html_node(xpath = '//*[@id="mainArea"]/ul/li[1]/ul/li/a') %>% 
  html_attr(name = "href") %>% 
  url_absolute(base = url) %>% 
  download.file("kihon.zip")
kihon <- unzip("kihon.zip") %>% 
  readr::read_csv(locale = readr::locale(encoding = "CP932")) %>%
  mutate_if(is.character, stringi::stri_trans_nfkc)

次に農薬の名称を文字数に変換する。比較のために他の文字列情報もいくつか一緒に変換しておく。

kihon %>% 
  select_if(is.character) %>% 
  select(-c(濃度, 用途, 剤型名, 登録年月日, 登録の有効期限)) %>%
  tidyr::gather() %>%
  unique() %>% 
  mutate(n = nchar(value)) -> kihon_mod

次にポアソン分布の確率質量を各変数の平均値から計算しておく。

# 平均値
lambda <- kihon_mod %>% 
  group_by(key) %>% 
  summarise(lambda = mean(n, na.rm = TRUE))
# 確率質量
pois <- data.frame(n = 1:50, key = rep(lambda$key, rep(50, nrow(lambda)))) %>% 
  left_join(lambda, by = "key") %>% 
  mutate(theory = dpois(n, lambda))

そしてプロットする。

library(ggplot2)
library(gghighlight)
kihon_mod  %>% 
  ggplot(aes(x = n)) +
  geom_bar(aes(y = ..prop..)) +
  geom_line(data = pois, 
            aes(x = n, y = theory)) +
  gghighlight(key == "農薬の名称") +
  facet_grid(key ~ ., scales = "free_y") +
  lims(x = c(0, 30)) +
  theme_bw(base_family = "IPAexGothic")

f:id:Rion778:20180902224525p:plain

農薬の名称はランダムに決まっている…？

ところで最初に雑に「農薬の名前」と言ってしまったが、農薬の名前には「農薬の名称（＝商品名）」の他に「種類名」と「化学名」があってややこしいのだが、今回ポアソン分布っぽい分布をしているのは商品名だ。

ちなみに今回プロット対象にした項目の実質的な意味は次の通り。

正式名称...販売しているメーカーの名前。
総使用回数における有効成分...使用回数のカウントの際に用いる有効成分名。なんと実際の有効成分名とイコールではなく厄介（分かりやすい例で言うとこちらでは光学異性体をまとめて1つの扱いにしているものがある）。
農薬の種類...種類名。
農薬の名称...商品名。
有効成分...有効成分名。

なぜもっと単純にならないのか…。

2018-05-11

歯を抜いたので状況を記録しておく

Diary

親知らずを抜いた。3週間ぶり2本目。次はないので、経過を記録しておきたい。随時追記していく所存。

抜歯にいたる経緯

昨年11月頃、歯が痛くなってしまったので歯医者に行った。昔詰め物をしたところがあったのだが、詰め物が取れてしまい、そのままほっといたら虫歯になってしまったのだ。虫歯の治療は割と痛かったし、セラミックの詰め物(というか被せ物)にしたら財布へのダメージも甚大だったし、歯医者には定期的に行くべきだ。

それはともかく、治療にあたってレントゲンを撮影した。すると、下の両奥歯の奥、明らかに角度を間違えた歯が2本ある。親知らずだ。横向きに生えている親知らずは抜くのが大変だとか聞いたことがあるが、抜かないといけないのだろうか。

「んーまあ左はちょっと頭出てるし、右もほぼ出てるようなもんだし、抜けるでしょう」

抜くことは確定らしい。

横向きに圧力がかかって歯並びに影響しているし、ちょっと出てるから虫歯のリスクも高く、さっさと抜いておくべきとのこと。

何か軽い感じで言ってるし、大した事ないのではないか、その時はそう思っていた。

それからしばらくして虫歯の治療が完了し、親知らずの抜歯の話になった。

「あっでは大学病院に紹介状書きますのでー」

そうだ。横向きの親知らずがそんなに簡単に抜けるわけがない。私は大学病院送りになった。

1本目の抜歯

予約した日に大学病院へ行くと、学生らしき若者に体温と血圧と脈拍を取られる。大学病院っぽさだ。誰だって最初は初心者なのだから仕方がないのだが、一向に脈が取れず、すごく不安感を煽られる。体温は37.8℃あった。どうも極度に緊張したり不安になったりすると発熱する体質らしいと気付いたのはこの時だった。

その日はレントゲン撮影やらなんやらを行った。レントゲン装置が顔の周りをぐるっと回る間、何故か装置からエリーゼのためにが流れる。誰だ無駄な機能を付けたのは。不意打ちはやめろ。そしてその日は日程調整をしただけで抜歯はしなかった。壮絶な決意で行ったのに。両方一度に抜くのではなく、一本抜いて落ち着いたらもう一本抜く、という方針で行くらしい。

日を改めて1本目の抜歯を行った。その日測定してもらった体温は37.6℃だった。親知らずの抜歯はちょっとした手術になるので、リスクの説明とか同意書へのサインとか諸々の儀式があり、不安感を煽られる。

「抜いてる時は最初麻酔でチクッとするだけで痛くないですが、その後は確実に痛いです。3日目くらいが腫れのピークです。」

救いはないのか。

そして手術が始まる。若干想定外だったのは「最初麻酔でチクッと」が1回ではなく3回くらい点だ。割とブスブス刺された気がする。さほど痛くはなかったが。麻酔が効いてしまうと確かに痛みはない。顔には覆いがされているので、後は音で想像するしかない。「水が出ますよ―」と何度か言われたが、正直水が出ているのかどうか、何のために水が出ているのかも良くわからなかった（多分削りカス的なものをアレしたりするためだろうけど想像するのはやめておいた）。

どれくらいの時間がたっただろうか。後で確認したら30分くらいだった。

「抜けましたからねー」

抜けたらしい。左の親知らずだったものは3分割されていた。

「持って帰る？」

いらない。

抜歯後の経過

痛み止めに処方してもらったロキソニンが効いてる間はさほど痛くはない（無痛ではないが、強い不快感があるほどでもない）が、ロキソニンが切れるとすぐ痛くなるという状態が2〜3日は続いた。4日目くらいに血餅が取れてビビったが、これくらい時間が経過していれば大丈夫らしい。

1週間後に再度大学病院に行って抜糸した。抜歯と違って秒単位で終わったが、意外と痛かった。もう一度あると思うと若干不安。この時2本目の抜歯の予約をした。（当然抜くよな）という圧を感じた。選択の余地は無かった。

2週間くらいすると違和感も大分なくなった。ただ穴はあいているので食べ物が入ったりして若干不快感はある。この穴は骨が再生して埋まるらしいが、数ヶ月くらいかかるらしい。

この間、割と大変だったような気もするのだが、喉元過ぎれば何とやらで、細かい部分は忘れてしまった。次回はちゃんと記録したいと思った。

2本目の抜歯

そして2本目の抜歯だ。前回の反省を活かしてこまめに経過を記録していきたい。

抜歯前日まで

前回抜歯時に調査した結果として、抜歯後の早い段階で血餅（かさぶた）が取れてしまうとドライソケットと呼ばれる状態になり、激痛に襲われるという情報を得ていた。しかも下の親知らずはドライソケットのリスクが高いらしい。そこで、少しでも可能性を減らすべく、食糧としてカロリーメイトのゼリータイプを買い込んだ。昔喉をやられたときに1週間くらいこれだけで生きたことがある。個人的にはりんご味が好きだ。

抜歯当日（2018/05/11）

2本目だろうが何本目だろうが嫌なものは嫌なので、めっちゃ嫌だなと思いながら大学病院へ向かう。

体温と血圧と脈拍を測定してもらう。やはり学生、スムーズには行かない。こういう実習も大切と思いつつも不安感が高まる。体温は37.3℃だった。

麻酔をして処置が始まる。前回30分くらいだったので、1800まで数えれば終るはずだと思い、数を数えることに集中する。やたら「破骨鉗子」という単語が聞こえる。何か破片を引っ張り出してるのだろうかと思うが想像するのはやめておく。

「抜けましたからねー」

と言われたのは320くらいまで数えた時だった。きちんと30分経過していたし、30分くらい経ってるよなという感じもあった。人は緊張するとまともに数字を数えられなくなるのだという知見を得た。

右の親知らずだったものは8つくらいの破片に分割されていた。なにか掴んで引っ張り出してるっぽかったのはこれか。

「持って帰る？」

いらない。

ロキソニンは家に在庫があったので、今回は化膿止めだけが処方された。あれ、前回はうがい薬的なものが出た気が…。

不安感を抱えつつその日はそのまま出社し、そっと仕事をした（家より会社のが近いのだ）。

やはりロキソニンが効いてるうちは何ということは無いのだが、薬が切れてくると痛い。すごく痛いという訳ではないのだが、ロキソニンの薬効は5時間くらいしか持たないので、夜中に1度は目が覚めてしまうのが辛いところ。

抜歯後の食事は全てカロリーメイトのゼリーで済ませた。

抜歯翌日（2018/05/12）

今朝は特に痛みがないと思ったが、そういえば夜中に一度痛みで目が覚めてロキソニンを飲んだのだった。時間を確認していなかったが明け方だったのかもしれない。

昨晩寝るまでの間は若干の出血を感じていたが、もう完全に止まっている。

抜歯した場所に触れないように注意しながら歯磨きをする。こういうときはタフトブラシが便利だ。前回歯を抜いた時はネオステリングリーンという洗口剤を処方されたのだが、今回はなぜか処方されなかったのでコンクールFといういつも使っている洗口剤でそっとうがいをした。この洗口剤も抜歯後に使用される場合があるもので、アルコールが入っていないので刺激が少ないが口の中はすっきりする。長期間の使用で着色しやすいという副作用があるらしい。

朝はそれほどでもなかったがやはり薬効が切れると多少痛みが出てきて、耐えられないほどではないものの不快感があるので今日は11時と19時にロキソニンを服用した。少し頬に腫れを感じる。

痛み止めさえ効いていれば全く痛みがなく、昨日抜歯したことすら忘れる。普通に食事できそうな気もしたが、固いものを食べるとうっかり抜歯した側で噛んでしまいそうなのでうどんやグラタンなどの柔らかいものを注意しながら食べた。

抜歯翌日というのは割と大丈夫なのだ。前回もそうだった気がする。

抜歯2日後（2018/05/13）

また痛みで目が覚めて明け方にロキソニンを飲んだ。昨日時刻を確認し忘れたので時計を見たところ3:50だった。10分くらいで痛みが引いて寝た気がする。

今日は朝から若干痛みというか不快感がある。前回抜歯したときはあまり気にならなかったが、今回は外観から分かる程度に腫れてきた。

今日は9時と19時にロキソニンを飲んだ。ロキソニンの影響なのか腹の具合が悪い。飲めば腹が痛いし飲まねば歯が痛いし割と辛い。やはり2〜3日後がピークなのか。

食事は薬さえ効いていれば結構普通にできるようになっているが、念のため長めに茹でたうどんを連投している。

抜歯3日後（2018/05/14）

やはり夜中に目が覚める。今回は2時頃だった。寝ている間にロキソニンを摂取するソリューションが求められる。

朝の違和感がないので昨日よりは多少マシ、一昨日と同等といった状況。頬はまだ腫れているが、触った時に軽い痛みを感じるエリアは狭まっており、多少は回復しているものと思われる。昨日がピークだったと思いたい。

今日は10時と15時にロキソニンを飲んだ。薬が切れるとやはり痛い。朝痛んできたのでロキソニン飲もうと思ってたら会議が思ったより長引いて会議どころではなかった。今考えると別に会議途中だろうがサッと薬飲めばよかった。

抜歯4日後（2018/05/15）

ロキソニンが切れると目が覚めるのであれば、目が覚める時間にロキソニンを飲めば早朝に起床できるのではないか。そう思って昨夜12時頃にロキソニンを飲んだら6時に目が覚めた。抜歯は生活習慣の改善に使えるのではないか。しかし痛いので結局7時にはロキソニンを飲んだ。

朝若干雑に歯を磨いたら少し出血してしまいビビる。まだ注意しないと。

昼を過ぎても痛みが出なかったのでコイツは行けるかと思ったけど3時を過ぎたらまた痛みだした。なんやかんやで前回も1週間くらいは痛かった気がする。

抜歯5日後（2018/05/16）

結局今日も寝る前、朝、3時頃と3回ロキソニンを飲んだ。激痛ではないが、ロキソニンあるなら飲むわという程度には痛い。薬さえ飲んでしまえば無痛。

明日抜糸する予定なのだが大丈夫なのか。よく考えると前回抜糸したときも（まだちょっと痛いけど大丈夫なのか…）と思っていた気はする。前回抜糸したとき、抜糸自体が割と痛かった気がするのであまり行きたくないのだが、虫歯に端を発する長かった歯医者通いも明日で一旦区切りと思うと感慨深いものがあるような気もするし、ないような気もする。

抜歯6日後（2018/05/17）

今朝も起きると歯(のあった場所)が痛い。寝ている間に歯ぎしりをしているような気もする。

抜糸もめっちゃ嫌なのでめっちゃ嫌だなと思いながら大学病院に向かった。若干早く着いてしまったところ若干早く案内をしてくれた。

処置を待つ間、隣の席の人の話が聞こえてきたが、月1回佐渡から通ってるとかなんとかですごく大変そうな雰囲気を感じた。

肝心の抜糸は「前のオペが長引いている」とかいう理由で、抜歯したのとは別の先生がやることになった。前回のときもそうだったのだが、抜歯した人が抜糸まで面倒をみる的な決まりになっているのだろうか。

抜糸はやっぱりちょっと痛かったが、前回よりは痛くなかった気がする。でも処置後にうがいしたら思ったより出血してて吐いた水が赤く染まっており若干引いた。わざとそうしてあるのかもしれないが、歯医者の診療台についてるうがい用の流し、血の色がすごく目立つ配色だから出血してるとビビる。処置の時間は1分もかかってないと思う。

まだ痛みが続いているので、ロキソニンを10回分処方してもらった。

大学病院で行う処置はこれで完了したということで、もともと紹介状を出してくれた行きつけの歯医者に治療終わりました連絡をするのでちょっとお値段かかりますという話だった。多少余分にかかった気はする。連絡も業務なのでしっかり金とってって欲しい。

ともかくこれで歯医者通いは終わりなので、何事もなければ次回は定期検診だ。めでたい。

抜歯7日後（2018/05/18）

一週間経った。

糸を抜いて大分良くなったなという気持ちではあったものの、今日もロキソニンは7時、12時、20時と3回飲んだ。

ただ、歯を磨いたら血が出るとか痛いとかそういうことはもうない。昨日くらいまでは歯肉がフワフワしていたような気がするが、大分しっかりしてきた。

薬飲んでれば痛くないし薬が切れると痛いという代わり映えのしない日々が続いていてボチボチ飽きてきた感じはある。

抜歯8日後（2018/05/19）

朝起きた時に割と痛かった。今日は結局朝と夜の2回ロキソニンを飲んだ。3回が2回になったので大分良くなっている印象がある。

もう食べ物は完全に自由に食べられるが、右も左も親知らずを抜いた穴が開いている状態なので、何をどう注意しても穴に食べ物が入ってしまう。そして後に抜いた方に食べ物が入ると少し痛む。前回もこんなに長い間痛みがあっただろうかと思ったが、前回は抜歯10日後くらいに実家に帰ったら割と辛いカレーが夕食に出てきて割と痛かったような気がした。まだ8日目なのでこんなもんだろう。

ところで左側の親知らずを抜歯してから今日で丁度1ヶ月だった。左側の方はというと痛みはもう完全に無いが、若干突っ張っているような違和感はある。

抜歯9日後（2018/05/20）

朝起きてもほとんど痛くなかった。朝食を普通に食べても大丈夫だったのでコイツはいけるだろと思っていたが昼食を食べた当たりで若干痛みのような違和感が出てしまい結局ロキソニンを飲んだ。今日は1回で済んだので、ぼちぼち薬も卒業できそう。

ところで右の抜歯跡に比べると左の抜歯跡は大分食べ物が入り込む事がなくなっている様に思う。1ヶ月も経つし埋まってきているのだろうか。

抜歯10日後（2018/05/21）

起床後少しの間痛いような気がしたが、しばらくそっとしていたら治まった。どうも寝ている間の歯ぎしりとかが原因のような気がする。

もうそんな痛くないしいけるやろとか油断してたら昼食後に痛くなってしまった。食べ物が入ると痛い感じがする。耐えられないほどではなかったが、別に我慢しても良いことも無いので結局痛み止めを飲んだ。

抜歯11日後（2018/05/22）

基本的には痛くないのだが、ふとした瞬間から痛みだすことがある。気になって舌で触ってしまうのが良くない気がする。一旦痛くなるとしばらく痛いので結局今日も1回ロキソニンを飲んだ。薬を飲んでしまえばその後痛むことはない。

抜歯12日後（2018/05/23）

変化が無いわけではないが、ゆっくりと回復していくだけなのでいい加減書くことが無い。とりあえず朝起きたら痛いとかそういうことはもう無い。

しかしやはりふと痛みだす時があって、今日も1回だけロキソニンを飲んだ。

抜歯13日後（2018/05/24）

電動歯ブラシが抜歯した跡に当たるとさすがに少し痛い、というのが昨日までの状況だったが、今朝はそれもあまりなかった。

抜歯跡に食べ物が入ると多少痛いが、すぐにおさまるのでやっと痛み止めを飲まずに1日過ごせた。2週間も経過すれば平気になるようだ。

抜歯14日後（2018/05/25）

もはや食べ物が抜歯跡に入らなければ全く痛くないのだが、抜歯した方で普通に噛むという行為をしばらくしていなかったのでうっかり頬の内側を噛んでしまって少し痛かった。

抜歯15日後（2018/05/26）

もう食べ物が抜歯跡に入っても痛くない。そしてこうなってしまうともう書くことがない。今となっては食べ物が詰まりやすい穴が奥歯の奥に空いているだけだ。前述のようにこの穴が埋まるのには数ヶ月とかかかるらしいので、記録はこのへんで終わりにしようと思う。

2018/06/18

左の奥歯を抜いて約70日が経過したが、だいぶ穴が歯肉で埋まってきているようで、舌先で触ると分かる程度になってきた（最初はまたなにか食べ物が詰まっているのかと思った）。右側の抜歯跡が若干知覚過敏のようになっているらしく、冷たいものを飲み食いするとしみるのが気になる。

2018-02-20

Rule of Three

Statistics

調べ物をしていて、このような法則があることをふと思い出した。

病害虫の発生率が $p$ %以下であるということを95%の信頼度で言うためには、3/ $p$ 個の対象を調べて病害虫が存在しないことを示せば良い。

この法則はRule of Threeあるいは3の法則と呼ばれている*1。

Rule of Threeは医薬品の稀な副作用を調査するという文脈の上で、次のように言われることもある。おそらくこちらの表現の方が有名だろう。

100人調べて1度も観測されなくても別の100人中3人に事象が生起する可能性があり、100 万人調べて0でもなお他の100万人中3人に生起する可能性は否定できない

導出

なんとなく「ポアソン分布から導出できる」みたいなアバウトな覚え方をしていて、いざやってみたら意外と手こずってしまったのでメモしておく。

まず前提条件を整理すると次のようになる。

事象が起こる確率 $p$ はすごく小さい(と予想している)
試行回数 $n$ は大きい
$n$ 回の調査で1回も事象は発生しなかった

$p$ は二項分布に従うと考えられるが、 $p$ が大きく $n$ が小さい二項分布は $np=\lambda$ をパラメータとするポアソン分布で近似できる。ポアソン分布において事象が $k$ 回起こる確率は

$P(k) = \frac{\lambda^k e^{-\lambda}}{k!}$

だが、今回 $k=0$ なので、

$P(0) = e^{-\lambda}$

となる。ここで、 $\lambda$ の上限を考えてみよう。 $\lambda$ が大きくなるほど $P$ は小さくなるだけで $\lambda$ はいくらでも大きくできることはできる。ただ、あまり小さな $P$ を想定するのは不自然なので、 $P=0.05$ に対応する $\lambda$ あたりを上限としよう。

わざとらしい感じになってしまったが要するに $\lambda$ の片側95%信頼区間を求めるということだ。

$0.05 = e^{-\lambda}$

を $\lambda$ について計算すると、

$\lambda = -\ln{0.05} \approx 3$

となる。要するに、n回の試行で事象の発生回数が0だったとしても、 $\lambda$ =n回の試行における事象の発生回数の期待値の95%信頼区間はほぼ0〜3ということになる。

これが

100 万人調べて0でもなお他の100万人中3人に生起する可能性は否定できない

に対応する。ここでさらに $np=\lambda$ として $n$ について整理すると

$n = \frac{3}{p}$

となり、これが