axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

ニュートン法で平方根の値を求める

きっかけ

電卓のルートってどうやって求めるんだっけ? → たぶんニュートン法だろうなぁ。ということで、Rで簡単に実装して計算してみました。

実装

ニュートン法

fnewton <- function(xn=1,p=2,a=2){
  return ( xn - ( xn^p - a )/( p * ((xn)^(p-1)) ) )
}

反復計算からグラフ作成

solve_with_fnewton <- function(a, p = 2, n, xn = 1 ){
  xs <- c(xn,xn+1:n)
  for(i in xs){
    xs[i+1] <- fnewton(xs[i],p,a)
  }
  last_xs <- xs[n]
  
  plot(xs,type="o")
  title(main=bquote(sqrt(.(a), .(p)) %~~% .(last_xs) ),sub = "")
  abline(h=a^(1/p) ,col="red")
  print( xs )
}

関数呼び出し

solve_with_fnewton(a = 2,n = 10,p = 2)

結果

f:id:axjack:20190430235037p:plain

f:id:axjack:20190430235121p:plain

感想

6回目の反復あたりでほぼ答えが電卓と一致して驚いた。

ソフトバンクからOCNモバイルONEにMNPで乗り換える

はじめに

平成もそろそろ終わりそうな今日この頃、ソフトバンクから格安SIMなOCNモバイルONEにMNPしながら乗り換える*1ことにした。端末も同時に切り替える。

この記事を書いている現在(3/19)はOCNモバイルONEからSIMカードが届くのをまだかまだかと待っている最中である。 OCNモバイルONEからSIMカードが届いたので追記する。

なお、本記事は手続きが進み次第随時追記していく。

手続き

諸元

手続き経過

日時 イベント
3月16日(土) Apple 公式サイトにてiPhone 8を購入する
3月18日(月) ソフトバンクMNP予約の電話をかける
3月18日(月) OCNモバイルONEに申し込む
3月18日(月) iPhone 8 が家に届く
3月25日(月) SIMカードを届けに来た旨の不在通知を佐川急便より受領
3月26日(火) SIMカードを佐川急便より受領
3月26日(火) 構成プロファイルをインストール
3月26日(火) MNP開通手続き

SIMカードが届いてから

構成プロファイルのインストール

構成プロファイルのダウンロードはこちらの案内に従い、http://s.ocn.jp/ltecfg へアクセス、構成ファイルをダウンロードした。この時、端末のSIMカードソフトバンクのままの状態。

MNP開通手続き

MNP開通手続きのご案内に従い、http://s.ocn.jp/mnp01 へアクセス、申し込み時に発行されたocnのメールアドレスとパスワードを入力し、MNP開通手続きを実行した。この時、端末のSIMカードソフトバンクのままの状態。

実行後、iPhoneアンテナマークの右にある回線事業者名:Softbankが消えた。(この時点でMNP開通手続きが完了した、と思われる。)

SIMカード入れ替え

iPhone8はナノSIMなのでSIMカードをナノサイズに切り取り、iPhone8からソフトバンクSIMカードを取り出して交換した。

すると、回線事業者名がdocomoに!

OCN関連のアプリをインストール

  • OCNモバイルONE
  • OCNでんわ

の2つのアプリをインストールした。

各種疎通テスト

以下を問題なく実施できた。

  • wifiを切断してOCNモバイルONEからインターネットに接続する
  • SMSメッセージを受信する
  • iMessageを受信する
  • iMessageを送信する
  • OCNでんわから電話をかける
  • 他社携帯から電話を受ける

おわりに

簡単だったこと

手続きがものすごく簡単だった。SIMカードが届くまで1週間掛かったものの、スマホは何事もなく使えたしなんら問題なく。

難しかったこと

SIMカードの取り出しが難しかった。慣れてないとコツをつかむまでに苦労しますね。まぁ、勢いよくぐいっと押せば蓋が開くのですが・・・。

