From: "Mauricio Fernández" <batsman.geo@...>
Date: 2003-04-21T00:14:23+09:00
Subject: Re: Biased weighted random?

On Sun, Apr 20, 2003 at 11:35:23PM +0900, Chris Pine wrote:
> ----- Original Message -----
> From: "Mauricio Fern�ndez" <batsman.geo@yahoo.com>
> 
> I have been thinking about this for a while and even though I understand
> the subjacent idea, I am not sure this is guaranteed to work at all.
> ----------------------------
> 
> Oh, I gave no guarantees!  That's why I ended my post with "Should work much
> of the time...".
> 
> What bugs me is that I can't even tell *when* it won't work.
> 
> 
> ----------------------------
> the latter, but it is obvious that this way (weights only) we're giving
> away several degrees of freedom that could come handy later :)
> ----------------------------
> 
> Absolutely... but then we are getting close to the question of *which*
> (among several) solutions is the most psychologically satisfying.  Earlier
> in this thread, I showed that their may be more than one transition matrix
> solution.  At least in that case, the "best" solution (I think I called it
> the most "locally random" solution) was the one which corresponded to using
> modified weights.  I couldn't figure out any other way to satisfy the
> psychological requirement.

I would feel psychologically satisfied if the entropy rate was maximum.
[I guess I am rather hard to satisfy ;-) ]

Anyway, here goes one piece of code that implements a naive version of
the algorithm I described in the past. If I get some time I might as
well go for the "best" one (ie., that using a linear system in the last
step); then the "best solution" will appear by imposing constraints on
the family of solutions found.

My code uses narray, which makes it somewhat faster (one order of
magnitude) than yours, but it is in fact much heavier as I'm working w/
the whole transition matrix instead of only the weights. The nice thing
is that it gives sometimes different (but valid!) results.

WARNING: this is kinda big and ugly

require 'narray'

class NVector
  alias_method :__old_inspect__, :inspect
  def inspect
    __old_inspect__.gsub!(/\n/, " ")
  end

  def to_na
    NArray[ self.to_a ]
  end

  def idx_of_max
    idx = 0
    max = 0
    size.times { |i|
      idx, max = i, self[i].abs if self[i].abs > max
    }
    return nil if max == 0
    idx
  end
end

class NArray
  def to_vector
    ret = NVector.float(size)
    idx = 0
    each { |i| ret[idx] = i; idx += 1 }
    ret
  end
end

class NMatrix
  def col(n)
    rows = self.shape()[1]
    c = self[n, 0..(rows - 1)]
    c.to_vector
  end

  def row(n)
    cols = self.shape()[0]
    c = self[0..(cols - 1), n]
    c.to_vector
  end

  def rows
    self.shape()[0]
  end

  def cols
    self.shape()[1]
  end
end

def transition_matrix(weights)
  n = weights.size
  cols = NMatrix.float(n,n)
  for x in 0...n
    for y in 0...n
      cols[x,y] = 0
      next if x == y
      cols[x,y] = weights[x]/(1-weights[y])
    end
  end
  cols
end

def stationary_distribution matrix
  n = matrix.shape()[0]

  10.times { matrix *= matrix }
  matrix.row(0)
end


def find_suitable_row(delta, pos, matrix, probs, exclude)
  column = matrix.col(pos)
  pcorrection = nil # possible correction
  column.collect! { |i| 1 - i } if delta > 0
  #puts "COLS: #{column.inspect}"
  pcorrection = (probs.to_na * column.to_na).to_vector
  exclude.each { |i| pcorrection[i] = 0 }
  pcorrection[pos] = 0
  #puts "Maximal correction: #{pcorrection.inspect}"
  idx = 0
  max = 0
  pcorrection.size.times { |i|
    idx, max = i, pcorrection[i] if pcorrection[i] > max
  }
  return nil if max == 0
  idx
end

def correct_matrix(mat, x, y, deltas, probs )
  row = mat.row(y)
  comp = (deltas[x] > 0)? proc{ |xx| xx < 1 } : proc { |xx| xx > 0 }
  elems = []
  row.size.times { |i| elems << i if comp.call(row[i]) }
  elems -= [y, x]
  #puts "Fixing columns #{elems.inspect} in row #{y}"
  raise "CANNOT FIX" if elems.size == 0
  # all the other diffs are of the same sign => cannot correct easily
  # moreover, it cannot be solved by changing the stat. dist. because
  # even in that case the changes must be of opposite sign
  #sumdeltas = 0
  #elems.each { |i| sumdeltas += deltas[i].abs }
  elems.each { |i|
    mat[i,y] -= deltas[y]/elems.size
    mat[i,y] = 0 if mat[i,y] <= 0
    mat[i,y] = 1 if mat[i,y] >= 1
  }
  mat[x,y] += deltas[x]
  sum = 0.0
  elems.each { |i| sum += mat[i,y] }
  elems.each { |i| mat[i,y] *= (1 - mat[x,y]).abs / sum }
