本番で解けなかった問題を解いてみます。今回は問１です。この問題で問われていることは
です。似たような問題ないかなぁと調べてみるとなんと、アクチュアリーの昭和39年数学1の問３にほぼ同じような問題があったので、これを見つつ今回の問題を解いてみることにしました。
以下、
を代入すると記述1〜記述4の解となります。
ポアソン分布は、
です。
の時、
が従う分布を求めます。
(はてなのmathjax記法がエラるので画像にて。。)
となるので、
が和の分布となる(ポアソン分布の再生性が確認できた)。よって、和の分布は
となる。
求めるのは、
という確率分布である。計算すると、
となるので、この条件付き確率分布は二項分布となる。
したがって、平均値は
となる。
テーマはロジスティック回帰。問題文に記載された出典を調べると、Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failureがヒットするので、このデータを用いれば再現が出来る(のだろうか)。ひとまずやってみると、
# Field Erosion or blowby #### yF <- c(0,1,0,0,0,0,0,0,1,1,1,0,0,2,0,0,0,0,0,0,1,0,1) TD <- ifelse(yF > 0, 1, 0) Temperature <- c(66,70,69,68,67,72,73,70,57,64,70,78,67,53,67,75,70,81,76,79,75,76,58) r <- glm(TD ~ Temperature, family = binomial(link = logit))
> summary(r)
Call:
glm(formula = TD ~ Temperature, family = binomial(link = logit))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.0665 -0.7679 -0.3844 0.4541 2.2047
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 14.9246 7.4235 2.010 0.0444 *
Temperature -0.2302 0.1087 -2.117 0.0343 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 28.267 on 22 degrees of freedom
Residual deviance: 20.511 on 21 degrees of freedom
AIC: 24.511
Number of Fisher Scoring iterations: 5
やや違うがおおよそ問題文の出力結果と合致している。
plot(Temperature, TD) points( Temperature,fitted(r), col='red')
y_i は 互いに独立な確率変数Y_i の実現値であり、Y_i は (ア) に従う。 π_i( 0 < π_i < 1 ) について、構造式(イ) を仮定する。
とある。ロジスティック回帰なのでY_iはポアソン分布に従い、構造式(リンク関数？)はロジットなのでlog( (π_i / ( 1 - π_i) ) =である。
0 = 15.0429 - 0.2322 x_i をx_i について解くと、x_i = 64.78423772
確率を求めるのでロジットではなくロジスティック関数に戻して解いてみる。すると、π_i = 1 / ( 1 - exp( - (15.0429 - 0.2322x_i) ) ) で x_i <- 31 としてπ_iを求めれば良い。結局、
> 1/(1+exp(-7.8447)) [1] 0.9996083
となる。
試験㆗、「聞いたことあるけど自分で解いたことないんだよなぁ、でも解けたら応用範囲が広そうだなぁ」と思っていた問題の１つ。
受験します。
ということでまた一年間みっちり勉強しようと思う所存です。とはいえずっと勉強だとアレなので、適宜Rを用いながらenjoyしようかとも思います。
今日で令和元年ゴールデンウィークも無事終了です。
前回に引き続きニュートン法で平方根を求めてみました。今回は、N = [1,1000] の範囲でニュートン法を用いて10回反復近似計算する時、真の値と近似値との残差はどんな感じになるのか、を調べてみます。
全てのNに対して初期値１から始めると、Nが大きくなればなるほど収束が遅れるような気配もありますが、まぁそこは部屋の片隅に置いておきましょう。
aの平方根を求めるには、初期値を設定し、
を順次反復計算します。今回はとおきました。
a = 2, の時、
初期値が1であっても、７回ほど反復計算しておけば問題なさそう。
#rm(list=ls())
# f(x) = x^2 - a をニュートン法で求める更新式
fnewton <- function(xn=1,a=2){
return ( xn - ( xn^2 - a )/( 2 * xn ) )
}
# 初期化
how_many_times_rep <- 10
N <- 1000
xs <- 1:N
root_xs <- sqrt(xs)
zeros <- rep(0,length(xs))
# 結果格納用行列
mat.x <- cbind(
matrix(c( xs,root_xs ),nrow = length(xs))
,matrix( rep(zeros, how_many_times_rep), nrow = length(zeros))
)
mat.x[,3] <- fnewton(xn = 1, a = mat.x[,1])
# ニュートン法にて更新
for(i in 4:(how_many_times_rep+2) ){
mat.x[,i] = fnewton(xn = mat.x[,i-1], a = mat.x[,1])
}
#mat.dif <- matrix(0, nrow = nrow(mat.x), ncol = ncol(mat.x))
mat.dif <- mat.x
# 真の値 - 近似値
for(i in 3:(how_many_times_rep+2)){
mat.dif[,i] <- mat.x[,2] - mat.x[,i]
}
# 結果確認
summary(mat.dif[,-c(1,2)])
par(family="Osaka")
boxplot(mat.dif[,-c(1,2)],
col="lightpink",
names=paste0(1:how_many_times_rep,"回目"),
las=2,
main=paste0("真の値と近似値との残差の箱ひげ図(N = 1 ~ ",N,")")
)
summary の値を書き出すと、
V1 V2 V3 V4 V5
Min. :-468.9 Min. :-219.6 Min. :-95.99 Min. :-36.103 Min. :-9.6226
1st Qu.:-348.2 1st Qu.:-161.4 1st Qu.:-69.00 1st Qu.:-24.697 1st Qu.:-5.8550
Median :-228.4 Median :-104.0 Median :-42.79 Median :-14.052 Median :-2.7105
Mean :-229.7 Mean :-105.3 Mean :-44.04 Mean :-15.158 Mean :-3.3997
3rd Qu.:-110.0 3rd Qu.: -48.1 3rd Qu.:-18.09 3rd Qu.: -4.824 3rd Qu.:-0.5632
Max. : 0.0 Max. : 0.0 Max. : 0.00 Max. : 0.000 Max. : 0.0000
V6 V7 V8 V9
Min. :-1.122493 Min. :-1.924e-02 Min. :-5.849e-06 Min. :-5.400e-13
1st Qu.:-0.515571 1st Qu.:-4.763e-03 1st Qu.:-4.140e-07 1st Qu.:-3.553e-15
Median :-0.146451 Median :-4.762e-04 Median :-5.069e-09 Median : 0.000e+00
Mean :-0.293577 Mean :-3.268e-03 Mean :-5.871e-07 Mean :-3.005e-14
3rd Qu.:-0.009673 3rd Qu.:-2.953e-06 3rd Qu.: 0.000e+00 3rd Qu.: 0.000e+00
Max. : 0.000000 Max. : 0.000e+00 Max. : 0.000e+00 Max. : 3.553e-15
V10
Min. :-3.553e-15
1st Qu.: 0.000e+00
Median : 0.000e+00
Mean : 3.242e-17
3rd Qu.: 0.000e+00
Max. : 3.553e-15
でした。
ですね。7回電卓ぽちぽちすれば平方根が得られそうです。
残差の範囲は反復回数が増す毎にどんどん小さくなっていきますね。５回目あたりで箱が見えなくなり、６回目あたりでヒゲもぺっちゃんこに見えます。
計算方法覚えておけばルート機能要らない気もしないでも無いですね。脳内記憶容量を節約するならやはりルート付き電卓を買いましょう。
電卓のルートってどうやって求めるんだっけ？ → たぶんニュートン法だろうなぁ。ということで、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)
6回目の反復あたりでほぼ答えが電卓と一致して驚いた。
平成もそろそろ終わりそうな今日この頃、ソフトバンクから格安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開通手続き |
構成プロファイルのダウンロードはこちらの案内に従い、http://s.ocn.jp/ltecfg へアクセス、構成ファイルをダウンロードした。この時、端末のSIMカードはソフトバンクのままの状態。
MNP開通手続きのご案内に従い、http://s.ocn.jp/mnp01 へアクセス、申し込み時に発行されたocnのメールアドレスとパスワードを入力し、MNP開通手続きを実行した。この時、端末のSIMカードはソフトバンクのままの状態。
実行後、iPhoneアンテナマークの右にある回線事業者名：Softbankが消えた。(この時点でMNP開通手続きが完了した、と思われる。)
iPhone8はナノSIMなのでSIMカードをナノサイズに切り取り、iPhone8からソフトバンクのSIMカードを取り出して交換した。
すると、回線事業者名がdocomoに！
の２つのアプリをインストールした。
以下を問題なく実施できた。
手続きがものすごく簡単だった。SIMカードが届くまで１週間掛かったものの、スマホは何事もなく使えたしなんら問題なく。
SIMカードの取り出しが難しかった。慣れてないとコツをつかむまでに苦労しますね。まぁ、勢いよくぐいっと押せば蓋が開くのですが・・・。
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}
で引っかかるのがデータレコード、となります。
(.{30})を\1\nと置換すると、30byte引っかかるごとに改行を挟むことになるので、結果的に分割できます。なお改行コードが気になる方は文字コード環境に合わせて\nなり\r\nなり¥r¥nなりに変更してください。
機械学習や人工知能の波が押し寄せようと来まいと、データフォーマットの変換がなくなることはないでしょう。固定長・CSV・XML・JSOなど、どんな何が来てもグリっと変換できるようにしておけば、平成の次の時代でも困ることはないはず(たぶん)。
唐突にテキストファイルだけ渡された時、しれっとフォーマット変換しなければいけない、だが手元にテキストエディタしかないのだが？な時にこの記事が役に立ったら幸いです。
\nをブランク(≠半角スペース)で置換すれば元に戻ります。
ggplotやqplotを日々練習しています。ggplot, geom_point, facet_gridを組み合わせると、
のように５変数くらい同時に、割と分かりやすく表示することができる。分割が煩わしいなら、facet_gridの項を消せば１つの散布図色付きのように表示することができる。
library(dplyr) library(ggplot2) ggplot(mpg) + geom_point(aes(cty,displ,colour=fl)) + facet_grid(year~cyl)
library(ISLR) a1 <- Auto %>% mutate(year_ = ifelse(year <= 80, 70, 80)) ggplot(a1) + geom_point(aes(displacement,horsepower,colour=weight)) + facet_grid(year_~origin)
