Nonuniform Random Number Generators IV

Io has a normally distributed random number generator:

bash $ io
Io 20090105
Io> Random gaussian
==> 1.0989883665373978
Io> Random gaussian(100,10)
==> 88.2743995098800269

The underlying source code uses the Box-Muller transform:

double RandomGen_gaussian(RandomGen *self, double m, double s)
{
		// http://www.taygeta.com/random/gaussian.html
	double x1, x2, w, y1, y2;

	do {
		x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
		x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
		w = x1 * x1 + x2 * x2;
	} while ( w >= 1.0 );

	w = sqrt( (-2.0 * log( w ) ) / w );
	y1 = x1 * w;
	y2 = x2 * w;

		// The following code resulted in a lot of nans being returned. The
		// following code *should* also be slower.
		/*
	double x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
	double x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
	double y1 = sqrt( - 2.0 * log(x1) ) * cos( M_PI_2 * x2 );
		*/

	return ( m + y1 * s );
}

Factor also has a normally distributed generator:

bash $ factor
( scratchpad ) USE: random
( scratchpad ) 0 1 normal-random-float

--- Data stack:
-1.651647700885478

The interesting story here is the Factor source is algorithmically identical to what is commented out in the Io source:

: normal-random-float ( mean sigma -- n )
    0.0 1.0 uniform-random-float
    0.0 1.0 uniform-random-float
    [ 2 pi * * cos ]
    [ 1.0 swap - log -2.0 * sqrt ]
    bi* * * + ;

Nonuniform Random Number Generators III

So a normally distributed random variable is easy to simulate, but what about the other distributions? It occurred to me that the inverse of the cumulative distribution function (qv quantile function) of an RV will convert a standard uniform random number into a sample from the RV.

Since most languages don’t have cumulative distribution functions, much less quantile functions, this converts the problem into an exercise in numerical analysis. I wondered if this was the technique used in the R source code. It wasn’t so for either rgamma or rnorm. rnorm underlyingly calls this routine, which consists of almost 300 lines of C:

/*
 *  REFERENCE
 *
 *    Ahrens, J.H. and Dieter, U.
 *    Extensions of Forsythe's method for random sampling from
 *    the normal distribution.
 *    Math. Comput. 27, 927-937.
 *
 *    The definitions of the constants a[k], d[k], t[k] and
 *    h[k] are according to the abovementioned article
 */
