2020-03-21

対角化

線形代数

示したいこと

正方行列Aに対してR‪⁻¹‬AR=Λと対角化できるような正則行列Rが存在すると仮定する。この時、ΛはAの固有値を並べた対角行列であることを示そう。即ちΛの対角成分λi はAu=auを満たす固有値aであることを示せばよい。ここでu≠0とする。

証明

Au=auよりAu-au=0⇔(A-aI)u=0. 
さてA-aIが逆行列を持つと仮定すると、
(A-aI)‪⁻¹‬(A-aI)u=0⇔u=0となり仮定に反する。

ゆえにA-aIは特異行列であるから|A-aI|=0を満たす。

さて、|R‪⁻¹‬|と|R|を両側から掛け算すると、
|R‪⁻¹‬||A-aI||R|=0より
|R‪⁻¹‬(A-aI)R|=0から
|R‪⁻¹‬AR-aR‪⁻¹‬R|=0⇔|Λ-aI|=0となる。

Λは対角行列なので行列式を展開すると、
(λ1-a)(λ2-a)...(λn-a)=0より
a = λ1,λ2,...,λnとなるから
Aの固有値aは対角行列の対角成分λiに等しいことが示せた
□

参考

対角化された行列の対角成分は固有値 - 理数アラカルト -

2020-03-16

集約関数

SQL

ちょっとgroup byと集約関数と仲良くする必要があったので練習してみた記録です。SQL FiddleのMS SQL Server 2017を使います。

その１

DDL

create table t (
  col1 int 
  ,col2 int
  ,col3 int
  ,col4 int
 )

insert into t values 
(1,1,1,0)
,(1,1,1,0)
,(1,1,1,0)
,(1,1,2,1)
,(1,1,2,0)
,(1,1,2,0)
,(1,1,3,0)
,(1,1,3,0)
,(1,1,3,0)
,(4,5,6,1)
,(4,5,6,-1)
,(5,5,6,1)
,(5,5,6,1)
,(6,5,6,0)
,(6,5,6,null)

DML

-- --------------------------------------------------------------------------------
-- http://mickindex.sakura.ne.jp/database/celko/celko_tis.html
-- col1、col2、col3 の三列でグループ化したときに
-- col4 がすべて 0 であるような行を一意に取得せよ。
-- --------------------------------------------------------------------------------
select 
  col1,col2,col3
from 
  t
group by 
  col1,col2,col3
having 
  max(col4) = min(col4) -- group 内全てのcol4が同一な値である
  and min(col4) = 0 -- 最小値が0である
  and count(*) = count(col4)  -- col4はnullを含んでいない
;

補足

最小値の抑え込み

  and min(col4) = 0 -- 最小値が0である

がないと、

,(5,5,6,1)
,(5,5,6,1)
,(6,5,6,0)
,(6,5,6,null)

が抽出されてしまう。

nullと集約関数

は中身がnullだとエラーになったり省略されたりする*1。なので、ISNULLやcoalesceでnullを抑え込むのもあり。

その２

count(*)とcount(X)とcount(distinct(X))の違いをみてみます。なお、count(distinct(X))は奥が深い問題*2のようです。

DDL

create table t(
  grp char(1) null 
   , val nvarchar(10) null
)

insert into t values 
('A',null)
,('B','bbb')
,('B','bbb')
,('B','ccc')
,('C','ddd')
,('C','eee')
,('D','fff')
,('D',null)
,('D',null)

DML

select
  grp
  ,max(val) as max_val
  ,min(val) as min_val
  ,count(*) as count_aster
  ,count(val) as count_val
  ,count(distinct(val)) as count_distinct_val
from t
  group by grp

実行結果

grp	max_val	min_val	count_aster	count_val	count_distinct_val
A	(null)	(null)	1	0	0
B	ccc	bbb	3	3	2
C	eee	ddd	2	2	2
D	fff	fff	3	1	1

*1:「警告 : NULL 値は集計またはその他の SET 演算で削除されました」が出たりする

*2:https://en.wikipedia.org/wiki/Count-distinct_problem

