Взвешенная выборка в Фортране

В программе на языке Fortran я хотел бы случайным образом выбрать конкретную переменную (в частности, ее индекс) с помощью весов. Веса будут представлены в отдельном векторе (элемент 1 будет содержать вес переменной 1 и так далее).

У меня есть следующий код, который выполняет работу без веса (mind является целочисленным вектором с индексом каждой переменной в исходном наборе данных)

call rrand(xrand)
j = int(nn * xrand) + 1
mvar = mind(j)

fortran weighted sampling
person Chris    schedule 02.09.2015    source источник
comment
Вы можете просто добавить петлю   -  person innoSPG    schedule 02.09.2015

Ответы (1)


arrow_upward
4
arrow_downward

Вот два примера. Первый

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

а второй взят с этой страницы

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

который по сути совпадает с первым. Оба выбирают случайный j из 1,2, ..., 5 с весом (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
person roygvib    schedule 02.09.2015
comment
спасибо ройдвиб! Я предполагаю, что j после выхода из цикла сохраняет индекс выбранной переменной. Правильно? - person Chris; 03.09.2015
comment
Да, значение j сохраняется после выхода из цикла. Таким образом, мы можем написать функцию, которая получает вес (:) в качестве фиктивного аргумента и возвращает j, полученный, например, из цикла DO. - person roygvib; 03.09.2015
comment
@roygvib спасибо, но я не уверен, что понимаю, как это работает. Кажется, я всегда получаю одно и то же j, когда запускаю это. Не могли бы вы рассказать подробнее? Кажется, каждый раз, когда я запускаю эту функцию, rand () дает одно и то же число. - person Herman Toothrot; 07.01.2020
comment
@HermanToothrot Привет, я прикрепил минимальный тестовый / рабочий код внизу своего ответа. Код немного изменен, чтобы подпрограмму random_index() можно было использовать автономно. Не могли бы вы попробовать этот код в своей среде, чтобы убедиться, что он работает должным образом? (Я почти забыл детали, но я предполагаю, что мы вычисляем кумулятивную вероятность из заданных весов и проверяем, где 'x' попадает в какую-то ячейку (с заданным весом).) - person roygvib; 08.01.2020
comment
@roygvib спасибо, да работает, я не понял что должно быть в шлейфе. Я фактически реализовал это с помощью второго примера, который вы предоставили. - person Herman Toothrot; 09.01.2020
comment
@roygvib, возможно, я нашел ошибку в этом коде, но хочу уточнить у вас. В случае, когда x == weight (j), ваш второй пример терпит неудачу и возвращает индекс, который находится вне цикла. Итак, «if (x‹ weight (j)) exit »в вашем примере должно быть« if (x ‹= weight (j)) exit», что, я думаю, вы используете в своем тестовом коде. - person Herman Toothrot; 09.09.2020
comment
@HermanToothrot Кажется, что первый и второй фрагменты кода (которые, я думаю, основаны на предыдущей странице вопросов и ответов) используют ‹, а третий фрагмент использует‹ =. Думаю, я боялся углового случая x, близкого к 1. Хотя в идеале x гарантированно меньше 1, x == 1.0 или даже x ›1.0 (немного выше) может произойти из-за ошибки округления (я предполагаю ). Итак, я думаю, что может быть лучше провести сравнение для idx ‹= nn -1 (которое может использовать‹ или ‹=), и, если неравенство не соблюдается, вернуть в качестве результата nn (т. Е. Мы рассматриваем x находиться в последнем интервале весов) ... - person roygvib; 14.09.2020
comment
@roygvib Я думал, что проблема была решена простым использованием x ‹= weight (j), но сбой кода вернул меня к реальности :). Итак, я думаю, вы правы, что цикл должен идти до nn-1, поэтому idx может быть не более nn. Я не уверен, что вы имеете в виду с idx ‹= nn -1, я не могу скомпилировать с do idx‹ = 1, size (weights) -1, поэтому я просто изменил, чтобы do idx = 1, size (weights) -1 . Следует ли нам изменить это в вашем ответе? Также я использую пример 2, поэтому мне интересно, может ли это случиться и в примере 1. - person Herman Toothrot; 05.10.2020
comment
@HermanToothrot, извините за поздний ответ ... (Я собирался подумать над кодом более внимательно в выходные, но не смог этого сделать, потому что в последнее время меня очень легко забывают ...) И да, если возможно более надежный код, я предлагаю вам добавить еще один ответ отдельно (в качестве улучшения). Кроме того, для дальнейшего обсуждения я полагаю, что один из reddit или Fortran Discourse или comp.lang.fortran будут очень полезны для получения дополнительных отзывов / мнений :) - person roygvib; 05.11.2020