*1:乗り換えた理由:[1]ソフトバンク、10年ぐらい使ったけど月額料金がやっぱり高い[2]ソフトバンクショップ店員の態度が気に食わなかったのが最大のきっかけとなった[3]OCNモバイルONEのトップ3かけ放題というサービスを使えば電話代もあまり気にしなくても良さそう

固定長ファイルをテキストエディタで分割する

はじめに

3月もあと10日ほどで終わりですね。

本日は、COBOLの時代からひっそりと伝わる固定長ファイル分割術を、記憶が薄れる前にここに書き留めておこうという趣旨の記事です。

で、Linux環境だったら伝家の宝刀ddコマンドを駆使するわけですが、Windows環境となるとそのような手軽なツールは準備されていませんね。ということで、テキストエディタでグリッとこなす技を記載します。

分割ってどういうこと?

こんな感じの一行ファイルを、

1HHHHHHHHHHHHHHHHHHHHHHHHHHHHH21DDDDDDDDDDDDDDDDDDDDDDDDDDDD22DDDDDDDDDDDDDDDDDDDDDDDDDDDD23DDDDDDDDDDDDDDDDDDDDDDDDDDDD24DDDDDDDDDDDDDDDDDDDDDDDDDDDD25DDDDDDDDDDDDDDDDDDDDDDDDDDDD26DDDDDDDDDDDDDDDDDDDDDDDDDDDD27DDDDDDDDDDDDDDDDDDDDDDDDDDDD28DDDDDDDDDDDDDDDDDDDDDDDDDDDD29DDDDDDDDDDDDDDDDDDDDDDDDDDDD8TTTTTTTTTTTTTTTTTTTTTTTTTTTTT9EEEEEEEEEEEEEEEEEEEEEEEEEEEEE

このように区切ることを、

1HHHHHHHHHHHHHHHHHHHHHHHHHHHHH
21DDDDDDDDDDDDDDDDDDDDDDDDDDDD
22DDDDDDDDDDDDDDDDDDDDDDDDDDDD
23DDDDDDDDDDDDDDDDDDDDDDDDDDDD
24DDDDDDDDDDDDDDDDDDDDDDDDDDDD
25DDDDDDDDDDDDDDDDDDDDDDDDDDDD
26DDDDDDDDDDDDDDDDDDDDDDDDDDDD
27DDDDDDDDDDDDDDDDDDDDDDDDDDDD
28DDDDDDDDDDDDDDDDDDDDDDDDDDDD
29DDDDDDDDDDDDDDDDDDDDDDDDDDDD
8TTTTTTTTTTTTTTTTTTTTTTTTTTTTT
9EEEEEEEEEEEEEEEEEEEEEEEEEEEEE

分割と称することにします。

想定するデータの形式

ありがちな、

  • ヘッダレコード
  • データレコード
  • トレーラレコード
  • エンドレコード

の全銀フォーマット的な形式とします。何がありがちなのかよくわからな人は全銀フォーマット 固定長でggってください。

文字種は半角アルファベット・半角数字、ファイルに改行は含まれず、各レコードは30byteとします。

やりかた

ヘッダ→エンド→トレイラー→データの順に分割します。

ヘッダレコード

ヘッダレコードはファイルの先頭に陣取るレコードなので、

^.{30}

正規表現を使って検索すれば、頭30byteが取得できます。

エンドレコード

エンドレコードはファイルの末尾に陣取るレコードなので、

.{30}$

正規表現を使って検索すれば、お尻30byteが取得できます。

トレイラーレコード

エンドレコードを切り取ってからもう一度エンドレコードと同じやり方で検索すると、トレーラレコードが取得できます。

データレコード

ヘッダ・エンド・トレーラを前述の手順で切り出してしまえば、 .{30} で引っかかるのがデータレコード、となります。

レコード種別に関係なく30byteずつ分割する

(.{30})\1\nと置換すると、30byte引っかかるごとに改行を挟むことになるので、結果的に分割できます。なお改行コードが気になる方は文字コード環境に合わせて\nなり\r\nなり¥r¥nなりに変更してください。

