Question

I don't know is this right room to ask this question of not. If not I am sorry for that. I am new user for the fortran and spending a lot of time for the following stuff.

I have constructed a function called "loglike" which returns a real number depending on two parameters. I want to use this function to construct a mcmc algorithm

which goes like this.

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

where propalpha = currentalpha + noise, and propbeta = currentbeta + noise, noises are random samples from some distribution.

Now I want to use this algorithm by calling previously constructed function "loglike".

1) how can I call the function 'loglike' for new program called main program
2) how can I use this for the subroutine?

Any help is very great for me. Thanks in advance

EDIT:

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
Was it helpful?

Solution

This is the general way it would look. Note that functions return a value by assigning the result to its own name. A return statement is not necessary.

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

OTHER TIPS

The easiest and most reliable technique is to place your functions and subroutines in a module and "use" that module from your main program. This can be done in one file. This method makes the interfaces of the procedures (functions and subroutines) known so that the compiler can check consistency between arguments in the call (actual arguments) and called (dummy arguments). Sketch:

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

ADDED: "implicit none" is used to disable implicit typing, which is dangerous in my opinion. So you will need to type all of your variables, including the function name in the function. You can call the subroutines and functions from other procedures of the module -- they will automatically be known. So you can use "func2" from "sub1", if you wish. For entities outside of the module, such as your main program, you must "use" the module.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top