將fortran程序轉換爲從R調用的子例程

我以前發佈了一個關於試圖從統計報告中重現樣本大小的問題here的問題。作者向我提供了Fortran程序中使用的代碼。我想在工作中使用我的團隊使用的R包來調用此代碼。爲此，我需要將此Fortran程序轉換爲子程序。這是我第一次使用Fortran，並且我似乎以某種方式破壞了該程序，並且在更改輸入時弄錯了不會改變的數字。我還應該注意到，我必須從作者提供的原始代碼中做一些小改動才能使程序運行。將fortran程序轉換爲從R調用的子例程

「原始」程序編譯時產生正確的結果，然後跑：

PROGRAM Exact_Mid_P_binomial_sample_size 
      real(8) probA,probB,part1,part2,part3,part4 
      real(8) totprA,totprB,factt, resp 
!  character resp 
1  format ('Enter proportion  ',$) 
2  format ('Enter error limit  ',$) 
3  format ('Enter confidence level ',$) 
4  format ('Calculated sample size is ',i6) 
5  format ('Exact mid-P with ',f7.5,' 2-tail probability') 
6  format ('Sorry, unable to mathmatically solve this problem.') 
7  format ('Reported sample size is not accuarate.') 
8  format ('Enter q to quit ',$) 
9  format ('Actual limits for distribution ',f5.3,' - ',f5.3) 
     print *, 'Exact sampleroportions' 
     print *, 'Using Mid-P methods' 
     print *, 'Geoff Fosgate DVM PhD' 
     print *, 'College of Veterinary Medicine' 
     print *, 'Texas A&M University' 
     print * 
10  print * 
     print 1 
      read *, prop1 
     print 2 
      read *,range 
     print 3 
      read *,conlev 
     print * 
!   Convert proportions less than 0.5 for algorithm 
     if (prop1 .lt. 0.5) then 
      prop = 1 - prop1 
      nprop = 1 
       else 
       prop = prop1 
       nprop = 0 
       end if 
     slimit = max ((prop - range) , 0.0001) 
     supper = min ((prop + range) , 0.9999) 
!   Probabilities cannot be calculated for p=0 and p=1 
     alpha = (1 - conlev) 
     if (alpha .gt. 1.0) go to 10 
     if (alpha .lt. 0.0) go to 10 
     if (prop .gt. 1.0) go to 10 
     if (prop .lt. 0.0) go to 10 
     numbr = (1/(1 - prop)) - 1 
!   Define and initialize variables 
!    Note names of variables based on Fortran 77 rules 
!    Starting sample size is based on estimated proportion 
!    Resulting sample size must be large enough to obtain this proportion 
100  numbr = numbr + 1 
     numx = (numbr * prop) + 0.001 
!    This is the number of binomial "successes" resulting in the proportion 
     if (numx .eq. numbr) go to 100 
     if (numx .lt. 1) go to 100 
     totprA = slimit**numbr 
     totprB = supper**numbr 
      do 130 loop1 = numx, (numbr - 1) 
!    Must initialize variables within loop 
      factt = 1.0 
      probA = 0.0 
      probB = 0.0 
      part1 = 0.0 
      part2 = 0.0 
      part3 = 0.0 
      part4 = 0.0 
!    Start loop to calculate factorial component of binomial probability 
!    Note that complete factorial calculations not necessary due to cancellations 
      do 110 loop2 = (loop1 + 1) , numbr 
       factt = factt * (loop2)/(numbr - (loop2 - 1)) 
110   continue 
!    Calculate probability for this particular number of successes 
!    Total probability is a running total 
!    Note that real variables must have high precision and be comprised 
!    of multiple bytes because factorial component can be very large 
!    and exponentiated component can be very small 
!    Program will fail if any component is recognized as zero or infinity 
      part1 = slimit**loop1 
       part2 = (1.0-slimit)**(numbr-loop1) 
      part3 = supper**loop1 
       part4 = (1.0-supper)**(numbr-loop1) 
      if (part1 .eq. 0.0) part1 = 1.0D-307 
      if (part2 .eq. 0.0) part2 = 1.0D-307 
      if (part3 .eq. 0.0) part3 = 1.0D-307 
      if (part4 .eq. 0.0) part4 = 1.0D-307 
      if (factt .gt. 1.0D308) factt = 1.0D308 
      probA = part1 * part2 * factt 
      probB = part3 * part4 * factt 
      if (loop1 .eq. numx) then 
       totprA = totprA + (0.5 * probA) 
       totprB = totprB + (0.5 * probB) 
       else 
        totprA = totprA + probA 
        totprB = totprB + probB 
       end if 
      if (probA .eq. 0.0) then 
        print 6 
        print 7 
        print * 
        go to 150 
       end if 
      if (probB .eq. 0.0) then 
        print 6 
        print 7 
        print * 
        go to 150 
       end if 