おわりに

機械学習人工知能の波が押し寄せようと来まいと、データフォーマットの変換がなくなることはないでしょう。固定長・CSVXML・JSOなど、どんな何が来てもグリっと変換できるようにしておけば、平成の次の時代でも困ることはないはず(たぶん)。

唐突にテキストファイルだけ渡された時、しれっとフォーマット変換しなければいけない、だが手元にテキストエディタしかないのだが?な時にこの記事が役に立ったら幸いです。

【追記】分割されたファイルから改行を消す

\nをブランク(≠半角スペース)で置換すれば元に戻ります。

ggplot, geom_point, facet_gridの練習

はじめに

ggplotやqplotを日々練習しています。ggplot, geom_point, facet_gridを組み合わせると、

  • グラフ内側のx軸(量的)
  • グラフ内側のy軸(量的)
  • グラフ外側の横側(質的)
  • グラフ外側の上側(質的)
  • 点の色(質的・量的)

のように5変数くらい同時に、割と分かりやすく表示することができる。分割が煩わしいなら、facet_gridの項を消せば1つの散布図色付きのように表示することができる。

コード

例1

library(dplyr)
library(ggplot2)
ggplot(mpg) + geom_point(aes(cty,displ,colour=fl)) + facet_grid(year~cyl)

f:id:axjack:20190228230847p:plain

例2

library(ISLR)
a1 <- Auto %>% mutate(year_ = ifelse(year <= 80, 70, 80))
ggplot(a1) + geom_point(aes(displacement,horsepower,colour=weight)) + facet_grid(year_~origin) 

f:id:axjack:20190228231005p:plain

データミックス(DataMix)のデータサイエンティスト育成コースページにあるスクール紹介動画[0:29]~[0:31]辺りに表示される書籍を出来る限り書き起こす

はじめに

インスタグラムで流れる企業宣伝動画を見ていたらデータミックス(DataMix)のデータサイエンティスト育成コースなるものが目に飛び込んだ。未経験から6か月間でデータサイエンティストを目指すと書いてあり、すごいなぁーと思いながらサイトを眺めていたところ、サイトの真ん中あたりでスクール紹介動画を発見。なるほどーと思いながら動画を見ると、統計関連の書籍が2,3秒ほどちらっと差し込まれていることに気づく。

ということで、これはメモる価値がありそうだなと思ったため、書籍のタイトルを書き起こす次第である。

データミックス(DataMix)とは?

以下のサイトを参照。 datamix.co.jp

書き起こし

当該箇所

以下である。 f:id:axjack:20190217042555p:plain 引用元:https://datamix.co.jp/data-scientist/

書籍リスト

左から順に

感想

割とガチ目の本が載っていてレベル高けぇーと思った。ハイレベル感。

補足

データサイエンティストとはそもそも何か、については下記も併せて読むとイメージが膨らむ。

第四次産業革命スキル習得講座

第四次産業革命スキル習得講座認定制度(METI/経済産業省)に指定されているということをサイトでは謳っていた。どうやらこのキーワードで調べればデータサイエンス講座を網羅的に発掘できそうだ。ということで、経済産業省のサイトを調べると、現在は認定が第三回まで実施されていて次回の認定は、平成31年10月以降に開講する講座を対象として、3月頃に申請受付をする予定とのこと。なお以下はこれまで認定を受けた講座一覧、が掲載されているページである。

不偏分散が不偏であることを確認する

元ネタ

統計の基礎母集団からランダムに取り出した標本でなくても,例えば10個の固定した値を4個と6個に分ける場合,両者を比較するには何で割った分散を使うかという問題についても,一貫して n − 1で割るほうが比較としては正しくなります: を別データに差し替えて確認する。

何を確認するのか

大きさN=9 の xdataを母集団とみなし、大きさn = {2, ... ,9} の標本を抽出する時の、それぞれの不偏分散を確認する。

コード