double norm_rand(void)
{

The paper mentioned in the comment doesn’t appear to be available online, but it is summarized here. It goes back to an idea von Neumann had in the 1950s. The technique still involves multiple uniform distribution samples, which are apparently cheaper than evaluating a quantile function. It uses far fewer samples on average than the 30 I used in rnorm.

Here is a relevant Knuth quote:

It is conceivable that someday somebody will invent a random number generator that produces one of these other random quantities directly, instead of getting it indirectly via the uniform distribution. But no direct methods have as yet proved to be practical, except for the “random bit” generator described in Section 3.2.2.

The Art of Computer Programming: Seminumerical Algorithms, p. 119

Nonuniform Random Number Generators II

According to the central limit theorem, a large enough sample from any distribution will converge to a normal distribution. So we might as well use the built in uniform distribution, which has variance 1/12, and the fact that the sum of independent random variables has variance equal to the sum of the variances.

def rnorm(mean=0, sd=1)
  n = 30
  sum = Array.new(n) { rand }.inject { |m,o| m+o }
  zscore = (sum - (n/2.0))/(n/12.0)**0.5
  zscore * sd + mean
end

How to make a quick ASCII histogram of the results:

h = Hash.new(0)
1000.times { h[(10*rnorm).to_i.to_f/10] += 1 }
(-35..35).each { |k| puts "#{k.abs/10.0} #{'#' * h[k/10.0]}" }

A javascript version:

function rnorm(mean,sd) {
  n=30;
  if (null == mean) { mean = 0; }
  if (null == sd) { sd = 1;}
  for(sum=0, i=0;i < n; i++) {
    sum += Math.random();
  }
  zscore = (sum - (n/2))/Math.pow(n/12,0.5);
  return zscore * sd + mean;
}

6/16 update: code fix: tripped up by the fact that 3/2 == 1 in ruby. In javascript 3/2 == 1.5

Nonuniform Random Number Generators

A nice R feature that I haven’t seen in any other language is random number generators for distributions other than the uniform distribution. The most useful is the normal distribution: here’s a crude simulation of rnorm in ruby

def rnorm(mean=0, sd=1)
  coin_flips=37
  binomial = Array.new(coin_flips) { rand(2) }.inject { |m,o| m+o }
  zscore = 2 * (binomial - (coin_flips/2))/coin_flips**0.5
  zscore * sd + mean
end

and in javascript:

function rnorm(mean,sd) {
  coin_flips=37;
  if (null == mean) { mean = 0; }
  if (null == sd) { sd = 1;}
  binomial=0;
  for(i=0;i < coin_flips; i++) {
    if ( 0.5 < Math.random() ) {
      binomial += 1;
    }
  }
  zscore = 2 * (binomial - (coin_flips/2))/Math.pow(coin_flips,0.5);
  return zscore  * sd + mean;
}

Unlike the R version, there are only 38 possible values that above functions can assume, which might not be acceptable if repeated trials are being made. To fix that, here is a ruby version which mixes in another random number.

def fuzz(width)
  [(rand-0.5)*width, ((width**2)/12.0)**0.5]
end

def rnorm(mean=0, sd=1)
  coin_flips=37
  fuzz, fuzz_sd = fuzz(1/coin_flips**0.5)
  binomial = Array.new(coin_flips) { rand(2) }.inject { |m,o| m+o }
  zscore = 2 * (binomial - (coin_flips/2))/coin_flips**0.5
  zscore * (sd**2-fuzz_sd**2)**0.5 + mean + fuzz
end

Transition Matrix Limits

I don’t remember whether I had a copy of Gamma World, but I remember the Artifact Use Charts. These were employed when trying to learn how to use some artifact from before the nuclear holocaust. They are state diagrams, and you advance through the diagram by rolling a die until you master the item or it explodes.

The puzzle which has been hanging over me since my childhood is to figure out what the probability of an explosion is, given that you will roll until success or failure. I think the easiest way is to convert the state diagram to a transition matrix and using software which performs matrix multiplication compute ever higher powers of the matrix until the entries approach limit values. The values in the first row will then represent the probabilities of your ultimate destination.

I also wrote some ruby code which computes an exact solution via a non-analytic technique. The idea is to choose a transient state and to remove it, replacing the n x n matrix with an equivalent n-1 x n-1 matrix. In the new matrix, the probability of going from i to j includes the probability of going from i to t to j, where t is the transient state that was removed. The solution is arrived at by removing all the transient states except the starting state.

require 'pp'
require 'rational'

class TransitionMatrix

  def initialize(rows)
    @rows = rows
    @n = @rows.size
  end

  def remove(row)
    raise "cannot remove terminal row #{row}" if @rows[row][row] == 1
    out = (0...(@n-1)).collect { |i| (0...(@n-1)).collect { |j| 0 } }
    (0...@n).each do |r|
      next if r == row
      (0...@n).each do |c|
        rout = r < row ? r : r - 1
        cout = c < row ? c : c - 1
        if c == row
          (0...@n).each do |i|
            next if i == row
            iout = i  0 and @rows[r][r] < 1
        m = 1/(1-@rows[r][r])
        (0...@n).each do |c|
          @rows[r][c] = c == r ? 0 : @rows[r][c] * m
        end
      end
    end
  end

  def print
    @rows.each { |r| pp r }
  end

  def dup
    TransitionMatrix.new(@rows.collect { |r| r.dup })
  end

  def remove_loops!
    (0...@n).each do |r|
      if @rows[r][r] > 0 and @rows[r][r] < 1
        m = 1/(1-@rows[r][r])
        (0...@n).each do |c|
          @rows[r][c] = c == r ? 0 : @rows[r][c] * m
        end
      end
    end
  end

  def check
    @rows.each_with_index do |r,i|
      raise "bad row #{i}: #{r.join(', ')}" if r.inject { |m,o| m+o } != 1
    end
  end

end

def r(n,d)
  Rational(n,d)
end

chartA = TransitionMatrix.new(
[[r(3,10),r(7,10),0,0,0,0,0,0,0],
 [0,r(2,10),r(3,10),r(5,10),0,0,0,0,0],
 [0,0,r(8,10),r(2,10),0,0,0,0,0],
 [0,0,0,r(4,10),r(3,10),0,r(3,10),0,0],
 [0,0,0,0,1,0,0,0,0],
 [r(2,10),r(2,10),0,0,0,r(6,10),0,0,0],
 [0,0,0,0,r(1,10),r(4,10),r(3,10),r(2,10),0],
 [0,r(4,10),0,0,0,r(3,10),0,0,r(3,10)],
 [0,0,0,0,0,0,0,0,1],
])

aout = chartA.remove(7).remove(6).remove(5).remove(3).remove(2).remove(1)
aout.remove_loops!
aout.print

Selling Books on Amazon

I pulled a few hundred books off my shelves. The most sensible way to get rid of books is to throw them away. That’s what I did with the Encyclopedia Britannica after Good Will refused to accept it. I’m trying to get the other books to somebody who wants them, though.

Few used bookstores in LA buy books. A coworker pointed me to a used book store in Pasadena that does. The owner gave me $100 for half my books. The other half he didn’t want, but he was willing to recycle them.

I also listed books on Amazon. I listed anything where the lowest available used book price was above $10. Listing the books was easy, but packing and shipping the books involved learning. For packing materials, I eventually settled on the durable but unpadded $.99 mailers sold by the post office. I got 200′ of bubble wrap at Staples for $.10 a foot for padding. Later I put comic book backer boards in there as well.

The first 9 books I sent parcel post, and the cost of shipping was well above the $3.99 Amazon collects for me. I would have given up on the enterprise, except that on the 2nd day the post office clerk pointed out to me it is cheaper to send books by media mail. That got average cost to ship a book down to about $5. My books were often heavier than 2 lbs.

When the clerk changed the shipping from parcel to media, she ending up charging me both parcel and media rate for one of the packages. I first noticed this when I was entering the data into a spreadsheet. The post office refused to refund me when I brought this to their attention, though.

I’ve received one email from my buyers so far:

I received it today. Because of the cheap packaging and the fact that it was not stamped DO NOT BEND the USPS successfully folded the entire book in half to stuff it into the mailbox pretty much destroying it. I don’t want money back because I have been a retailer of books at Amazon myself. I just wanted to let you know so maybe you can improve your packaging processes.

Thanks

The book that got folded in half was a hardcover. I refunded the guy, and Amazon had the decency to refund their fee to me. Hopefully there won’t be any more bad news. I collected my numbers into a spreadsheet, and they look good.

Of the 57 books I initially listed, 36 have sold, 10 are still listed, and 11 got pulled off because my price got undercut below $10. The Amazon API makes it easy to see when somebody is undercutting you, and makes it easy to undercut them in turn. I checked my inventory every couple of days, and tried to remain the price leader.

Of the books the sold, I scatterplotted the log of the Amazon sales rank against days to sale. The correlation was 0.52, not as high as I thought. Also I graphed the linear regression line onto the scatterplot, and the slope was barely positive. If I were to hazard a guess for the best fit line, I would have drawn it through the origin.

Errata for do Carmo’s “Differential Geometry of Curves and Surfaces”

Here is a PDF of errata for Manfredo do Carmo’s Differential Geometry of Curves and Surfaces.

And here are a few additions to the errata sheet:

p.35 In the 4th equation there is a missing minus sign: x/y’ = – y/x’, and the following fraction of radicals should have a plus or minus.
p.38 The tangent indicatrix of the first example is not quite right. The curve should extend beyond the positive x-axis to the vector drawn in the first quadrant before turning around.
p.194 In 4th line: replace x = kxy with z = kxy
p.234 After Theorema Egregium, y should be consistently in bold font
p.329 In the statement of Proposition 2: parametried -> parametrized

Idiot notes:

p.140 Making use of remark 1 on p. 72 in proof of Proposition 1.
p.194 For (6), recall that EG-F² = | x_t ^ x_u|² (see p.156)

Linear Algebra Homework

There was a homework that has been on the back of my mind for several years.

It contained a problem which asks the reader to prove that a finite vector space has a well-defined dimension. I would have proved this by building a matrix from the coefficients you get by expressing x₁,…,x_m as linear combinations of y₁,…,y_n and doing row operations on the augmented matrix. But the problem instructs the reader to use a different approach:

The TA’s comments were unclear to me. I assumed I made an error and moved on to the next assignment. But reading my solution again and comparing it to the proof in Ch. 2 of Axler makes me think my proof is correct and the TA just failed to see that I specified an iterative process.

For reference:

Math Story I: Memorizing Pi

toothpastefordinner.com

When I was in 9th grade I attended a junior high, but I got up early every morning, as did my friend Lance, so we could take a geometry class at the high school. I showed off to one of the sophomores by writing out 14 digits of pi. He looked at what I had written and wrote a 4 at the end. “Really?”, I said, amazed. He was just bluffing. Still, the lesson is clear: never give out all your digits of pi.

Math Reading List

These texts give practice in solving hard problems like on the Putnam exam or at a mathematics olympiad. They only require knowledge of high school math:

Herman, Kucera, & Simsa, Equations and Inequalities
Engel, Problem-Solving Strategies
Herman, Kucera, and Simsa, Counting and Configurations

These books cover core math topics at the upper undergraduate or 1st year graduate level:

Mendelson, Number Systems and the Foundations of Analysis¹
Patty, Foundations of Topology, ch. 1-7 ²
Rudin, Principles of Mathematical Analysis, ch. 3-8 ³
Axler, Linear Algebra Done Right⁴
Edwards, Advanced Calculus of Several Variables, ch. II-VI
de Souza & Silva, Berkeley Problems in Mathematics
Gamelin, Complex Analysis, ch. 1-11
Dummit & Foote, Abstract Algebra
Folland, Real Analysis, ch. 1-6
Hatcher, Algebraic Topology
Lee, Smooth Manifolds

These are on supplementary topics at the upper undergraduate or 1st year graduate level:

Kincaid & Cheney, Numerical Analysis: The Mathematics of Scientific Computing
Folland, Fourier Analysis and Its Applications
Brauer & Nohel, The Qualitative Theory of Ordinary Differential Equations
McOwen, Partial Differential Equations
Hogg & Craig, Introduction to Mathematical Statistics ⁵
Athreya & Lahiri, Measure Theory and Probability Theory, ch. 6-15
Apostol, Introduction to Analytic Number Theory^5a
do Carmo, Differential Geometry of Curves and Surfaces
Cox, Little, & O’Shea, Ideals, Varieties, and Algorithms
Hall, Lie Groups, Lie Algebras, and Representations
Milnor, Topology from the Differentiable Viewpoint
Hartshorne, Geometry: Euclid and Beyond⁶
Silverman & Tate, Rational Points on Elliptic Curves
Moschovakis, Notes on Set Theory

These are 2nd year graduate texts:

Bump, Lie Groups
Hormander, The Analysis of Linear Partial Differential Operators I
Stein, Introduction to Fourier Analysis on Euclidean Spaces
Mac Lane, Categories for the Working Mathematician
Eisenbud, Commutative Algebra with a View Toward Algebraic Geometry
Shafarevich, Basic Algebraic Geometry 1: Varieties in Projective Space ⁷
Shafarevich, Basic Algebraic Geometry 2: Schemes and Complex Manifolds
Lax, Functional Analysis

¹This book teaches basic logic and set theory, and then constructs the real numbers via the counting numbers, integers, and rationals. Cauchy sequences are used for the construction of the reals. This seems to me like a good introduction to higher math, but the book has fallen out of print and I’m not aware of any equivalent texts.
²Another book which undertakes teaching basic mathematical reasoning and which has fallen out of print. Maybe Topology by Munkres is the standard text. Introduction to Topology by Gamelin & Greene is a good value.
³Famously concise and overpriced. Some prefer the more expansive treatment of Mathematical Analysis by Apostol.
⁴Effecient treatment which avoids matrices and determinants when possible. Linear Algebra by Curtis is not so finicky and has additional topics: rational canonical form, tensors, principal axis theorem, Hurwitz’s theorem.
⁵Statistics is not math, but it feels like it and everybody should know some.
^5aIncludes proof of prime number theorem. For a brief, non-analytic course to prepare one for algebra, consider The Theory of Numbers, sections 1.2-1.4, 2.1-2.3,2-6-2.9,3.1-3.7
⁶Meant to be read with a copy of Euclid’s Elements. Contains interesting advances on Greek mathematics: Gauss’s construction of the regular 17-sided polygon. Use of modern algebra to show that constructions such as trisecting an angle, duplicating a cube, and squaring a circle are impossible with just straight edge and compass. Hilbert’s extra axioms to make Euclidean geometry completely rigorous. Non-Euclidean geometry.
⁷Eventually one wants to read Algebraic Geometry by Hartshorne but this text is more accessible.