130  continue 
140  if ((totprA + (1 - totprB)) .gt. alpha) go to 100 
!    go to beginning and increase sample size by 1 if have not 
!    reached specified level of confidence 
150   if (nprop .eq. 1) then 
       print 4,numbr 
       print 9, (1-supper),(1-slimit) 
      else 
       print 4,numbr 
       print 9, slimit,supper 
      end if 
      if (totprA+(1-totprB) .lt. alpha) print 5,(totprA+(1-totprB)) 
      print * 
      print 8 
      result = resp 
!   print * 
!  if (resp .ne. 'q') go to 10 
     print * 
     print * 
999 end

上面變成一個子程序的程序：

 subroutine midpss(prop1, range, conlev, numbr) 

     integer numbr, nprop 
     real(8) prop1, range, conlev, prop, slimit, supper 
     real(8) probA,probB,part1,part2,part3,part4,factt 
     real(8) totprA,totprB, resp 
c   character resp 

c   Convert proportions less than 0.5 for algorithm 
      if (prop1 .lt. 0.5) then 
       prop = 1 - prop1 
       nprop = 1 
      else 
       prop = prop1 
       nprop = 0 
      end if 
     slimit = max ((prop - range) , 0.0001) 
     supper = min ((prop + range) , 0.9999) 

      numbr = (1/(1 - prop)) - 1 
c   Define and initialize variables 
c    Note names of variables based on Fortran 77 rules 
c    Starting sample size is based on estimated proportion 
c    Resulting sample size must be large enough to obtain this proportion 
100  numbr = numbr + 1 
      numx = (numbr * prop) + 0.001 
c    This is the number of binomial "successes" resulting in the proportion 
      if (numx .eq. numbr) go to 100 
      if (numx .lt. 1) go to 100 
      totprA = slimit**numbr 
      totprB = supper**numbr 
      do 130 loop1 = numx, (numbr - 1) 
c    Must initialize variables within loop 
      factt = 1.0 
      probA = 0.0 
      probB = 0.0 
      part1 = 0.0 
      part2 = 0.0 
      part3 = 0.0 
      part4 = 0.0 
c    Start loop to calculate factorial component of binomial probability 
c    Note that complete factorial calculations not necessary due to cancellations 
     do 110 loop2 = (loop1 + 1) , numbr 
     factt = factt * (loop2)/(numbr - (loop2 - 1)) 
110 continue 
c    Calculate probability for this particular number of successes 
c    Total probability is a running total 
c    Note that real variables must have high precision and be comprised 
c    of multiple bytes because factorial component can be very large 
c    and exponentiated component can be very small 
c    Program will fail if any component is recognized as zero or infinity 
      part1 = slimit**loop1 
     part2 = (1.0-slimit)**(numbr-loop1) 
      part3 = supper**loop1 
     part4 = (1.0-supper)**(numbr-loop1) 
      if (part1 .eq. 0.0) part1 = 1.0D-307 
      if (part2 .eq. 0.0) part2 = 1.0D-307 
      if (part3 .eq. 0.0) part3 = 1.0D-307 
      if (part4 .eq. 0.0) part4 = 1.0D-307 
      if (factt .gt. 1.0D308) factt = 1.0D308 
     probA = part1 * part2 * factt 
      probB = part3 * part4 * factt 
      if (loop1 .eq. numx) then 
       totprA = totprA + (0.5 * probA) 
       totprB = totprB + (0.5 * probB) 
      else 
       totprA = totprA + probA 
       totprB = totprB + probB 
      end if 

130  continue 
140  if ((totprA + (1 - totprB)) .gt. alpha) go to 100 


     return 
     end

要編譯原來的程序，我運行這個在命令行中：

gfortran midpSS_original.f

要編譯子程序，我在c中運行這個ommand行：

R CMD SHLIB midpSS_subroutine.f