#データを用意し
xdata <- c(11,22,33,44,55,66,77,88,99)

#不偏分散を自力で計算
huhenbunsan <- (1/( length(xdata) - 1)) * sum( ( xdata - mean(xdata) )^2)

#あらかじめ用意されている不偏分散の関数
var(xdata)
#> var(xdata)
#[1] 907.5

#答え合わせ
var(xdata) - huhenbunsan
#> var(xdata) - huhenbunsan
#[1] 0

#xdataから任意の2つの組み合わせで不偏分散を取り
#それらの平均値を求める
mean( combn(xdata,2,var) )
#> mean( combn(xdata,2,var) )
#[1] 907.5

#xdataから任意の3つの略
mean( combn(xdata,3,var) )
#> mean( combn(xdata,3,var) )
#[1] 907.5

#xdataから任意の4つの略
mean( combn(xdata,4,var) )
#> mean( combn(xdata,4,var) )
#[1] 907.5

#xdataから任意の5つの略
mean( combn(xdata,5,var) )
#> mean( combn(xdata,5,var) )
#[1] 907.5

#xdataから任意の6つの略
mean( combn(xdata,6,var) )
#> mean( combn(xdata,6,var) )
#[1] 907.5

#xdataから任意の7つの略
mean( combn(xdata,7,var) )
#> mean( combn(xdata,7,var) )
#[1] 907.5

#xdataから任意の8つの略
mean( combn(xdata,8,var) )
#> mean( combn(xdata,8,var) )
#[1] 907.5

#xdataから任意の9つの略
mean( combn(xdata,9,var) )
#> mean( combn(xdata,9,var) )
#[1] 907.5

おわりに

感想

標本の不偏分散の平均値が、母集団の不偏分散と全部一致していて気持ちいい(?)です。

combnについて

combn {utils}    R Documentation
Generate All Combinations of n Elements, Taken m at a Time

Description

Generate all combinations of the elements of x taken m at a time. 
If x is a positive integer, returns all combinations of the elements of 
seq(x) taken m at a time. If argument FUN is not NULL, 
applies a function given by the argument to each point. If simplify is FALSE, 
returns a list; otherwise returns an array, typically a matrix. ... 
are passed unchanged to the FUN function, if specified.

所沢市の平成29年版統計書-2 人口(その1)-1.人口の推移から、人口の総数=男+女となっていることを確認する

はじめに

動機

昨今不正統計が話題である。ニュースを見て、我が市の統計データは大丈夫なのかと不安になったので、基礎中の基礎である人口の推移データに不備がないかを調べてみた。

やること

難しい分析をする気は毛頭無之、各年に於ける 人口の総数 = 男性の人口 + 女性の人口をRで確認するだけである。

ところで、市が提供してるのは、

  • htmlファイル
  • Excelファイル

の2形式であるが、私の手元にはExcelがないので、

  • htmlファイルから該当箇所をコピー
  • sublime で加工
  • Rで検算

を試みる。

なお、htmlファイルが公開されているのになぜ{rvest}を使わないのか、という疑問に対する回答として、htmlファイルが構造化されていないことおよびclassidが振られていないことを挙げる。htmlファイルのソースを見ればわかる通り、XPathを頑張って書くのが面倒臭い系のhtmlである。

データソース

www.city.tokorozawa.saitama.jp

コード

sublimeでの加工

1.人口の推移明治9年から平成29年行の一世帯当たり人員2.2までをコピる。続いて、以下のように検索・置換をしてCSVファイルを作成する。

1. 総数-男-女の列を取得する

検索:(.*?)   (.*?)   (.*?)   (.*?)   (.*?)   (.*?)   (.*?)   (.*?)   (.*)
置換:"\4","\5","\6"

2. 数値表現をカンマ有りからカンマ無しへと置換する

検索:(\d),
置換:\1

3. ダブルクオーテーションを削除する

検索:"
置換:

これで総数-男-女の列を数値表現カンマ区切りで加工できたので、CSVとして保存する。

