Re: Help with generating a random matrix with row and column constraint

"George Marsaglia" <geo@xxxxxxxxxxxx> wrote in message
news:dcqdnXn3Odc9_DXeRVn-tw@xxxxxxxxxxxxxx

> you might generate 1500 independent
> truncated Poisson variates with
> parameter 1/3, so that the
> expected col sums will be 100.
> Then about 4% of the time, the sum
> will be exactly 100.

I tried a couple of the methods suggested in this thread:

! File: ran_mat.f90
! Public domain 2005 James Van Buskirk

module col_funcs
implicit none
contains
function col_marsaglia(len_in,sum_in,max_in)
integer, intent(in) :: len_in
integer, intent(in) :: sum_in
integer, intent(in) :: max_in
integer col_marsaglia(len_in)
real probs(0:max_in)
integer i
real lambda
integer j
real temp

if(len_in*max_in < sum_in .OR. sum_in < 0) then
write(*,'(a)') 'Error in col_marsaglia: impossible sum'
stop
end if
lambda = real(sum_in)/len_in
probs(0) = 1
do i = 1, max_in
probs(i) = probs(i-1)*lambda/i
end do
probs = probs/sum(probs)
do i = 1, max_in
probs(i) = probs(i-1)+probs(i)
end do
probs(max_in) = 1
do
do j = 1, len_in
call random_number(temp)
do i = 0, max_in
if(temp <= probs(i)) exit
end do
col_marsaglia(j) = i
end do
if(sum(col_marsaglia) == sum_in) exit
end do
end function col_marsaglia

function col_rusin(len_in,sum_in,max_in)
integer, intent(in) :: len_in
integer, intent(in) :: sum_in
integer, intent(in) :: max_in
integer col_rusin(len_in)
type bin
integer position
integer offset
end type bin
type(bin) bins(len_in*(max_in+1))
integer j
integer i
real temp

if(len_in*max_in < sum_in .OR. sum_in < 0) then
write(*,'(a)') 'Error in col_rusin: impossible sum'
stop
end if
bins%position = ((/(i,i=1,size(bins))/)-1)/(max_in+1)+1
bins%offset = mod((/(i,i=1,size(bins))/)-1,max_in+1)
do j = size(bins), size(bins)-sum_in+1, -1
call random_number(temp)
i = j*temp+1
if(i < j) bins((/i,j/)) = bins((/j,i/))
end do
col_rusin = 0
do j = j+1, size(bins)
col_rusin(bins(j)%position) = col_rusin(bins(j)%position)+1
end do
end function col_rusin

subroutine print_results(title, ran_mat, sum_in, max_ent, max_row)
character(*), intent(in) :: title
integer, intent(in) :: ran_mat(:,:)
integer, intent(in) :: sum_in
integer, intent(in) :: max_ent
integer, intent(in) :: max_row
integer i

write(*,'(a)') 'Name of procedure: '//title
write(*,'(a,i0,"X",i0)') 'Shape of matrix = ', &
shape(ran_mat)
write(*,'(a,i0)') 'Number of bad entries = ', &
count(ran_mat > max_ent)
write(*,'(a,i0)') 'Number of bad columns = ', &
count(sum(ran_mat,1) /= sum_in)
write(*,'(a,i0)') 'Number of bad rows = ', &
count(sum(ran_mat,2) > max_row)
write(*,'(a)') 'Value Frequency'
do i = 0, max_ent
write(*,'(i5,i10)') i, count(ran_mat == i)
end do
end subroutine print_results
end module col_funcs

program test
use col_funcs
implicit none
integer, parameter :: max_ent = 5
integer, parameter :: target = 100
integer, parameter :: max_row = 5
integer ran_mat(300,5)
integer j

call random_seed()
do
do j = 1, size(ran_mat,2)
ran_mat(:,j) = col_marsaglia(size(ran_mat,1), target, max_ent)
end do
if(all(sum(ran_mat,2) <= max_row)) exit
end do
call print_results('Marsaglia',ran_mat,target,max_ent,max_row)
do
do j = 1, size(ran_mat,2)
ran_mat(:,j) = col_rusin(size(ran_mat,1), target, max_ent)
end do
if(all(sum(ran_mat,2) <= max_row)) exit
end do
write(*,'()')
call print_results('Rusin',ran_mat,target,max_ent,max_row)
end program test

Sample results:

Name of procedure: Marsaglia
Shape of matrix = 300X5
Number of bad entries = 0
Number of bad columns = 0
Number of bad rows = 0
Value Frequency
0 1087
1 337
2 66
3 9
4 1
5 0

Name of procedure: Rusin
Shape of matrix = 300X5
Number of bad entries = 0
Number of bad columns = 0
Number of bad rows = 0
Value Frequency
0 1067
1 371
2 57
3 5
4 0
5 0

It seems that these methods don't generate very many 5's,
though.

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end


.

Relevant Pages

  • Re: has this been done?
    ... of entries in the matrix subject to row sum and column ... 2x2 matrix with row and column sums A and B> 0. ... nonnegative integer values. ... average entry is 1, so that an entry ...
    (sci.math)
  • Re: How to sum the last 7 number in a row?
    ... search this discussion board but cant find anything. ... Anyways I was wondering if there is a formula that sums the last 7 entries I ... number at the back and do not want to change the range everytime I do so. ...
    (microsoft.public.excel.misc)