| Nestor Grion 2005-11-19, 9:58 pm |
| Hello
In another post, V. Gantovnik wrote:
"How to get random number from the range 0 <=3D x <=3D 1, including 1?"
I was thinking and I believe that I found a solution. Because it is
more general I write a new post.
Problem description:
- The generator can produce N objects different (number, string...),
with equal probability (N can be unknown)
- The equal operator is defined.
- The value of an object, x0, is known (0.0, for example)
- xn is the value to add (1,0, for example)
- Each possible object has a probability of 1/(N+1)
Solution:
1) Simple Algorithm.
In order to obtain a random object (xi), two random objects are
generated (a and b). The algorithm fails 1/N^2 of the times. In this
case, two new objetcts are generated, until achieving a success.
Generate a and b.
if (a=3D=3Db) then
if(a=3D=3Dx0) then
Fail
else
xi =3D xn
endif
else
xi =3D a
endif
Then:
Pr(fail) =3D 1/N^2
Pr(xi=3Dxn) =3D (1/N)*((N-1)/N) =3D (N-1)/N^2
Pr(xi=3Dx*) =3D ((N-1)/N)/N =3D (N-1)/N^2; (x*: xi, xi /=
=3D
xn)
Considering only the valid attempts, each value has equal probability!.
2) Algorithm that makes the failure improbable.
Are desired M random values (vector X). Two vectors, A and B, are
generated (of length M).
When a =3D b, its value it is not used (xi =3D fail or xn)
Those values are put in vectors C and D (with uniform distribution!).
Repeat the simple algorithm using C and D. Note that, in average, it is
necessary to replace M/N^2 values using M/N values. Next, the failure
probability is almost null (for big M).
3) Efficient Algorithm
To obtain M values from M values (only A is generated).
Conparison:
A(i)=3D=3DA(mod(2*i, M+1))
Use algorithm 2) to solve failures.
The generated numbers have uniform distribution (no correlation between
the elements, I believe). Although some sequences it does not have the
same probability. This can be important in some cases. But I am not an
expert in these things.
In the attached program, I implement and verify the algorithm 3). The
generator produce the values 1, 2 and 3. The value 4 is added.
N=E9stor C. Grion
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
3D=3D=3D=3D
! La subrutina random_add implementa un algoritmo para agregar
! un nuevo valor a los que produce el generador de numeros
! aleatorios.
!
! El ejemplo es con enteros. Pero el algoritmo se puede aplicar
! sobre cualquier tipo de datos, mientras est=E9 definida la
! operaci=F3n de igualdad.
!
! Autor: N=E9stor C. Grion
! Fecha: 19 Nov 2005
module randext
implicit none
contains
! x: vector con random number (cantidad par de elementos)
! know: valor conocido
! new: valor a agregar
! satus: 0 si ok
subroutine random_add(x,know,new,status)
implicit none
integer :: x(:),know,new,status
integer, allocatable :: xaux(:) ! type(xaux)=3Dtype(x)
integer :: p(size(x)) ! pointer
integer :: n,i,ip,j1,j2,fail
status =3D 1
n =3D size(x,1)
! Pointer to fail and new
ip =3D 1
fail =3D 0
do i=3D1,n
if(x(i)=3D=3Dx(mod(2*i,n+1))) then
if(x(i)=3D=3Dknow) then
fail =3D fail+1
p(ip) =3D -i ! x(i) =3D ?
else
p(ip) =3D i ! x(i) =3D new
endif
ip =3D ip+1
endif
enddo
ip =3D ip-1
! Resolver no validos y escribir new
allocate(xaux(ip))
xaux =3D x(abs(p(1:ip)))
if(fail>0) then
if(ip=3D=3D1) return ! Fall=F3 el algoritmo
j1 =3D 0
j2 =3D 1
do i=3D1,ip
if(p(i)<0) then ! x(i) =3D ?
in: do
j1 =3D j1+1
j2 =3D mod(2*j1,ip+1)
if(j1>ip) return ! Fallo el algoritmo
if(xaux(j1)=3D=3Dxaux(j2)) then
if(p(j1)>0) then
x(abs(p(i))) =3D new ! x(i) =3D new
exit in
endif
else
x(abs(p(i))) =3D xaux(j1) ! x(i) =3D no new
exit in
endif
enddo in
else
x(p(i)) =3D new
endif
enddo
endif
status =3D 0
end subroutine
! vector de frecuencias
subroutine vfrec(y,vf)
implicit none
integer i,y(:)
real vf(:)
vf =3D 0.0
do i=3D1,size(y,1)
vf(y(i)) =3D vf(y(i))+1
enddo
vf =3D vf/size(y,1)
end subroutine
! Transition matrix
subroutine mcorrc(y,mc,k)
implicit none
integer i,y(:),k,n,m
real mc(:,:)
n =3D size(y,1)
m =3D size(mc,1)
mc =3D 0.0
do i=3D1,n-k
mc(y(i),y(i+k)) =3D mc(y(i),y(i+k))+1
enddo
do i =3D1,k
mc(y(n-k+i),y(i)) =3D mc(y(n-k+i),y(i))+1
enddo
mc =3D mc/spread(sum(mc,2)/m,2,m)
end subroutine
subroutine printm(x)
real x(:,:)
integer i
do i=3D1,size(x,1)
write(*,*) x(i,:)
enddo
end subroutine
end module
! Ejemplo con enteros
program main
use randext
implicit none
integer, parameter :: n =3D 2000000, mv =3D 3
real :: x(n),vf(mv+1),mc1(mv+1,mv+1),mc2(mv+1,mv
+1)
integer :: know, new
integer :: y(n),i,j,status
know =3D 1
new =3D mv+1
call random_number(x)
y =3D ceiling(mv*x)
call random_add(y,know,new,status)
write(*,*) 'status =3D',status
call vfrec(y,vf)
write(*,*) ' Frequency'
write(*,*) vf
do j=3D1,4
call mcorrc(y,mc1,j)
write(*,*) 'Transition matrix, k =3D',j
call printm(mc1)=20
enddo
end program
|