2020-02-24

コーシー分布

復習確率勉強

準備：微分の公式など

$\displaystyle (\frac{f}{g})' = \frac{f'g - fg'}{g^2} \\ \displaystyle (\sin \theta)' = \cos \theta \\ \displaystyle (\cos \theta)' = -\sin \theta \\ \displaystyle (\tan \theta)' = (\frac{\sin \theta}{\cos \theta} )' = \frac{(\sin \theta)'\cos \theta - \sin \theta(\cos \theta)'}{\cos^{2} \theta} = \frac{ 1 }{ \cos^{2} \theta } \\ \displaystyle ( 1 + \tan^{2} \theta)^{-1} = ( \frac{1}{\cos^2 \theta} )^{-1} = \cos^{2} \theta$

コーシー分布の積分

$\displaystyle f(x) = \frac{1}{ \pi ( 1 + x^2) }$ を $x \in \mathbb{R}$ で積分すると、 $\displaystyle x = \tan \theta$ と置換し、

$\displaystyle \int^{\infty}_{-\infty} \frac{ 1 }{ \pi( 1 + x^2 ) } dx = \int^{\frac{\pi}{2}}_{-\frac{\pi}{2}} \frac{1}{\pi} \cos^{2} \theta \frac{1}{\cos^2 \theta } d\theta = \int^{\frac{\pi}{2}}_{-\frac{\pi}{2}} \frac{1}{\pi} d\theta = \frac{1}{\pi} \left[ \theta \right]^{\frac{\pi}{2}}_{-\frac{\pi}{2}} = \frac{1}{\pi} ( \frac{\pi }{ 2} - \frac{-\pi}{2}) = 1$

となる。

2019-12-29

2019年大掃除大会まとめ

揃えるもの

液体固体

カビキラー
バスマジックリン
油汚れマジックリン ← 心強い
パイプユニッシュ
いつもの台所洗剤

ゴシゴシ系

歯ブラシ
メラミンスポンジ ← ボロボロになりがち
ぞうきん ← 便利。乾拭き重要
スポンジ
アルコール除菌ペーパー

防具

マスク
ゴム手袋
新聞紙的な何か

心得ておくこと

ワンクール1時間を限度に適宜休むこと
長時間のしゃがみは腰に来る
換気はしっかり
油汚れは油汚れマジックリンが最強
重曹は割とうまくいかない。上級者向けの技である為、安易に使用しないこと。
まずは表面を軽く落とす
てか大掃除は冬より夏の方が良いと思う( ˘ω˘ )

2019-12-23

『確率４万分の１、県の入札くじに６回連続当選「奇跡的」』なのかどうかを調べよう

R 調べ物統計

www.asahi.com

動機

面白いニュースを見かけたので計算してみましたなポストです。記事のこの部分に注目しました。

６回のくじにはそれぞれ制限価格で入札した３～８の鑑定業者が参加した。６回連続で当たる確率を計算すると約４万分の１となる。

奇跡ってどのくらいの割合で起きるのかな？もうすぐクリスマスだしってことでRでやってみます。

まず引用箇所を表にすると、

i回目	鑑定業者数 X_i
1	X_1
2	X_2
3	X_3
4	X_4
5	X_5
6	X_6

ここでX_i ∈ [3,8]の整数で、「６回連続で当たる確率を計算すると約４万分の１」を満たすX_iであって欲しいってことでX_iを全部掛け算したら4万になれば良いでしょう。

とりあえずそのような鑑定業者数の組み合わせを力技で探すと、

kgyousya <- c(3,4,5,6,7,8)

kaisuu <- data.frame(
  x1=kgyousya
  ,x2=kgyousya
  ,x3=kgyousya
  ,x4=kgyousya
  ,x5=kgyousya
  ,x6=kgyousya
)

kaisuu.d <- expand.grid(kaisuu)
kaisuu.m <- as.matrix(kaisuu.d)

kaisuu.prod <- apply(kaisuu.m,MARGIN = 1,prod)

table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] )

#> table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] )
#
#39200 40000 40320 40960 
#180    15   720    60

