平均、分散、標準偏差

2016年8月19日 sikino コメントする

(i) 全データが存在する場合

　例１：４０人クラスにおいてテスト点数の解析を行いたい(４０人全員の点数を知っている)。
　例２：過去１年間、365日で体重がどのくらい増減したかの解析をしたい(365日分のデータが存在する)。
平均$m$
$
\displaystyle m=\frac{1}{N}\sum_{i=1}^N x_i
$
分散$v$
$
\displaystyle v=\frac{1}{N}\sum_{i=1}^N (x_i-m)^2
$
標準偏差$\sigma$
$
\displaystyle \sigma=\sqrt{v}
$

平均値から$\pm 66\%$の誤差に収まっている範囲は、
平均m ± 標準偏差σ　と書ける。
例1 のクラスの点数の場合は、
平均m ± 標準偏差σ [点]　　となる。

(ii)一部分のデータしか存在しない場合

例１：４０人クラスにおいて２３人しかテストの点数を知らないが、その２３人からクラス全体の点数の解析を行いたい（２３人の点数は分かる。残り１７人の点数は不明）。
　例２：過去１年間、365日で体重の増減を解析したいが、200日分のデータしかない。だけど１年の解析を行いたい（165日分のデータは不明、紛失した）。
　　　　例３：物理量の測定（測定回数が1000回だとしても、物理量は無限回の測定をしなければ真値は不明だ、と考える。）

この場合、平均、分散、標準偏差は
平均$m’$
$
\displaystyle m’=\frac{1}{N}\sum_{i=1}^N x_i
$
分散$v’$
$
\displaystyle v’=\frac{1}{N}\sum_{i=1}^N (x_i-m’)^2
$
標準偏差$\sigma’$
$
\displaystyle \sigma’=\sqrt{v’}
$
平均値から$\pm 66\%$の不確かさに収まっている範囲は、
平均m’ ± 標準偏差σ’　と書ける。
例1 のクラスの点数の場合は、
平均m’ ± 標準偏差σ’ [点]　　となる。

誤差と不確かさ　は違うので注意しましょう。
誤差 → 完璧なデータのときに使う。
不確かさ → 本当の値を知らないときに使う。

プログラミングと数値計算

√xの数値積分(端点で特異点が存在する場合)

2016年8月18日 sikino コメントする

$\sqrt{x}$を$[0\sim b]$の範囲で数値的に積分したいとします。
この時に問題となるのは$\sqrt{x}$が$x=0$で特異点(この場合、一階微分が発散する)を持ってしまうことです。
解決するには変数変換を行い、
$
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{b^{1/2}} 2x^2 dx
$
として右辺を計算するのが良いでしょう。

端点特異点を持つ場合、数値積分は工夫しなければなりません。
ここでは$\sqrt{x}$を$0$から$b$まで積分することを考えます。

問題点

台形則、シンプソン則などのニュートン・コーツ型の数値積分は点と点をラグランジュ補間(多項式による補間)を行って積分します。
そのため、多項式で表せない関数を数値積分するときは著しく積分精度が落ちてしまいます。
ガウス・ルジャンドル求積法は次数が高いだけでニュートンコーツ型と同じラグランジュ補間を行うため精度が良くありません。
こういった時には台形則を組み合わせて精度評価が出来るロンバーグ積分法による方法が推奨されますが、結局は台形則なので、多項式の表現から抜け出ることはできず、被積分関数の評価回数がかなり多くなってしまいます。

ガウス求積法の一部にラゲール多項式を用いるガウス・ラゲール求積法があります。
これは
$x=0$で$x^{\alpha},(0\lt \alpha)$の振る舞い
かつ
$x=\infty$で被積分関数が$e^{-x}$で減衰する関数
に対して非常によい方法ですが、これを使う際の積分範囲は$[0\sim \infty]$であり、今回の場合うまく用いることが出来ません。

最善の方法は二重指数関数型数値積分公式(DE公式)という積分法でしょう。
補間するのではなく、変数変換によって積分する空間を変えて積分する方法で、端点特異点がある場合に強い方法です二重指数関数型数値積分(理論とプログラム)。
端点でDE公式も使えないほどの非常に強い特異性を持つ場合、超函数法(hyperfunction method)と呼ばれる数値積分法が良いそうです。

ある程度の精度で簡単に計算するのならば、積分区間を分割します。
被積分関数が特異性を発揮するのは$x=0$付近ですから、積分
$
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{\Delta}\sqrt{x} dx + \int_{\Delta}^b\sqrt{x} dx
$
にし、右辺第1項を低次の方法(例えば台形則)で分点数を多くとり計算します。第2項には特異性を含まないので荒い積分方法で良いでしょう。

