Reputation: 33
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
Upvotes: 1
Views: 306
Reputation: 29381
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.
Upvotes: 4
Reputation: 57774
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
Upvotes: 2