平均、分散、標準偏差

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


 例1 : 40人クラスにおいてテスト点数の解析を行いたい(40人全員の点数を知っている)。
 例2 : 過去1年間、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)一部分のデータしか存在しない場合


例1 : 40人クラスにおいて23人しかテストの点数を知らないが、その23人からクラス全体の点数の解析を行いたい(23人の点数は分かる。残り17人の点数は不明)。
 例2 : 過去1年間、365日で体重の増減を解析したいが、200日分のデータしかない。だけど1年の解析を行いたい(165日分のデータは不明、紛失した)。
    例3 : 物理量の測定(測定回数が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の数値積分(端点で特異点が存在する場合)

\(\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によるヒープソートとバブルソート

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