匿名希望のおでんFortranツヴァイさん太郎

生き物、Fortran、川について書く

分散の大きさの感覚を養う

Fortran 統計学

何事も実際にやってみるのは大切です．経験することで直観が養われます．

本記事は分散の大きさを感覚的に理解することが目的です．分散は平均とともに統計学で最初に習うものですから，定義を知っている人は多いでしょう，しかし，散布図を見て分散がどのくらいかを感覚的に答えられる人は少ないのではないでしょうか．

実験してみましょう．コンパイラはgfortran4.9.3です．

xは [0, 2\times3.1415)の一様乱数から発生させます．Fortranの組み込み手続きであるrandom_numberを用います．なお，random_numberで用いられている [0, 1)からの一様乱数の発生アルゴリズムは，gfortranではxorshift1024*です*1．今回はxは別に乱数でなくてもよいのですが，ランダムにとったデータぽくみせる演出です．

yは平均\sin x，分散\sigma^2の正規乱数から発生させます．正規乱数のアルゴリズムはZiggurat法です．実装方法は四辻（2010）*2を参照してください．

\sigma^2= 1.0, 0.1, 0.01, 0.001の四つの場合で結果を比較しましょう．


f:id:y_tis:20181204222216j:plain:w400
\sigma^2= 1.0


f:id:y_tis:20181204222222j:plain:w400
\sigma^2= 0.1


f:id:y_tis:20181204222258j:plain:w400
\sigma^2= 0.01


f:id:y_tis:20181204222315j:plain:w400
\sigma^2= 0.001

正規分布なら各xにおいて\sin x \pm 2\sigmaにそのxにおけるデータの約95%が入ります．求めるものにもよりますが，分散1.0というのは結構大きく感じます．精密な較正曲線が求められる場面でこんなグラフを出したら怒られそうです．

最後に，今回用いたコードを以下に示します．

program sin_normal
implicit none
double precision x(100)
double precision y_1
double precision y_01
double precision y_001
double precision y_0001
integer n
integer i

n = size(x)

open(10, file="sin_normal.csv" , status="replace")
write(10, *) "x, y_1.0, y_0.1, y_0.01, y_0.001"
call random_number(x) 
x = 2.0d0 * 3.1415d0 * x
do i = 1, n
	y_1 = sin(x(i)) + random_normal(0.0d0, 1.0d0)
	y_01 = sin(x(i)) + random_normal(0.0d0, 1.0d-1)
	y_001 = sin(x(i)) + random_normal(0.0d0, 1.0d-2)
	y_0001 = sin(x(i)) + random_normal(0.0d0, 1.0d-3)
	write(10, *) x(i), ",", y_1, ",", y_01, ",", y_001, ",", y_0001
end do
close(10)

contains

function random_normal(mu, sigma2) result(rslt)
	double precision, intent(in)::mu
	double precision, intent(in)::sigma2 
	integer, parameter::k = 8
	integer, parameter::n = 2**k
	double precision, parameter::r = 3.6541528853610088d0
	double precision, parameter::v = 0.00492867323399d0
	double precision x(0:n)
	double precision y
	double precision u1
	double precision u_1
	double precision u__1
	double precision u___1
	double precision ux
	double precision u2
	double precision gu
	double precision gl
	integer sgn
	integer i
	integer j
	integer b
	integer l
	double precision rslt
	
	x(n) = v * exp(r * r * 0.5d0)
	x(n-1) = r
	do l = n-2, 1, -1
		x(l) = sqrt(-2.0d0 * log(exp(-x(l+1) * x(l+1) * 0.5d0) + v / x(l+1)))
	end do
	x(0) = 0.0d0
	
	do 
		!    ステップ 1
		i = 0
		j = 0
		b = 1
		call random_number(u1)
		
		do 
			!    ステップ2
			u_1 = u1 + u1
			if (u_1 < 1.0d0) then
				i = i + b
				u__1 = u_1
			else if(u_1 >= 1.0d0) then
				u__1 = u_1 - 1.0d0
			end if
			!    ステップ3
			j = j + 1
			if (j < k) then
				b = b + b
				u1 = u__1
				cycle ! ステップ2へ戻る
			end if
			exit !ステップ4へ進む
		end do

		!    ステップ4
		u___1 = u__1 + u__1
		if (u___1 < 1.0d0) then
			sgn = 1
			ux = u___1 * x(i+1)
		else
			sgn = -1
			ux = (u___1 - 1.0d0) * x(i+1)
		end if
		
		!    ステップ5
		if (ux < x(i)) then
			y = ux * dble(sgn)
			exit !ステップ7へ進む
		end if
		
		!    ステップ6-a
		if (i == n - 1) then
			y = random_normal_tail(r) * dble(sgn)
			exit !ステップ7へ進む
		!    ステップ6-b1
		else
			gu = exp(-(x(i) * x(i) - ux * ux) * 0.5d0)
			gl = exp(-(x(i+1) * x(i+1) - ux * ux) * 0.5d0)
			!   ステップ6-b2
			call random_number(u2)
			if (u2 * (gu - gl) <= 1.0d0 - gl) then
				y = ux * dble(sgn)
				exit !ステップ7へ進む
			end if
			cycle !ステップ1に戻る
		end if
	end do
	
	!    ステップ7
	rslt = mu + sqrt(sigma2) * y

end function

function random_normal_tail(b) result(rslt)
	double precision, intent(in)::b
	double precision d
	double precision x
	double precision u1
	double precision u2
	double precision rslt
	
	!    ステップ1
	d = b * b
	do
		!    ステップ2
		call random_number(u1)
		call random_number(u2)
		x = sqrt(d - (log(1.0d0 - u1) + log(1.0d0 - u1)))
		!    ステップ3
		if (x * u2 <= b) exit
		cycle
	end do
	rslt = x
end function

end program

*1:詳しいs仕様は以下のリンクを参照してください． gcc.gnu.org

*2:四辻哲章（2010）計算機シミュレーションのための確率分布乱数生成法．プレアデス出版．

計算機シミュレーションのための確率分布乱数生成法

計算機シミュレーションのための確率分布乱数生成法