915ケースありました。

次に、掛け算して39,200〜40,960となるような鑑定業者数の組み合わせ：915ケースで、各ケースを10000回シミュレーションする。そして各ケースで確率4万分の1以下が10000回中何回起きるかを調べる。(←奇跡の割合）最後に奇跡の割合の分布を観察しましょう。

# 候補のケース数: 915
#> nrow(kaisuu.sub)
#[1] 915
mycase <- nrow(kaisuu.sub)
results <- numeric(mycase)

# シミュレーション回数
ntimes <- 1:10000
result <- numeric( length(ntimes) )

for(j in 1:mycase){
  for(i in ntimes){
    np <- sample(c(0,1),size=6,replace = TRUE, prob=c(.5,.5))
    result[i] <- prod( abs( np  -  (1/kaisuu.sub[j,]) ) )
  }
  results[j] <- sum( result < 1/40000 )/length(result)
}

結果

度数分布表

> table(results)
results
    0 0.005 0.006 0.007 0.008 0.009  0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 
  180     1     3     7     8    17    21    56    63    66    64    54    74    74    49 
0.019  0.02 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 
   46    45    30    18    19     7     4     4     3     1     1

要約値

> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.01000 0.01400 0.01223 0.01700 0.02900

どんなケースでも平均すれば1%ぐらいは起きるのでは？という見方で良いのかなぁ・・・。

プロット

f:id:axjack:20191223212612p:plain

f:id:axjack:20191223212627p:plain

感想

こういうの式だけ立てて計算できるようになりたいですなぁ。

2019-11-19

2018年11月統計検定2級問9

勉強復習

この問題をなぜ解くのか
問題文の概要
解答
- 問(1)
- 問(2)

この問題をなぜ解くのか

前にこんなこと

統計検定であって確率検定ではないので、サイコロを７回投げる時２以下の目が出る確率を解けるようになってもなんだかなぁ
統計検定２級に合格したので勉強法やら参考書などを書いておきます。 - axjack's blog

を書きました。がしかし、統計検定準1級の勉強をしていく中「やっぱ確率もちゃんと解けないとダメじゃん😢」と気持ちが一転したので改めてちゃんと解いてみようと思った次第です。

問題文の概要

どの面も等確率で出る6面のサイコロを7回投げる。2以下の目が出る回数をXとする。

問(1)：P(X = x+1)とP(X = x)の比、i.e. $\frac{ P(X=x+1) }{ P(X=x) }$ を求めよ
問(2)：P(X = x) が最大となるxは?

解答

問(1)

$\displaystyle X_i$ を出目が2以下の時に1、そうでない時を0とすると $\displaystyle X_i$ はベルヌーイ分布：Bin(1,p)に従う。ここでp=1/3である。この時、 $\displaystyle X = \sum X_i$ は二項分布：Bin(n,p)に従う。よって、確率質量関数は $\displaystyle P(X=k)= {}_n \mathrm{ C }_k p^{k} (1-p)^{n-k}$ となる。ここで、n=7, p=1/3を代入するとP(X=x)とP(X=x+1)はそれぞれ、

$\displaystyle P(X=x) = {}_7 \mathrm{ C }_x (1/3)^x (2/3)^{7-x}$
$\displaystyle P(X=x+1) = {}_7 \mathrm{ C }_{x+1} (1/3)^{x+1} (2/3)^{7-(x+1)}$

となる。すると、

$\displaystyle \frac{P(X=x+1)}{P(X=x)} \\ \displaystyle = \frac{ {}_7 \mathrm{ C }_{x+1} (1/3)^{x+1} (2/3)^{7-(x+1)} } { {}_7 \mathrm{ C }_x (1/3)^x (2/3)^{7-x} } \\ \displaystyle = \frac{ {}_7 \mathrm{ C }_{x+1} (1/3)^{x} (1/3) (2/3)^{7-x} (2/3)^{-1} } { {}_7 \mathrm{ C }_x (1/3)^x (2/3)^{7-x} } \\ \displaystyle = \frac{ \frac{ 7! }{ (7-(x+1))! (x+1)! } (1/3)^{x} (1/3) (2/3)^{7-x} (2/3)^{-1} }{ \frac{7!}{(7-x)!x!} (1/3)^x (2/3)^{7-x} }\\ \displaystyle = \frac{ \frac{ 7! }{ (6-x)! (x+1)x! } (1/3)^{x} (1/3) (2/3)^{7-x} (2/3)^{-1} }{ \frac{7!}{(7-x)(6-x)!x!} (1/3)^x (2/3)^{7-x} } \\ \displaystyle = \frac { \frac{ 1 }{ x+1 } (1/3) (3/2) }{ \frac{1}{7-x} } = \frac{ 7-x }{ 2x+2 }$

となる。よって、 $\displaystyle \frac{P(X=x+1)}{P(X=x)} = \frac{ 7-x }{ 2x+2 }$

問(2)

P(X=x)が最大となるxを $\displaystyle \frac{P(X=x+1)}{P(X=x)} = \frac{ 7-x }{ 2x+2 } = Q(x)$ を用いて考えてみると、Q(x)が1より大であればP(X=x+1) > P(x)であると言えるし1未満であればP(X=x+1) < P(x)と言える。具体的にxへ0から6を代入すると、

x	Q(x)	判定
0	7/2 > 1	P(1) > P(0)
1	6/4 > 1	P(2) > P(1)
2	5/6 < 1	P(3) < P(2)
3	4/8 < 1	P(4) < P(3)
4	3/10 < 1	P(5) < P(4)
5	2/12 < 1	P(6) < P(5)
6	1/14 < 1	P(7) < P(6)

となる。関数の増減表のように判定の列を眺めれば、P(x)はP(2)がピークであることからP(X = x) が最大となるxは2である。

2019-11-04

超幾何分布

統計勉強復習 R

何度やっても忘れるのでブログに書いて覚えよう。

超幾何分布の確率質量関数

$\displaystyle P(X=x) = f(x; N,M,n) = \frac{ \binom{M}{x} \binom{N-M}{n-x} } { \binom{N}{n} }$

計算例

その1

男50人女50人から10人を選ぶ。10人のうち男3人女7人となる確率pは、超幾何分布を用いて

$\displaystyle p = \frac{ \binom{50}{3} \binom{50}{7} }{ \binom{100}{10}}$

その2

統計検定準一級2017年問10より引用。表の80人から30人を無作為抽出する。その30人のうち男性かつ就職している人数Xは超幾何分布に従う。

＼	就職	非就職	計
男性	38	3	41
女性	30	9	39
計	68	12	80

80人から30人を選ぶ
「男性かつ就職」である人：38
「男性かつ就職」でない人：80-38 = 42

と整理して、
$\displaystyle P(X=x) = \frac{ \binom{38}{x} \binom{42}{30-x} }{ \binom{80}{30} }$

となる。

その3

同じ表を用いて、80人から20人を無作為抽出する。その20人のうち「女性かつ非就職」な人数Xは超幾何分布に従う。

80人から20人を選ぶ
「女性かつ非就職」である人：9
「女性かつ非就職」でない人：80-9 = 71

と整理して、
$\displaystyle P(X=x) = \frac{ \binom{9}{x} \binom{71}{20-x} }{ \binom{80}{20} }$

となる。

Rで計算する

choose関数で立てた式が超幾何分布の質量関数dhyperと同じことを確認。

# その1
> choose(50,3)*choose(50,7)/choose(100,10)
[1] 0.1130964
> dhyper(3,50,50,10)
[1] 0.1130964

dhpyer(x,m,n,k)の引数*1は、

x: vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
m: the number of white balls in the urn.
n: the number of black balls in the urn.
k:the number of balls drawn from the urn.

*1:https://www.rdocumentation.org/packages/stats/versions/3.6.1/topics/Hypergeometric

2019-09-26

連の検定(Runs Test for Detecting Non-randomness)

統計調べ物 R

はじめに

連の検定についてよく分からなかったのでRで実装して確かめてみました。

連とは？

２値{A,B}を取る系列ABAABBAAABBBABがあった時、これをA | B | AA | BB | AAA | BBB | A | Bと連続する同じ値ごとに分割できるように見える。この時同じ文字または記号のひと続きを連(run)*1と言い、またひと続きの個数を連の長さ*2と言う。この場合、Aの連が4つ(= {A,AA,AAA,A})とBの連が4つ(= {B,BB,BBB,B})あることが分かる。AAAAAAAAABBBBBBBBBBBであれば、Aの連が1つでBの連が1つとなる。

確かめたこと

Aの数(= Na)とBの数(= Nb)が与えられた時、N = Na+Nbの大きさを持つ系列は、

もしA・Bがランダムに現れるなら(←帰無仮説)連の数はどのような分布になるのか　→　正規分布
連の数の平均値・分散が以下の式*3どおりとなるか

f:id:axjack:20190926230103p:plain

を確かめた。

ソース

# 連の検定####
rm(list=ls())

# --------------------
# 関数: generateRun ####
# --------------------
# 引数
#   Na_: Aの個数
#   Nb_: Bの個数
# 戻り値
#   連の数を返す関数
# --------------------
generateRun <- function(Na_, Nb_){
  return(
      function(dummy){
        N <- Na_ + Nb_
        Arep <- rep(TRUE,Na_)
        Brep <- rep(FALSE,Nb_)
        AB <- c(Arep,Brep)
        xx <- sample(AB, N, replace = FALSE)
        yy <- xx[1:(N-1)]
        zz <- xx[2:N]
        return( sum( yy != zz)  + 1 )
    }
  )
}

# シミュレーション ####
## 各種変数設定
Na <- 9 # Aの個数
Nb <- 11 # Bの個数
sim_N <- 10000 # シミュレーション回数

## 平均値と分散
( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値
( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散

# シミュレーション実行
Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) )

# 要約
table(Result.Run)
summary( Result.Run )
var( Result.Run )

# グラフ表示
hist(Result.Run,breaks = seq(min(Result.Run),max(Result.Run),1))
par(new=T)

# 正規分布を重ね合わせ
curve(dnorm(x,t_mean,sqrt(t_var)),min(Result.Run),max(Result.Run),xlab="",ylab="",xaxt="n", yaxt="n",col="red",bty="n")

結果の確認

数値の確認

> # シミュレーション ####
> ## 各種変数設定
> Na <- 9 # Aの個数
> Nb <- 11 # Bの個数
> sim_N <- 10000 # シミュレーション回数
> ## 平均値と分散
> ( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値
[1] 10.9
> ( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散
[1] 4.637368
> # シミュレーション実行
> Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) )
> # 要約
> table(Result.Run)
Result.Run
   2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
   1    3   13   28  167  312  852 1218 1725 1699 1699 1144  704  279  121   30    5 
> summary( Result.Run )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.00    9.00   11.00   10.89   12.00   18.00 
> var( Result.Run )
[1] 4.666794

プロット

f:id:axjack:20190926230756p:plain

結果の考察

計算した平均値・分散がシミュレーションの結果とだいたい合っていることがわかった。
ヒストグラムとカーブプロットが大まかにあっている。

おわりに

帰無仮説の下での平均値・分散の式が奇妙で覚えられないし式を見てもあまり実感が湧かなかったものの、シミュレーションしてみると確かに平均値・分散がそれっぽく合うのでこの式の妥当性を確認することができた。疑わしい系列があってもこれで確かめることができそうだ。

参考

連の検定関連

R言語のplot関連

49. 図の重ね合わせ
- 図の重ね合わせ方法を参考にしました。
正規分布のグラフをRで描く。【curve()の使い方】
- 正規分布をcurveで書く方法を参考にしました。
R 使い方軸・ラベルの調整（向き・サイズ・色など）グラフの描き方
- plotを重ねる時に軸やラベルを消す方法を参考にしました。

*1:[連の検定]より引用

*2:[連の検定]より引用

*3:式は[連の検定]より引用