Rで読み込み

先ほどのCSVtoko_jinkou.csvとしてdatフォルダに格納し、列のヘッダを"total","man","woman"と指定して読み込む。

toko_ji <- read.csv("dat/toko_jinkou.csv",col.names = c("total","man","woman"))

結果確認

以下の通り、各年に於ける人口の総数男性の人口 + 女性の人口に等しいことが確認できた。

> toko_ji$total - ( toko_ji$man + toko_ji$woman )
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[52] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

おわりに

結果確認の箇所で0じゃない要素が出てきたらどうしようと少し焦ったものの、まずは一安心。データの前処理方法は沢山あるのでいつでもさらっとスパッと書けるようにしたい。

次の元号をRで予想してみる〜前処理編〜

はじめに

平成もあと3ヶ月足らずで終わってしまいますね。色んな方が次の元号を予想している流れがあるので、自分でも予想してみようと思います。今回は前処理編と称して、ひとまず大化から平成までの元号を取得します。取得するとは、webスクレイピングしてみるということです。どうやってwebスクレイピングするかというと、Rの{rvest}を使ってみます。

コード

rvestの準備

installしてlibraryします。

install.packages("rvest")
library(rvest)
library(dplyr) #テスト用で一瞬使います

{rvest}が使えるかテストします。

> read_html("http://example.com/") %>% html_nodes("p") %>% html_text
[1] "This domain is established to be used for illustrative examples in documents. You may use this\n    domain in examples without prior coordination or asking for permission."
[2] "More information..."   

大丈夫そうなので次に進みます。

rvestでwebスクレイピング→テキスト処理(1文字ずつに分解)

元号国立公文書館 デジタルアーカイブから取得します。が、このサイトは和暦元号・中国元号・朝鮮元号それぞれ記載されているので、後ほど力技で和暦のみ抽出します。 www.digital.archives.go.jp

以下、デジタルアーカイブのサイトから和暦元号を取得 → テキスト処理(1文字ずつに分解)するまで、を一気に書きます。

# table td周りをごそっと取得
src_url <- 'https://www.digital.archives.go.jp/DAS/meta/era#1'
gengou_html <- read_html(src_url)
gengou_html_nd <- html_nodes(gengou_html,"table td")

# 元号のみ取得し
gengou_ <- html_text(gengou_html_nd)[1:5 %% 5 == 1]

# 和暦で最も古い「大化」の添字を見つける
which(gengou_ == "大化")
#[1] 247


# 和暦だけ抽出
gengou_j <- gengou_[1:247]

  
# 確認
head(gengou_j)
#[1] "平成" "昭和" "大正" "明治" "慶応" "元治"

tail(gengou_j)
#[1] "和銅" "慶雲" "大宝" "朱鳥" "白雉" "大化"

# 1文字ずつにパースするが、
gengou_jParse <- strsplit(gengou_j,"")

# パースするとリスト形式となるので、
head(gengou_jParse)[1]
#[[1]]
#[1] "平" "成"

# リストをベクトル形式に変換する
gengou_jParse2 <- unlist(gengou_jParse)

結果確認

いい感じに元号を取得できました。

> gengou_jParse2
  [1] "平" "成" "昭" "和" "大" "正" "明" "治" "慶" "応" "元" "治" "文" "久" "万"
 [16] "延" "安" "政" "嘉" "永" "弘" "化" "天" "保" "文" "政" "文" "化" "享" "和"
 [31] "寛" "政" "天" "明" "安" "永" "明" "和" "宝" "暦" "寛" "延" "延" "享" "寛"
 [46] "保" "元" "文" "享" "保" "正" "徳" "宝" "永" "元" "禄" "貞" "享" "天" "和"
 [61] "延" "宝" "寛" "文" "万" "治" "明" "暦" "承" "応" "慶" "安" "正" "保" "寛"
 [76] "永" "元" "和" "慶" "長" "文" "禄" "天" "正" "元" "亀" "永" "禄" "弘" "治"
 [91] "天" "文" "享" "禄" "大" "永" "永" "正" "文" "亀" "明" "応" "延" "徳" "長"
