Algorithme en fortran 90 pour créer un échantillon aléatoire avec poids

Dans cet article on va voir comment mettre en oeuvre simplement un échantillonnage aléatoire et pondéré en se basant sur le papier de recherche suivant: Weighted random sampling with a reservoir. Pour cela, considérons le problème simple suivant: on dispose de 100 données avec des poids différents (ici, un poids de 1 pour les données de 1 à 79 et un poids de 4 pour les données de 80 à 100) et on souhaite par exemple prélever un échantillon aléatoirement de taille 10 (les données avec un poids plus important doivent être associées avec une probabilité plus grande d'être sélectionnées). Pour cela, le papier Weighted random sampling with a reservoir propose un algorithme avec réservoir qui nécessite de parcourir seulement une fois l'ensemble des données dont la logique est présentée sur la figure ci-dessous (pour une démonstration détaillée veuillez consulter l'article). Voici un exemple de code en fortran 90 basé sur cette logique:

Weighted random sampling with a reservoir)">

program WeightRandomSamplingReservoir

implicit none

integer :: i,j
integer*4 timeArray(3) ! Holds the hour, minute, and second

integer, dimension(100,2) :: V ! Population
real, dimension(10,2) :: R ! Reservoir
integer, dimension(1) :: MinIndex
integer, dimension(2) :: R_dim

real :: ui,ki
real :: ki_min

call itime(timeArray) ! Generate sequence of random number

i = rand ( timeArray(1)+timeArray(2)+timeArray(3) )
!do i = 1, 10
!print *, rand(0)
!end do

R_dim = shape(R)

!----- Population V: (vi,wi) -----!

V = 0 
do i = 1, 100
  V(i,1) = i
  if( i < 80 )then
    V(i,2) = 1
  else
    V(i,2) = 4
  end if
  !write(6,*) 'Data, Weight: ',V(i,1),V(i,2)
end do

!----- Weighted Random Sampling -----!

ui = rand(0)
R(1,1) = V(1,1)
R(1,2) = ui**(1.0/V(1,2))
ki_min = R(1,2)
do i = 2, R_dim(1)
  ui = rand(0)
  R(i,1) = V(i,1)
  R(i,2) = ui**(1.0/V(i,2))
  if( R(i,2) < ki_min ) ki_min = R(i,2)
end do


do i = R_dim(1)+1, 100
  ui = rand(0)
  ki = ui**(1.0/V(i,2))

  !write(6,*) V(i,1),V(i,2),ki
  !write(6,*) (R(j,1), j=1,10)
  !write(6,*) (R(j,2), j=1,10)
  !write(6,*) minloc(R(:,2)),'dimension',shape(minloc(R)),MinIndex(1)
  !write(6,*) "------------"

  MinIndex = minloc(R(:,2))

  if( ki > ki_min )then
    j = MinIndex(1)
    R(j,1) = V(i,1)
    R(j,2) = ki

    MinIndex = minloc(R(:,2))
    j = MinIndex(1)
    ki_min = R(j,2)
  end if

  !if( rand(0) < 0.2 )then
  !write(6,*) D(i,1),D(i,2)
  !end if
end do

!----- Print Random Sampling -----!

do i=1,R_dim(1)
write(6,*) R(i,1), R(i,2)
end do

end program WeightRandomSamplingReservoir

Pour compiler le code: "fortran WeightRandomSamplingReservoir.f90 -o WeightRandomSamplingReservoir". Le code ci dessus donne par exemple l'échantillon de taille 10 suivant (colonne 1: donnée sélectionnée; colonne 2: clé ki associée):

59.0000000      0.966274738    
2.00000000      0.959636211    
92.0000000      0.940751970    
39.0000000      0.957770824    
5.00000000      0.954653263    
89.0000000      0.963813066    
96.0000000      0.942957997    
45.0000000      0.960242271    
80.0000000      0.968827009    
88.0000000      0.988696277

Recherches associées