mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Fortran - 最小二乗法!

[ プログラミング, 数学 ] [ Fortran ]

こんばんは。

前々回、前回は、C++, Ruby による「最小二乗法」のアルゴリズムを紹介しました。

今回は、同じアルゴリズムを Fortran 95 で実現してみました。アルゴリズムについては、上記リンクの記事を参照してください。

以下、Fortran によるサンプルコードです。

0. 前提条件

  • Linux Mint 13 Maya (64bit) での作業を想定。
  • GNU Fortran (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3
  • 最小二乗法についての説明は割愛。(「C++ - 最小二乗法!」を参照)

1. Fortran ソースコード作成

今回作成した Fortran ソースコードは以下の通りです。

least_squares_method.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
!***********************************************************
! 最小二乗法
!***********************************************************
!
!***********************************************************
! 定数モジュール
!***********************************************************
module constants
  implicit none  ! 暗黙の型指定を使用しない
  ! データ数
  integer, parameter :: N = 7
  ! 予測曲線の次数
  integer, parameter :: M = 5
  ! 測定データ
  integer, parameter :: X(7) = (/ -3, -2, -1,  0, 1, 2, 3 /)
  integer, parameter :: Y(7) = (/  5, -2, -3, -1, 1, 4, 5 /)
end module constants

!***********************************************************
! 計算モジュール
!***********************************************************
module calc
  use constants       ! 定数モジュール
  implicit none       ! 暗黙の型指定を使用しない
  private             ! モジュール外非公開
  public :: calc_lsm  ! calc_lsm のみモジュール外公開

contains

  !*********************************************************
  ! 最小二乗法
  !*********************************************************
  subroutine calc_lsm()
    ! 配列定義
    real :: a(M + 2, M + 1)
    real :: s(2 * M + 1) = 0, t(M + 1) = 0

    ! s[], t[] 計算
    call calc_st(s, t)

    ! a[][] に s[], t[] 代入
    call ins_st(a, s, t)

    ! 掃き出し
    call sweap_out(a)

    ! y 値計算&結果出力
    call display(a, s, t)
  end subroutine calc_lsm

  !*********************************************************
  ! 以下、 private subroutine
  !*********************************************************
  ! s[], t[] 計算
  subroutine calc_st(s, t)
    real, intent(inout) :: s(:), t(:)
    integer :: i, j

    do i = 1, N
      do j = 1, 2 * M + 1
        s(j) = s(j) + x(i) ** (j - 1)
      end do
      do j = 1, M + 1
        t(j) = t(j) + x(i) ** (j - 1) * y(i)
      end do
    end do
  end subroutine calc_st

  ! a[][] に s[], t[] 代入
  subroutine ins_st(a, s, t)
    real, intent(inout) :: a(:,:)
    real, intent(in)    :: s(:), t(:)
    integer :: i, j

    do i = 1, M + 1
      do j = 1, M + 1
        a(j, i) = s(i + j - 1)
      end do
      a(M + 2, i) = t(i)
    end do
  end subroutine ins_st

  ! 掃き出し
  subroutine sweap_out(a)
    real, intent(inout) :: a(:,:)
    integer :: i, j, k
    real    :: p, d

    do k = 1, M + 1
      p = a(k, k)
      do j = k, M + 2
        a(j, k) = a(j, k) / p
      end do
      do i = 1, M + 1
        if (i /= k) then
          d = a(k, i)
          do j = k, M + 2
            a(j, i) = a(j, i) - d * a(j, k)
          end do
        end if
      end do
    end do
  end subroutine sweap_out

  ! 結果出力
  subroutine display(a, s, t)
    real, intent(in) :: a(:,:)
    real, intent(in) :: s(:)
    real, intent(in) :: t(:)
    integer :: i, px
    real    :: p

    do i = 1, M +1
      write (*, '(a,i1,a,f10.6)'), "a", i - 1, " = ", a(M + 2, i)
    end do
    write (*, '(a)'), "    x    y"
    do px = -3 * 2, 3 * 2
      p = 0
      do i = 1, M + 1
        p = p + a(M + 2, i) * ((px / 2.0) ** (i - 1))
      end do
      write (*, '(f5.1,f5.1)'), (px / 2.0), p
    end do
  end subroutine display
end module calc

!***********************************************************
! 主処理
!***********************************************************
program least_squares_method
  use constants  ! 定数モジュール
  use calc       ! 計算モジュール
  implicit none  ! 暗黙の型指定を使用しない

  ! ==== 最小二乗法
  call calc_lsm()
end program least_squares_method

2. Fortran ソースコードコンパイル

以下のコマンドで Fortran ソースコードをコンパイルする。

1
$ f95 least_squares_method.f95 -o least_squares_method

何も出力されなければ成功。

3. 実行

実際に、実行してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
$ ./least_squares_method
a0 =  -1.259740
a1 =   2.100002
a2 =   0.424242
a3 =  -0.083334
a4 =   0.030303
a5 =  -0.016667
    x    y
 -3.0  5.0
 -2.5  0.3
 -2.0 -2.1
 -1.5 -2.9
 -1.0 -2.8
 -0.5 -2.2
  0.0 -1.3
  0.5 -0.1
  1.0  1.2
  1.5  2.6
  2.0  3.9
  2.5  4.9
  3.0  5.0

C++ 版、Ruby 版と同じ結果になるはず。


数学的な内容だったので Fortran 95 でも実装してみた次第です。

以上。

Comments