[106] "享" "文" "明" "応" "仁" "文" "正" "寛" "正" "長" "禄" "康" "正" "享" "徳"
[121] "宝" "徳" "文" "安" "嘉" "吉" "永" "享" "正" "長" "応" "永" "明" "徳" "康"
[136] "応" "嘉" "慶" "至" "徳" "永" "徳" "康" "暦" "永" "和" "応" "安" "貞" "治"
[151] "康" "安" "延" "文" "文" "和" "観" "応" "貞" "和" "康" "永" "暦" "応" "元"
[166] "中" "弘" "和" "天" "授" "文" "中" "建" "徳" "正" "平" "興" "国" "延" "元"
[181] "建" "武" "正" "慶" "元" "弘" "元" "徳" "嘉" "暦" "正" "中" "元" "亨" "元"
[196] "応" "文" "保" "正" "和" "応" "長" "延" "慶" "徳" "治" "嘉" "元" "乾" "元"
[211] "正" "安" "永" "仁" "正" "応" "弘" "安" "建" "治" "文" "永" "弘" "長" "文"
[226] "応" "正" "元" "正" "嘉" "康" "元" "建" "長" "宝" "治" "寛" "元" "仁" "治"
[241] "延" "応" "暦" "仁" "嘉" "禎" "文" "暦" "天" "福" "貞" "永" "寛" "喜" "安"
[256] "定" "嘉" "禄" "元" "仁" "貞" "応" "承" "久" "建" "保" "建" "暦" "承" "元"
[271] "建" "永" "元" "久" "建" "仁" "正" "治" "建" "久" "文" "治" "元" "暦" "寿"
[286] "永" "養" "和" "治" "承" "安" "元" "承" "安" "嘉" "応" "仁" "安" "永" "万"
[301] "長" "寛" "応" "保" "永" "暦" "平" "治" "保" "元" "久" "寿" "仁" "平" "久"
[316] "安" "天" "養" "康" "治" "永" "治" "保" "延" "長" "承" "天" "承" "大" "治"
[331] "天" "治" "保" "安" "元" "永" "永" "久" "天" "永" "天" "仁" "嘉" "承" "長"
[346] "治" "康" "和" "承" "徳" "永" "長" "嘉" "保" "寛" "治" "応" "徳" "永" "保"
[361] "承" "暦" "承" "保" "延" "久" "治" "暦" "康" "平" "天" "喜" "永" "承" "寛"
[376] "徳" "長" "久" "長" "暦" "長" "元" "万" "寿" "治" "安" "寛" "仁" "長" "和"
[391] "寛" "弘" "長" "保" "長" "徳" "正" "暦" "永" "祚" "永" "延" "寛" "和" "永"
[406] "観" "天" "元" "貞" "元" "天" "延" "天" "禄" "安" "和" "康" "保" "応" "和"
[421] "天" "徳" "天" "暦" "天" "慶" "承" "平" "延" "長" "延" "喜" "昌" "泰" "寛"
[436] "平" "仁" "和" "元" "慶" "貞" "観" "天" "安" "斉" "衡" "仁" "寿" "嘉" "祥"
[451] "承" "和" "天" "長" "弘" "仁" "大" "同" "延" "暦" "天" "応" "宝" "亀" "神"
[466] "護" "景" "雲" "天" "平" "神" "護" "天" "平" "宝" "字" "天" "平" "勝" "宝"
[481] "天" "平" "感" "宝" "天" "平" "神" "亀" "養" "老" "霊" "亀" "和" "銅" "慶"
[496] "雲" "大" "宝" "朱" "鳥" "白" "雉" "大" "化"

終わりに

次が果たしてあるのか無いのかいつになるのか、は分かりませんが・・・この情報を元に新元号が発表される4月1日より前になんらかの予想ができたらなぁと。

参考にしたサイト

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.