./math/rands.rb

download original
# -*- ruby -*-

# some random number generators


# equally distributed random numbers between 0 and 1
$standardRand = proc { rand }

# gaussian random numbers, mean 0, std. deviation 1
# uses "Box-Muller transformation", see http://www.taygeta.com/random/gaussian.html
$gaussRand =
  proc {
    begin
      x1 = 2.0 * rand - 1.0;
      x2 = 2.0 * rand - 1.0;
      w = x1 * x1 + x2 * x2;
    end until ( w < 1.0 );
  w = Math.sqrt( (-2.0 * Math.log( w ) ) / w );
  x1 * w;
}

# return a generator that "cuts" the output numbers of generator rand
# to fit into (lowerBound, upperBound)
def boundedRand(rand, lowerBound, upperBound)
  proc {
      begin
	r=rand.call()
      end until (r>lowerBound) and (r<upperBound)
    r
  }
end


# return a generator that adds shift to the
# output numbers of generator rand
def shiftedRand(rand, shift)
  proc {
      rand.call() + shift
  }
end


# return a generator that multiplies the
# output numbers of generator rand by scale
def scaledRand(rand, scale)
  proc {
      rand.call() * scale
  }
end

# unify two generators rand1 and rand2, with "preference"
# of rand1 given by rand1amount (0..1)
def unionRand(rand1, rand2, rand1amount)
  proc {
      r = rand
      if r<rand1amount
      then rand1.call
      else rand2.call
      end
  }
end

# add two generators rand1 and rand2
def sumRand(rand1, rand2)
  proc {
      rand1.call + rand2.call
  }
end


# generate and return  histogram of generator rand
# output format [[x1,y1],[x2,y2],...,[xn,yn]]
def calcRandHistogram2(rand, lowerBound, upperBound, nSections=100, nRuns=10000)
  bRand = boundedRand(rand, lowerBound, upperBound)
  sectionWidth = (upperBound - lowerBound) / nSections

  result = Array.new(nSections) {|i| [lowerBound + i*sectionWidth, 0]}

  nRuns.times { |i|
    r = bRand.call
    result[(r-lowerBound)/sectionWidth][1] += 1
  }

  result
end

def getRandMeanAndStdDev(rand, nRuns=1000)
  sum = 0.0; sqSum = 0.0
  nRuns.times {
    r = rand.call
    sum += r
    sqSum += r * r
  }
  mean = sum / nRuns
  stdDev = Math.sqrt(sqSum / nRuns - mean * mean)
  [mean, stdDev]
end

def getRandMean(rand, nRuns=1000)
  getRandMeanAndStdDev(rand,nRuns)[0]
end

def getRandStdDev(rand, nRuns=1000)
  getRandMeanAndStdDev(rand,nRuns)[1]
end

def getRandVariance(rand, nRuns=1000)
  stddev = getRandMeanAndStdDev(rand,nRuns)[1]
  stddev * stddev
end

# generate and return  histogram of generator rand
# output format [[x1,x2,...,xn],[y1,y2,...yn]] (suitable for plotRandHistogram)
def calcRandHistogram(rand, lowerBound, upperBound, nSections=100, nRuns=10000)
  bRand = boundedRand(rand, lowerBound, upperBound)
  sectionWidth = (upperBound - lowerBound) / nSections

  xs = (0...nSections).collect {|i| lowerBound + i*sectionWidth}
  ys = [0]*nSections

  nRuns.times { |i|
    r = bRand.call
    ys[(r-lowerBound)/sectionWidth] += 1
  }

  [xs,ys]
end


# plot histgram as generated by calcRandHistogram
# requires Ruby Gnuplot
def plotRandHistogram(hist,out=nil)
  require 'gplot/Gnuplot'
  include Gnuplot

  xs,ys = hist

  plot = Plot.new(out)
  plot.title("random number generator histogram")
  plot.xlabel("domain")
  plot.ylabel("amount")

  plot.draw(ArrayDataSet.new(ys, "with"=>"point 3 3", "title"=>nil, "xgrid"=>xs))
end


if __FILE__ == $0
  myrand=unionRand(shiftedRand(scaledRand($standardRand,5.0),
                               -2.0),
                   shiftedRand(scaledRand($gaussRand,0.7),
                               1.0),
                   0.5)
  plotRandHistogram(calcRandHistogram(myrand,-4.0,4.0), (ARGV[0] or nil))
end


############################

$quadRand =
  proc {
  (2.0 * rand - 1.0) * (2.0 * rand - 1.0)
}


def calcExactQuadRandHistogram(lowerBound, upperBound, nSections=100, nRuns=10000)
  def P(x,dx)
    0.5*(dx-Math.log(x+dx)*(x+dx)+x*Math.log(x))
  end

  sectionWidth = (upperBound - lowerBound) / nSections

  xs = (0...nSections).collect {|i| lowerBound + i*sectionWidth}
  ys = (0...nSections).collect {|i|
    x = lowerBound + i*sectionWidth
    if x.abs>1.0
      0.0
    else
      if x.abs<1e-7
        nRuns * P(1e-4,sectionWidth)
      else
        if x>0.0
          nRuns * P(x,sectionWidth)
        else
          nRuns * P(-x,sectionWidth)
        end
      end
    end
  }

  [xs,ys]
end
back to math