実際に$\sqrt{x}$の数値積分が困難であることを数値計算で示します。
$
\displaystyle \int_0^1 \sqrt{x} dx = \frac{2}{3} = 0.66666666666\cdots
$
を各々の方法で計算しますと、

0.666172080968984 短冊近似　　　(1000分割)
0.666660362218984 台形則　　　　(1000分割)
0.666664189108662 シンプソン則　(1000分割)
0.666664574081247 ロンバーグ則　(1000分割,相対誤差12桁一致)
0.666666768098704 ガウス・ルジャンドル求積(100点)

という結果を得ます。
いかに特異点が存在する際の数値積分が難しいか分かるでしょう。

解決法

変数変換を行います。
これにより被積分関数を多項式に置き換えることが出来ます。
変数変換により、
$
\displaystyle \int_0^{b}\sqrt{x} dx = \int_0^{b^{1/2}} 2x^2 dx
$
と、簡単になります。
実際、右辺を数値的に計算してみますと、

0.665690422058105 短冊近似　　　(1000分割)
0.666666984558105 台形則　　　　(1000分割)
0.666666666666667 シンプソン則　(1000分割)
0.666666666666667 ロンバーグ則　(1000分割,相対誤差12桁一致)
0.666666666666654 ガウス・ルジャンドル求積(100点)

となります。短冊近似、台形則で精度が悪いのは$x^2$の多項式はこの二つの方法では元々出来ないものなのでこれは問題ありません。数値計算法の想定と一致する結果が導かれています。

プログラミングと数値計算

fortran90によるヒープソートとバブルソート

2016年8月8日 sikino コメントする

fortran90で配列のソートを行うプログラムを載せます。
クイックソート(ここでは載せません)はヒープソートよりも若干遅くなってしまいました。
プログラムの最適化が出来ていない、Fortranの性質（ルーチンの呼び出しが遅い）からかもしれません。

どういうソート方法なのかはwikipediaをご参照ください（ヒープソート、バブルソート）。

計算時間はヒープソートの方が圧倒的に早いです。

クイックソートの方がヒープソートよりも早いです。コードは最速のクイックソートへ。

倍精度型の、サイズNの配列arrayを昇順(小さい数から大きい数へ)に並べ替えます。

バブルソート

平均計算量$O(N^2)$

subroutine bubblesort(N,array)
!sikinote, 2016/08/08
implicit none
integer,intent(in)::N
double precision,intent(inout)::array(1:N)
integer::i,j
double precision::t

do i=1,N-1
do j=i+1,N
if(array(i) .gt. array(j))then
t=array(i)
array(i)=array(j)
array(j)=t
end if
end do
end do

return
end subroutine bubblesort

ヒープソート

平均計算量$O(N\log N)$

subroutine heapsort(n,array)
implicit none
integer,intent(in) :: n
double precision,intent(inout) :: array(1:n)

integer ::i,k,j,l
double precision :: t

if(n.le.0)then
write(6,*)"Error, at heapsort"; stop
endif
if(n.eq.1)return

l=n/2+1
k=n
do while(k.ne.1)
if(l.gt.1)then
l=l-1
t=array(L)
else
t=array(k)
array(k)=array(1)
k=k-1
if(k.eq.1) then
array(1)=t
exit
endif
endif
i=l
j=l+l
do while(j.le.k)
if(j.lt.k)then
if(array(j).lt.array(j+1))j=j+1
endif
if (t.lt.array(j))then
array(i)=array(j)
i=j
j=j+j
else
j=k+1
endif
enddo
array(i)=t
enddo

return
end subroutine heapsort

[adsense1]

配列arrayをどのように入れ替えたのか？という情報が欲しければ

バブルソート

平均計算量$O(N^2)$

subroutine bubblesort2(N,data,turn)
!sikinote, 2016/08/08
integer,intent(in)::N
integer,intent(out)::turn(1:N)
double precision,intent(inout)::data(1:N)
integer::i,j,ti
double precision::tmp

do i=1,N
turn(i)=i
enddo

do i=1,N-1
do j=i+1,N
if(data(i) .gt. data(j))then
tmp=data(i)
data(i)=data(j)
data(j)=tmp

ti=turn(i)
turn(i)=turn(j)
turn(j)=ti
end if
end do
end do

return
end subroutine bubblesort2

ヒープソート

平均計算量$O(N\log N)$

subroutine heapsort2(n,array,turn)
implicit none
integer,intent(in)::n
integer,intent(out)::turn(1:n)
double precision,intent(inout)::array(1:n)

integer::i,k,j,l,m
double precision::t

if(n.le.0)then
write(6,*)"Error, at heapsort"; stop
endif
if(n.eq.1)return

do i=1,N
turn(i)=i
enddo

