2020-04-29

行簡約行列をRで

R 線形代数

pracmaのrrefを使って行簡約行列を出してみます。これで適当な行列を手計算で簡約化して答え合わせができますね。

ソース

# install.packages('pracma')
# library('pracma')

## 行ベクトルを4本
a1 <- c(1,2,0)
a2 <- c(2,4,0)
a3 <- c(0,1,3)
a4 <- c(1,3,3)

## 行列にする
M <- matrix(c(a1,a2,a3,a4),nrow = 3, ncol=4)

## 行列を表示
print(M)

## 行簡約行列を表示 
rref(M)

## 行列の階数を表示
Rank(M)

実行結果

> a1 <- c(1,2,0)
> a2 <- c(2,4,0)
> a3 <- c(0,1,3)
> a4 <- c(1,3,3)
> M <- matrix(c(a1,a2,a3,a4),nrow = 3, ncol=4)
> print(M)
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    1
[2,]    2    4    1    3
[3,]    0    0    3    3
> Rank(M)
[1] 2
> rref(M)
     [,1] [,2] [,3] [,4]
[1,]    1    2    0    1
[2,]    0    0    1    1
[3,]    0    0    0    0
>

2020-04-17

update文

SQL

とある１レコードの値で本番のデータを全件更新してしまうというヤラカシをやってしまったので、次回は気をつける意を込めた記録です。 SQLはcompile sql server onlineを使用します。

仕様

tbl1.val1をtbl2.val2へと更新する
更新の条件
- tbl1.idとtbl2.idを結合し、tbl2.val2 = 'foo'の時だけ更新する
- tbl2.val2 <> 'foo'の時は更新しない

tbl1

f:id:axjack:20200417235435p:plain

tbl2

f:id:axjack:20200417235458p:plain

ということでこの場合、id = {10,30,80}が更新対象となります。

問題のコード(ヤラカシあり)

-- データ作成 ここから
create table tbl1  (
 id int primary key
  , val1 nvarchar(10) null
)
;

insert into tbl1 values
  (10,'a')
  ,(20,'b')
  ,(30,'c')
  ,(40,'d')
  ,(50,'e')
  ,(60,'f')
  ,(70,'g')
  ,(80,'h')
  ,(90,'i')
;

create table tbl2 (
 id int primary key
  , val2 nvarchar(10) null
  , constraint fk_id foreign key(id) references  tbl1
)
;

insert into tbl2 values
  (10,'foo')
  ,(30,'foo')
  ,(50,'baz')
  ,(60,'baz')
  ,(80,'foo')
;
-- データ作成 ここまで



--  一時テーブルに tbl1 のレコードを格納
select * into #t1 from tbl1
;


-- tbl1 更新前
select * from tbl1
;

-- 更新処理
update tbl1
set val1 = tbl2.val2
from #t1 
  inner join tbl2 
    on #t1.id = tbl2.id 
    and tbl2.val2 = 'foo'
;


-- tbl1 更新後
select * from tbl1
;

実行結果

上：更新前のtbl1, 下：更新後のtbl1なのですが、見ての通り残念な感じに更新されてしまいました。。

f:id:axjack:20200417234647p:plain

better なコード

★の箇所のように、updateで更新対象となるテーブルをinner join してあげるとinner join で抽出される結果通りに更新されます。

update tbl1
set val1 = tbl2.val2
from #t1 
  inner join tbl2 
    on #t1.id = tbl2.id 
    and tbl2.val2 = 'foo'
  inner join tbl1 on tbl1.id = #t1.id -- ★
;

better なコードの実行結果

仕様を満たした更新結果となりました。 f:id:axjack:20200417235104p:plain

まなび

更新対象のテーブルをちゃんとjoin しましょう。

2020-04-05

カーネルは線型空間

線形代数

示したいこと

線型空間V, W に対して線型写像 $\displaystyle T:V \rightarrow W$ が与えられているとする。
この時、カーネル $\displaystyle \rm Ker \ T = \{ x \in V | T(x) = 0_w \in W \}$ は線型空間であることを示す。

証明

$\displaystyle \forall x, y \in \rm Ker \ T$ および $\displaystyle \forall a, b \in \mathbb{R}$ をそれぞれ取る。
すると、 $\displaystyle ax+by \in V$ より、 $\displaystyle T(ax+by) = T(ax) + T(by) = aT(x) + bT(y) = a0_w + b0_w = 0_w + 0_w = 0_w \in W$
となるので、Ker T の任意の元はベクトル和とスカラー倍に関して閉じていることが分かった。よってKer T は線型空間であることが示せた。

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

感想

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