Fortran 中的加权采样
Weighted sampling in Fortran
在 Fortran 程序中,我想通过使用权重随机选择一个特定变量(特别是它的索引)。权重将在单独的向量中提供(元素 1 将包含变量 1 的权重,依此类推)。
我有以下代码在没有权重的情况下完成工作(mind
是一个整数向量,具有原始数据集中每个变量的索引)
call rrand(xrand)
j = int(nn * xrand) + 1
mvar = mind(j)
这里有两个例子。第一个是
integer, parameter :: nn = 5
real :: weight( nn ), cumsum( nn ), x
weight( 1:nn ) = [ 1.0, 2.0, 5.0, 0.0, 2.0 ]
do j = 1, nn
cumsum( j ) = sum( weight( 1:j ) ) / sum( weight( 1:nn ) ) !! cumulative sum
enddo
x = rand()
do j = 1, nn
if ( x < cumsum( j ) ) exit
enddo
第二个取自this page
real :: sum_weight
sum_weight = sum( weight( 1:nn ) )
x = rand() * sum_weight
do j = 1, nn
if ( x < weight( j ) ) exit
x = x - weight( j )
enddo
这与第一个基本相同。两者都从 1,2,...,5 中随机抽取一个 j
,权重为 (j)。 100000 次试验给出的分布类似于
j : 1 2 3 4 5
count : 10047 19879 50061 0 20013
编辑:下面附有最小测试代码(使用 gfortran-8/9 测试):
program main
implicit none
integer j, num( 5 ), loop
real weights( 5 )
weights(:) = [ 1.0, 2.0, 5.0, 0.0, 2.0 ]
num(:) = 0
do loop = 1, 100000
call random_index( j, weights )
num( j ) = num( j ) + 1
enddo
do j = 1, size( weights )
print *, j, num( j )
enddo
contains
subroutine random_index( idx, weights )
integer :: idx
real, intent(in) :: weights(:)
real x, wsum, prob
wsum = sum( weights )
call random_number( x )
prob = 0
do idx = 1, size( weights )
prob = prob + weights( idx ) / wsum !! 0 < prob < 1
if ( x <= prob ) exit
enddo
end subroutine
end program
在 Fortran 程序中,我想通过使用权重随机选择一个特定变量(特别是它的索引)。权重将在单独的向量中提供(元素 1 将包含变量 1 的权重,依此类推)。
我有以下代码在没有权重的情况下完成工作(mind
是一个整数向量,具有原始数据集中每个变量的索引)
call rrand(xrand)
j = int(nn * xrand) + 1
mvar = mind(j)
这里有两个例子。第一个是
integer, parameter :: nn = 5
real :: weight( nn ), cumsum( nn ), x
weight( 1:nn ) = [ 1.0, 2.0, 5.0, 0.0, 2.0 ]
do j = 1, nn
cumsum( j ) = sum( weight( 1:j ) ) / sum( weight( 1:nn ) ) !! cumulative sum
enddo
x = rand()
do j = 1, nn
if ( x < cumsum( j ) ) exit
enddo
第二个取自this page
real :: sum_weight
sum_weight = sum( weight( 1:nn ) )
x = rand() * sum_weight
do j = 1, nn
if ( x < weight( j ) ) exit
x = x - weight( j )
enddo
这与第一个基本相同。两者都从 1,2,...,5 中随机抽取一个 j
,权重为 (j)。 100000 次试验给出的分布类似于
j : 1 2 3 4 5
count : 10047 19879 50061 0 20013
编辑:下面附有最小测试代码(使用 gfortran-8/9 测试):
program main
implicit none
integer j, num( 5 ), loop
real weights( 5 )
weights(:) = [ 1.0, 2.0, 5.0, 0.0, 2.0 ]
num(:) = 0
do loop = 1, 100000
call random_index( j, weights )
num( j ) = num( j ) + 1
enddo
do j = 1, size( weights )
print *, j, num( j )
enddo
contains
subroutine random_index( idx, weights )
integer :: idx
real, intent(in) :: weights(:)
real x, wsum, prob
wsum = sum( weights )
call random_number( x )
prob = 0
do idx = 1, size( weights )
prob = prob + weights( idx ) / wsum !! 0 < prob < 1
if ( x <= prob ) exit
enddo
end subroutine
end program