l=n/2+1
k=n
do while(k.ne.1)
if(l.gt.1)then
l=l-1
t=array(l)
m=turn(l)
else
t=array(k)
m=turn(k)
array(k)=array(1)
turn(k)=turn(1)
k=k-1
if(k.eq.1) then
array(1)=t
turn(1)=m
exit
endif
endif
i=l
j=l+l
do while(j.le.k)
if(j.lt.k)then
if(array(j).lt.array(j+1))j=j+1
endif
if (t.lt.array(j))then
array(i)=array(j)
turn(i)=turn(j)
i=j
j=j+j
else
j=k+1
endif
enddo
array(i)=t
turn(i)=m
enddo

return
end subroutine heapsort2

というプログラムで実装できます。

プログラム例1 ソートする

実際に実行してみます。
乱数で振った配列aを昇順に並べ替えるには

program main
implicit none
integer::i,N
double precision,allocatable::a(:)
double precision::t

N=10
allocate(a(0:N-1)); a=0d0

do i=0,N-1
call random_number(t)
a(i)=t
write(6,*)a(i)
enddo

call heapsort(size(a),a)
do i=0,N-1
write(6,*)i,a(i)
enddo

stop
end program main

subroutine heapsort(n,array)
implicit none
integer,intent(in) :: n
double precision,intent(inout) :: array(1:n)

integer ::i,k,j,l
double precision :: t

if(n.le.0)then
write(6,*)"Error, at heapsort"; stop
endif
if(n.eq.1)return

l=n/2+1
k=n
do while(k.ne.1)
if(l.gt.1)then
l=l-1
t=array(L)
else
t=array(k)
array(k)=array(1)
k=k-1
if(k.eq.1) then
array(1)=t
exit
endif
endif
i=l
j=l+l
do while(j.le.k)
if(j.lt.k)then
if(array(j).lt.array(j+1))j=j+1
endif
if (t.lt.array(j))then
array(i)=array(j)
i=j
j=j+j
else
j=k+1
endif
enddo
array(i)=t
enddo

return
end subroutine heapsort

とすればよいです。実行結果は

$ gfortran main.f90
$ ./a.out
0.99755959009261719
0.56682470761127335
0.96591537549612494
0.74792768547143218
0.36739089737475572
0.48063689875473148
7.3754263633984518E-002
5.3552292777272470E-003
0.34708128851801545
0.34224381607283505
0 5.3552292777272470E-003
1 7.3754263633984518E-002
2 0.34224381607283505
3 0.34708128851801545
4 0.36739089737475572
5 0.48063689875473148
6 0.56682470761127335
7 0.74792768547143218
8 0.96591537549612494
9 0.99755959009261719
$

となります。

[adsense2]

プログラム例2 どのように並べ替えたかの順番を知る

program main
implicit none
integer::i,N
integer,allocatable::turn(:)
double precision,allocatable::a(:)
double precision::t

N=10
allocate(a(0:N-1)); a=0d0
allocate(turn(0:N-1)); turn=0

do i=0,N-1
call random_number(t)
a(i)=t
write(6,*)i,a(i)
enddo
write(6,*)

call heapsort2(size(a),a,turn)
do i=0,N-1
write(6,*)turn(i),a(i)
enddo

stop
end program main

subroutine heapsort2(n,array,turn)
implicit none
integer,intent(in)::n
integer,intent(out)::turn(1:n)
double precision,intent(inout)::array(1:n)

integer::i,k,j,l,m
double precision::t

if(n.le.0)then
write(6,*)"Error, at heapsort"; stop
endif
if(n.eq.1)return

do i=1,N
turn(i)=i
enddo

l=n/2+1
k=n
do while(k.ne.1)
if(l.gt.1)then
l=l-1
t=array(l)
m=turn(l)
else
t=array(k)
m=turn(k)
array(k)=array(1)
turn(k)=turn(1)
k=k-1
if(k.eq.1) then
array(1)=t
turn(1)=m
exit
endif
endif
i=l
j=l+l
do while(j.le.k)
if(j.lt.k)then
if(array(j).lt.array(j+1))j=j+1
endif
if (t.lt.array(j))then
array(i)=array(j)
turn(i)=turn(j)
i=j
j=j+j
else
j=k+1
endif
enddo
array(i)=t
turn(i)=m
enddo

return
end subroutine heapsort2

実行例

$ gfortran main.f90
$ ./a.out
0 0.99755959009261719
1 0.56682470761127335
2 0.96591537549612494
3 0.74792768547143218
4 0.36739089737475572
5 0.48063689875473148
6 7.3754263633984518E-002
7 5.3552292777272470E-003
8 0.34708128851801545
9 0.34224381607283505

8 5.3552292777272470E-003
7 7.3754263633984518E-002
10 0.34224381607283505
9 0.34708128851801545
5 0.36739089737475572
6 0.48063689875473148
2 0.56682470761127335
4 0.74792768547143218
3 0.96591537549612494
1 0.99755959009261719