1

この質問をするのにこの適切な部屋があるかどうかはわかりません。そうでない場合は、申し訳ありません。私はFortranの新規ユーザーであり、次のことに多くの時間を費やしています。

2つのパラメータに応じて実数を返す「loglike」という関数を作成しました。この関数を使用してmcmcアルゴリズムを構築したい

これはこのようになります。

psi = min(0, loglike(propalpha,propbeta) - loglike(currentalpha,currentbeta))

ここpropalpha = currentalpha + noiseでpropbeta = currentbeta + noise、、、、のノイズは、ある分布からのランダムなサンプルです。

ここで、以前に作成した関数「loglike」を呼び出して、このアルゴリズムを使用したいと思います。

1）メインプログラムと呼ばれる新しいプログラムの関数「loglike」を呼び出すには
どうすればよいですか？2）サブルーチンにこれを使用するにはどうすればよいですか？

どんな助けも私にとってとても素晴らしいです。前もって感謝します

編集：

module mcmc

implicit none

contains

    subroutine subru(A,B, alphaprop, betaprop)
        real, intent(in)::A,B
        real, intent(out)::alphaprop, betaprop
    end subroutine subru

    real function aloglike(A,B)
        real:: A,B,U, aloglike
        aloglike = U
    end function aloglike

end module mcmc


program  likelihood
    use mcmc

    implicit none

    real :: alpha,beta,dist1,dist2,prob1,prob2
    real:: psi,a(1000),b(1000), u1, u2,u, alphaprop, betaprop
    real, dimension(1:1):: loglike
    integer :: t,i,j,k,l
    real, dimension(1:625):: x
    real, dimension(1:625):: y
    integer, dimension(1:625):: inftime

    alpha = 0.5
    beta = 2.0

    open(10, file = 'epidemic.txt', form = 'formatted')
    do l = 1,625
       read(10,*,end = 200) x(l), y(l), inftime(l)
    enddo
200 continue

    loglike = 0.0
    do t=1,10
        do i=1,625
            if(inftime(i)==0 .or. t < (inftime(i)-1)) then
                dist1 = 0.0
                do  j = 1, 625
                    if(t >= inftime(j) .and. inftime(j)/=0)then
                        dist1 = dist1 + sqrt((x(i) - x(j))**2 + (y(i) - y(j))**2)**(-beta)
                    endif
                enddo
                prob1 = 1 - exp(-alpha * dist1)
                loglike = loglike + log(1 - prob1)
            endif

            if(inftime(i) .eq. (t+1)) then
                dist2 = 0.0
                    do k=1, 625
                    if(t >= inftime(k) .and. inftime(k)/=0) then
                        dist2 = dist2 + sqrt((x(i) - x(k))**2 + (y(i) - y(k))**2)**(-beta)
                    endif
                enddo
                prob2 = 1 - exp(-alpha * dist2)
                loglike = loglike + log(prob2)
            endif
        enddo
    enddo

    do i = 2, 1000
        a(1)= 0.0
        b(1) = 0.0
        call subru(a(i),b(i), alphaprop, betaprop)
        call random_number(u1)
        call random_number(u2)

        alphaprop = a(i-1) + (u1*0.4)-0.2
        betaprop= b(i-1) + (u2*0.4)-0.2

        if(alphaprop> 0 .and. alphaprop < 0.2 .and. betaprop > 0 .and. betaprop < 0.2)then
            psi = min(0.0,aloglike(alphaprop,betaprop)- aloglike(a(i-1),b(i-1)))
            call random_number(u)

            if(u < psi)then
                a(i)= alphaprop
                b(i) = betaprop
            else
                a(i) = a(i-1)
                b(i) = b(i-1)
            endif
        endif
    enddo

    do j = 1, 1000
        print *, A(j), A(j), LOGLIKE
    enddo

end program
4

2 に答える 2

4

最も簡単で信頼性の高い手法は、関数とサブルーチンをモジュールに配置し、メイン プログラムからそのモジュールを "使用" することです。これは、1 つのファイルで実行できます。このメソッドは、プロシージャー (関数およびサブルーチン) のインターフェースを認識できるようにするため、コンパイラーは、呼び出し (実引数) と呼び出される (ダミー引数) の引数の間の整合性をチェックできます。スケッチ：

module mysubs

implicit none

contains

subroutine sub1 (xyz)
   declarations
   code
end subroutine sub1

function func2 (u)
   declarations
   code
   func2 = ...
end func2

end module mysubs

program myprog

use mysubs

implicit none

declarations...

call sub1 (xyz)

q = func2 (z)

end program myprog

追加: 「implicit none」は暗黙の型付けを無効にするために使用されますが、これは私の意見では危険です。そのため、関数内の関数名を含め、すべての変数を入力する必要があります。モジュールの他のプロシージャからサブルーチンと関数を呼び出すことができます。それらは自動的に認識されます。したがって、必要に応じて、「sub1」から「func2」を使用できます。メイン プログラムなど、モジュールの外部にあるエンティティについては、モジュールを「使用」する必要があります。

于 2012-02-17T22:18:00.513 に答える
2

これが一般的な外観です。関数は、結果を独自の名前に割り当てることによって値を返すことに注意してください。returnステートメントは必要ありません。

C "loglike" is renamed "aloglike" so implicit typing uses REAL
C A type declaration could be used instead

    function aloglike (alpha, beta)
        aloglike = ... (expression which computes value)
    end

    program whateveryouwanttocallit
       ...
       propalpha = ...
       propbeta = ...
       currentalpha = ...
       currentbeta = ...
       ...
       psi = min(0, aloglike(propalpha,propbeta) - aloglike(currentalpha,currentbeta))
       print *, 'psi = ', psi

    end
于 2012-02-17T21:03:14.233 に答える