如何使 fortran 子程序结果在 R 上可重现？

问：

我编写了一个 Fortran 子例程：

我想在 R 上实现它，所以我将进行编译和所有操作，但我有一个问题，因为我的 fortran 子例程包含随机数生成过程。

这是我的 fortran 子程序，我使用 [https://jblevins.org/mirror/amiller/random.f90][1] 中的 random.f90 和 [https://people.sc.fsu.edu/~jburkardt/f_src/cdflib/cdflib.html][2] 中的 cdflib

我在 transform_samples_module中使用 random.f90（使用来自 random.f90 的随机正态样本，将它们转换为生成具有特定均值和 sd 的正态样本）和来自 cdflib 的 cdfnor 子例程 getprob1 子例程

subroutine getprob1(n1, mu0, mu1, cv0, cv1, nsim, prob, seed) bind(C, name="getprob1")
    use iso_c_binding
    use random
    use transform_samples_module
    
    implicit none
    ! ---- arg types -----------------------
    integer(c_int), intent(in) :: n1, nsim
    real(c_double), intent(in) :: mu0, mu1, cv0, cv1
    real(c_double), intent(out) :: prob(nsim)
    integer(c_int), intent(in) :: seed

    ! ---------------------------------------
  

    real(c_double) :: sigma0, sigma1
    real(c_double) :: mu_2arm(2), sigma_2arm(2), ybar_2arm(2), var_2arm(2)
    real(c_double), dimension(n1) :: nsamples 
    real(c_double), dimension(n1, 2) :: sample_j
    integer(c_int) :: ite, jj, ii
    real(c_double) :: p, q, x, bound, mean, sd
    integer(c_int) :: status, which
    real(c_double) :: cdf(2), delta(2)

    !!!!!!!!!!!!!!!!!

    integer(c_int), dimension(8) :: random_seed_seq

    ! Set the seed
    random_seed_seq = [seed, seed, seed, seed, seed, seed, seed, seed]
    call random_seed(PUT=random_seed_seq)

        
    
    sigma0 = sqrt(log(cv0*cv0 + 1))
    sigma1 = sqrt(log(cv1*cv1 + 1))
    mu_2arm = [mu0, mu1]
    sigma_2arm = [sigma0, sigma1]
    
    do ite = 1, nsim
        
        ybar_2arm = 0.0
        var_2arm = 0.0
  
        
        do jj = 1, 2

    
            call transform_samples(n1, mu_2arm(jj), sqrt(sigma_2arm(jj)), nsamples)
            sample_j(:,jj ) = nsamples    
            ybar_2arm(jj ) = sum(sample_j(:,jj )) / real(n1)
            var_2arm(jj ) = sum((sample_j(:,jj ) - ybar_2arm(jj) )**2) / real(n1 - 1)
          
           
        end do   
        which = 1
        delta = [0.223 , -0.223]
        mean = ybar_2arm(1) - ybar_2arm(2)
        sd = sqrt(var_2arm(1) / real(n1) + var_2arm(2) / real(n1))
        cdf = 0.0
        

        do ii = 1, 2
            x = delta(ii)
            call cdfnor ( which, p, q, x, mean, sd, status, bound )
            if (status /= 0) then
                print *, "Error in cdfnor: status=", status
                stop
            end if
            cdf(ii) = p 
    
        end do
    
        prob(ite) = cdf(1) - cdf(2)


    end do    
    

end subroutine getprob1

这是我transform_samples_module的样子：

module transform_samples_module
    use iso_c_binding
    use random
    implicit none

contains

    subroutine transform_samples(n_samples, desired_mean, desired_stddev, output_samples) bind(C, name="transform_samples")
        integer(c_int), intent(in) :: n_samples
        real(c_double), intent(in) :: desired_mean, desired_stddev
        real(c_double), dimension(n_samples), intent(out) :: output_samples
        real(c_double) :: random_numbers(n_samples)
        integer(c_int) :: i
        
        call random_number(random_numbers)
        random_numbers = sqrt(-2.0d0 * log(random_numbers)) * cos(2.0d0 * 3.141592653589793d0 * random_numbers)
        
        do i = 1, n_samples
            output_samples(i) = desired_mean + desired_stddev * random_numbers(i)
        end do
    end subroutine transform_samples

end module transform_samples_module

所以我在R上进行了编译，使用：

library(dotCall64)
# Load the shared library
dyn.load("mylibrary.so")
getprob1 <- function(n1, mu0, mu1, cv0, cv1, nsim, seed) {
  seed <- .Random.seed[-1]  # Exclude the first element of the R random seed
  .C64("getprob1",
       SIGNATURE = c("integer", "double", "double", "double", "double", "integer", "double", "integer"),
       as.integer(n1),
       as.double(mu0),
       as.double(mu1),
       as.double(cv0),
       as.double(cv1),
       as.integer(nsim),
       double(nsim),
       as.integer(seed)
  )[[7]]  # This assumes that the `prob` array is the 7th argument in your Fortran subroutine
}

好消息是它可以运行！但我想做的是，我希望它能够在 R 中使用种子，例如，如果我在 R 中设置了 seed（1），它应该能够重现相同的结果。因为现在我只能这样做：

> set.seed(1)
> getprob1(34, 1.0, 0.777, 0.5, 0.5, 10, seed = 1)
 [1] 0.8458423 0.7302202 0.6892795 0.4100851 0.7436128 0.5248045 0.5963506 0.6377372 0.2405343 0.3257491
> getprob1(34, 1.0, 0.777, 0.5, 0.5, 10, seed = 1)
 [1] 0.05997150 0.37370831 0.67539986 0.28370543 0.46473538 0.77224209 0.01801777 0.60191086 0.69541114 0.69710732
> set.seed(1)
> getprob1(34, 1.0, 0.777, 0.5, 0.5, 10, seed = 1)
 [1] 0.44087374 0.07380424 0.77093088 0.19258925 0.10651923 0.31367524 0.69142104 0.41027685 0.70321205 0.75460458
> getprob1(34, 1.0, 0.777, 0.5, 0.5, 10, seed = 1)
 [1] 0.52343660 0.83195322 0.68004610 0.19952573 0.33348669 0.09619295 0.27260786 0.62125872 0.28953613 0.71272206
> ## the results will be different even if I set.seed


> ## the below will be the result I want
> set.seed(1)
> get_prob_within(34, 1.0, 0.777, 0.5, 0.5, 10)
 [1] 0.5273274 0.4073686 0.9364763 0.8362440 0.1895381 0.7563707 0.2480652 0.2958896 0.4755769 0.3697152
> get_prob_within(34, 1.0, 0.777, 0.5, 0.5, 10)
 [1] 0.44855057 0.17260087 0.17335768 0.86675369 0.32480827 0.08470480 0.10164869 0.45593050 0.80664073 0.01580556
> set.seed(1)
> get_prob_within(34, 1.0, 0.777, 0.5, 0.5, 10)
 [1] 0.5273274 0.4073686 0.9364763 0.8362440 0.1895381 0.7563707 0.2480652 0.2958896 0.4755769 0.3697152
> get_prob_within(34, 1.0, 0.777, 0.5, 0.5, 10)
 [1] 0.44855057 0.17260087 0.17335768 0.86675369 0.32480827 0.08470480 0.10164869 0.45593050 0.80664073 0.01580556

我想知道我是否可以做到这一点，我知道 Fortran 和 R 使用不同的 RNG，但我真的需要找到一种方法来组合 fortran 和 R。

谢谢你们的帮助！！ [1]：https://jblevins.org/mirror/amiller/random.f90 [2]：https://people.sc.fsu.edu/~jburkardt/f_src/cdflib/cdflib.html