，然後從R控制檯，我運行此：

> dyn.load("midpSS_subroutine.dll") 
> is.loaded("midpss") 
[1] TRUE 
> .Fortran("midpss", prop1=as.numeric(0.9), range=as.numeric(0.1), conlev=as.numeric(0.90), numbr=as.integer(0)) # numbr should be 29 
$prop1 
[1] 0.9 

$range 
[1] 0.1 

$conlev 
[1] 0.9 

$numbr 
[1] 2091 

> .Fortran("midpss", prop1=as.numeric(0.9), range=as.numeric(0.1), conlev=as.numeric(0.95), numbr=as.integer(0)) # numbr should be 47 
$prop1 
[1] 0.9 

$range 
[1] 0.1 

$conlev 
[1] 0.95 

$numbr 
[1] 2091

在第一次調用的子程序，所述numbr應該是29。在第二個，numbr應47.編譯「原始」Fortran程序並運行相同的參數時，結果是正確的。我無法弄清楚這裏發生了什麼。任何幫助將不勝感激。

來源

2016-01-05 panterasBox

你能找到一種方法將編譯選項傳遞給R嗎？讓編譯器警告你說，你沒有在你的子程序中設置'alpha'會很有幫助。 – francescalus

我建議你寫一個新的小Fortran main調用你的新子程序。確保首先工作，然後在'R'界面上工作。 – agentp

@agentp我只是做了這個。當通過R調用它時，我得到的輸出錯誤...感謝。我現在*知道*我把它轉換成子程序時搞砸了。 – panterasBox

以下塊檢查輸入數據：

if (alpha .gt. 1.0) go to 10 
if (alpha .lt. 0.0) go to 10 
if (prop .gt. 1.0) go to 10 
if (prop .lt. 0.0) go to 10

此塊使代碼去入口位置再次讀取輸入數據。您可能會得到錯誤的結果，因爲輸入到子例程中的數據可能超出範圍。請先檢查一下！

另一個丟失塊是

if (probA .eq. 0.0) then 
    print 6 
    print 7 
    print * 
go to 150 
end if 

if (probB .eq. 0.0) then 
print 6 
print 7 
print * 
go to 150 
end if

和標籤150爲好。你確定省略這兩個塊沒有改變原代碼的準確性和功能嗎？

將程序轉換爲子程序非常簡單。主要步驟如下：

你的子程序
的把名字後轉換的語法program爲subroutine
轉換語法end program到end subroutine
將一個括號「（）」的所有輸入和輸出程序中的括號'（）'

這也是一個好習慣implicit none後program或subroutine

我將程序轉換成我自己的子程序。請注意，子例程末尾的result = resp可能是錯誤的，因爲據我所知，命令'result'是一個函數而不是子例程。

在子程序中你必須指定輸出。由於我不知道哪個變量是輸出，因此我沒有在子程序名稱後的括號'（）中包含該變量。請在那裏寫下。我的子程序如下所示：

! PROGRAM Exact_Mid_P_binomial_sample_size 
    subroutine midpss (prop1, range, conlev, output) 
    implicit none 
      real(8) probA,probB,part1,part2,part3,part4 
      real(8) totprA,totprB,factt, resp 
      integer numbr, nprop,loop1,loop2 
      real(8) prop1,prop,range,conlev,slimit,supper,alpha,numx,output 

!  character resp 
!1  format ('Enter proportion  ',$) 
!2  format ('Enter error limit  ',$) 
!3  format ('Enter confidence level ',$) 
4  format ('Calculated sample size is ',i6) 
5  format ('Exact mid-P with ',f7.5,' 2-tail probability') 
6  format ('Sorry, unable to mathmatically solve this problem.') 
7  format ('Reported sample size is not accuarate.') 
8  format ('Enter q to quit ',$) 
9  format ('Actual limits for distribution ',f5.3,' - ',f5.3) 
     print *, 'Exact sampleroportions' 
     print *, 'Using Mid-P methods' 
     print *, 'Geoff Fosgate DVM PhD' 
     print *, 'College of Veterinary Medicine' 
     print *, 'Texas A&M University' 
     print * 