end

targets = NVector[ 0.4, 0.2, 0.2, 0.1, 0.001, 0.099 ]
puts "Trying to approximate #{targets.inspect}"
weights = targets.dup
exclude = {}
probs = deltas = nil
tmatrix = transition_matrix( weights )

iter = 0
catch(:done) do
  200.times do |iter|
    probs   = stationary_distribution tmatrix
    deltas = targets - probs
    #puts "Got transition matrix: #{tmatrix.inspect}"
    #puts "Corresponding stationary dist.: #{probs.inspect}"
    #puts "Delta: #{deltas.inspect}"
    #puts "Max diff: #{deltas.max}"
    throw :done if deltas.max < .00001
    i = deltas.idx_of_max
    #puts "Trying to correct #{i}"
    n = find_suitable_row(deltas[i], i, tmatrix, probs, [])
    #puts "Proceed on #{n}"
    correct_matrix(tmatrix, i, n, deltas, probs)
  end
end
puts "Needed #{iter} iterations."

puts "Got transition matrix: #{tmatrix.inspect}"

puts "Stationary distribution probabilities"
puts probs.to_a.collect { |p| "#{p} " }.join






batsman@kodos:~/germany2/src$ ruby markov2.rb
Trying to approximate NVector.float(6):  [ 0.4, 0.2, 0.2, 0.1, 0.001,
0.099 ]
Needed 112 iterations.
Got transition matrix: NMatrix.float(6,6):
[ [ 0.0, 0.339637, 0.339724, 0.161189, 0.0, 0.15945 ],
  [ 0.749724, 0.0, 0.122931, 0.0633426, 0.00125589, 0.0627467 ],
  [ 0.749559, 0.12309, 0.0, 0.0633408, 0.00126638, 0.0627433 ],
  [ 0.503809, 0.196678, 0.196678, 0.0, 0.00252385, 0.100312 ],
  [ 0.4004, 0.2002, 0.2002, 0.1001, 0.0, 0.0990991 ],
  [ 0.498523, 0.198466, 0.198466, 0.101994, 0.00255094, 0.0 ] ]
Stationary distribution probabilities
0.3999909783 0.199990125 0.1999928201 0.1000090866 0.001009401954
0.09900758802

Your program would give 
0.5888367575
0.1409106736
0.1409106736
0.0647154815
0.0006056404603
0.0640207732

which corresponds to the matrix
NMatrix.float(6,6):
[ [ 0.0, 0.342712, 0.342712, 0.157396, 0.00147299, 0.155706 ],
  [ 0.68542, 0.0, 0.164023, 0.0753303, 0.00070498, 0.0745217 ],
  [ 0.68542, 0.164023, 0.0, 0.0753303, 0.00070498, 0.0745217 ],
  [ 0.62958, 0.150661, 0.150661, 0.0, 0.000647547, 0.0684506 ],
  [ 0.589194, 0.140996, 0.140996, 0.0647547, 0.0, 0.0640596 ],
  [ 0.629113, 0.150549, 0.150549, 0.069142, 0.000647066, 0.0 ] ]

somewhat different. as you can see.


 
> In any case, as I mentioned in #69742, I think this is simply the wrong (if
> more fun) approach.
  ========

That is the most important factor for my psychological satisfaction :-)
 

-- 
 _           _                             
| |__   __ _| |_ ___ _ __ ___   __ _ _ __  
| '_ \ / _` | __/ __| '_ ` _ \ / _` | '_ \ 
| |_) | (_| | |_\__ \ | | | | | (_| | | | |
|_.__/ \__,_|\__|___/_| |_| |_|\__,_|_| |_|
	Running Debian GNU/Linux Sid (unstable)
batsman dot geo at yahoo dot com

No, that's wrong too.  Now there's a race condition between the rm and
the mv.  Hmm, I need more coffee.
	-- Guy Maor on Debian Bug#25228