!10  print * 
!  print 1 
!   read *, prop1 
!  print 2 
!   read *, range 
!  print 3 
!   read *, conlev 
!  print * 
!   Convert proportions less than 0.5 for algorithm 
     if (prop1 .lt. 0.5) then 
      prop = 1 - prop1 
      nprop = 1 
       else 
       prop = prop1 
       nprop = 0 
       end if 
     slimit = max ((prop - range) , 0.0001) 
     supper = min ((prop + range) , 0.9999) 
!   Probabilities cannot be calculated for p=0 and p=1 
     alpha = (1 - conlev) 
!  if (alpha .gt. 1.0) go to 10 
!  if (alpha .lt. 0.0) go to 10 
!  if (prop .gt. 1.0) go to 10 
!  if (prop .lt. 0.0) go to 10 
     numbr = (1/(1 - prop)) - 1 
!   Define and initialize variables 
!    Note names of variables based on Fortran 77 rules 
!    Starting sample size is based on estimated proportion 
!    Resulting sample size must be large enough to obtain this proportion 
100  numbr = numbr + 1 
     numx = (numbr * prop) + 0.001 
!    This is the number of binomial "successes" resulting in the proportion 
     if (numx .eq. numbr) go to 100 
     if (numx .lt. 1) go to 100 
     totprA = slimit**numbr 
     totprB = supper**numbr 
      do 130 loop1 = numx, (numbr - 1) 
!    Must initialize variables within loop 
      factt = 1.0 
      probA = 0.0 
      probB = 0.0 
      part1 = 0.0 
      part2 = 0.0 
      part3 = 0.0 
      part4 = 0.0 
!    Start loop to calculate factorial component of binomial probability 
!    Note that complete factorial calculations not necessary due to cancellations 
      do 110 loop2 = (loop1 + 1) , numbr 
       factt = factt * (loop2)/(numbr - (loop2 - 1)) 
110   continue 
!    Calculate probability for this particular number of successes 
!    Total probability is a running total 
!    Note that real variables must have high precision and be comprised 
!    of multiple bytes because factorial component can be very large 
!    and exponentiated component can be very small 
!    Program will fail if any component is recognized as zero or infinity 
      part1 = slimit**loop1 
       part2 = (1.0-slimit)**(numbr-loop1) 
      part3 = supper**loop1 
       part4 = (1.0-supper)**(numbr-loop1) 
      if (part1 .eq. 0.0) part1 = 1.0D-307 
      if (part2 .eq. 0.0) part2 = 1.0D-307 
      if (part3 .eq. 0.0) part3 = 1.0D-307 
      if (part4 .eq. 0.0) part4 = 1.0D-307 
      if (factt .gt. 1.0D308) factt = 1.0D308 
      probA = part1 * part2 * factt 
      probB = part3 * part4 * factt 
      if (loop1 .eq. numx) then 
       totprA = totprA + (0.5 * probA) 
       totprB = totprB + (0.5 * probB) 
       else 
        totprA = totprA + probA 
        totprB = totprB + probB 
       end if 
      if (probA .eq. 0.0) then 
        print 6 
        print 7 
        print * 
        go to 150 
       end if 
      if (probB .eq. 0.0) then 
        print 6 
        print 7 
        print * 
        go to 150 
       end if 
130  continue 
140  if ((totprA + (1 - totprB)) .gt. alpha) go to 100 
!    go to beginning and increase sample size by 1 if have not 
!    reached specified level of confidence 
150   if (nprop .eq. 1) then 
       print 4,numbr 
       print 9, (1-supper),(1-slimit) 
      else 
       print 4,numbr 
       print 9, slimit,supper 
      end if 
      if (totprA+(1-totprB) .lt. alpha) print 5,(totprA+(1-totprB)) 
      print * 
      print 8 
      !result = resp 
!   print * 
!  if (resp .ne. 'q') go to 10 
!  output=***  !write the output variable instead of *** 
     print * 
     print * 
999 end subroutine midpss

來源

2016-01-06 07:21:36 Fortraneer

關於'結果= resp'， 'result'只是一個變量的名稱，並且與函數的結果子句無關。 [這也是一個變量，在這個任務之外似乎沒有被引用，但這是另一個問題。] – francescalus

'隱式無作用'是什麼？ – panterasBox

我用你的第一個建議檢查輸入數據。當我從一個fortran程序中調用子程序時（現在它似乎被掛起來調用與R相同的子程序），現在就可以工作了。謝謝。 – panterasBox

將fortran程序轉換爲從R調用的子例程

回答

相關問題