Wednesday, May 30, 2012

Non-Uniform Quadrature Points in Octave/Matlab

Suppose you are given a bunch of points in the form of a vector x (between "a" and "b") and an accompanying density function p(x)>0. For concreteness, consider

x = [0:0.1:1.0]'; px = 0.3 - (x-0.5).^2;

The top subfigure depicts the density function p(x), and the cumulative density function C(x). The bottom subfigure depicts the  quadrature points "z" (circles), and the lines mark the end of the intervals, using N=11.
This density function is shown in the picture above.

We want to create a N*1 vector "z" of quadrature points that are distributed according to the density p(x). We also want a vector  "dz" corresponding to the width of strip belonging to a particular quadrature point. As an additional constraint, we assume that we want z(1) = a, and z(N)  = b to match the end-points.

The program attached at the end is able to do this relatively efficiently. In fact, I used the following command to generate the picture above:

[z dz] = GridDensity(x,px,11,1);

GNU Octave Code:

%
%  PROGRAM:
% Takes in a PDF or density function, and spits out a bunch of points in
%       accordance with the PDF
%
%  INPUT:
%       x  = vector of points. It need *not* be equispaced
%       px = vector of same size as x: probability distribution or
%            density function. It need not be normalized but has to be positive.
%   N  = Number of points >= 3. The end points of "x" are included
%       necessarily
%       Pt = Optional argument. If present, then some plotting.
%  OUTPUT:
%       z  = Points distributed according to the density
%       hz = width of the "intervals" - useful to apportion domain to points
%            if you are doing quadrature with the results, for example.
%
%  (*) Sachin Shanbhag, March 5, 2012
%  GNU Public License
%
function [z h] = GridDensity(x,px,N,Pt)

  npts = 100;                              % can potentially change
  xi   = linspace(min(x),max(x),npts)';    % reinterpolate on equi-spaced axis
  pint = interp1(x,px,xi,'spline');        % smoothen using splines
  ci   = cumtrapz(xi,pint);                
  pint = pint/ci(npts);
  ci   = ci/ci(npts);                      % normalize ci

  alfa = 1/(N-1);                          % alfa/2 + (N-1)*alfa + alfa/2
  zij  = zeros(N,1);                       % quadrature interval end marker
  z    = zeros(N,1);                       % quadrature point

  z(1)  = min(x);  
  z(N)  = max(x); 

%
% ci(Z_j,j+1) = (j - 0.5) * alfa
%
  beta  = [0.5:1:N-1.5]'*alfa;
  zij   = [z(1); interp1(ci,xi,beta,'spline'); z(N)];
  h        = diff(zij);
  clear beta;
%
% Quadrature points are not the centroids, but rather the center of masses
% of the quadrature intervals
%
  beta     = [1:1:N-2]'*alfa;
  z(2:N-1) = interp1(ci,xi,beta,'spline');

%
% Some plotting if required
%
  if(nargin>3)

    subplot(2,1,1)

    plot(xi,ci,'b-','LineWidth',2,xi,pint,'r-','LineWidth',2);
    axis([min(xi),max(xi)]);
    legend('CDF','PDF')
    title('Visualization')
    xlabel('x');
    ylabel('CDF(x)/PDF(x)')

    subplot(2,1,2)
    plot(z,0.5*ones(size(z)),'ro');
    axis([min(xi),max(xi)]);
    xlabel('z');
    hold on;
    for i = 2:N
      X = [zij(i); zij(i)];
      Y = [0; 1];
      plot(X,Y,'b')
    endfor
    hold off;

  endif

endfunction

Friday, May 25, 2012

MathJax: LaTeX on Blogger - finally!

I am easily excited when I find yet another place where I can use LaTeX syntax to typeset mathematics (Google Docs, OpenOffice, again, presentations using BeamerInkscape).

Finally, it seems, it is not incredibly clunky to write math in Blogger. In fact, it is as convenient as writing it in a native LaTeX document.

The latest avataar of MathJax, while not really new, is new to me. Here are a few links which describe how to go about empowering your Blogger account.

The easiest permanent fix is to add some script in your Blogger template file, as described here.

If you don't like messing with your template, you can take the approach outlined here. Basically, you add some lines to each post in the "html" view, which has math in it.

Here is a mandatory test:
\[\frac{df}{dx} = \frac{1}{2} \left(\frac{ab \textrm{ sech}^2(b \sqrt{x})}{x} - \frac{a \tanh(b \sqrt{x})}{x^{3/2}} \right)\]
Beautiful!

Thursday, May 24, 2012

Google v/s Oracle

In spite of how close it looked like in the media (thanks in part to the FUD raised by paid shills), Oracle essentially had to walk away with peanuts - if that. As usual Groklaw was among the best places to follow the trial.

There were many fascinating parts. One of my favorites was the exchange between Judge Alsup (who could program) and the lead counsel for Oracle, Boies, on "rangeCheck". You could not write better Hollywood screenplay:
Alsup: I have done, and still do, a significant amount of programming in other languages. I've written blocks of code like rangeCheck a hundred times before. I could do it, you could do it. The idea that someone would copy that when they could do it themselves just as fast, it was an accident. There's no way you could say that was speeding them along to the marketplace. You're one of the best lawyers in America, how could you even make that kind of argument?
Boies: ... I want to come back to rangeCheck.
Alsup: rangeCheck! All it does is make sure the numbers you're inputting are within a range, and gives them some sort of exceptional treatment. That witness, when he said a high school student could do it--
Boies: I'm not an expert on Java -- this is my second case on Java, but I'm not an expert, and I probably couldn't program that in six months.
Priceless.

Of course, Oracle being Oracle, they will probably appeal. Here's what Linus Torvalds has to say (full disclosure: on his Google+ account)
Prediction: instead of Oracle coming out and admitting they were morons about their idiotic suit against Android, they'll come out posturing and talk about how they'll be vindicated, and pay lawyers to take it to the next level of idiocy.

Sometimes I really wish I wasn't always right. It's a curse, I tell you.

Monday, May 21, 2012

Some interesting cartoons at Spiked Math

1. Epsilon-Delta: reminded me of times long gone

2. My big fat Greek letter

3.  Hilarious IQ test

Check out the other stuff there too!

Thursday, May 17, 2012

The Checklist Manifesto

I've been a big fan of checklists since as long as I can remember. So much so that a lot of friends and family mock my obsession with them.

In high-school, I had an exam-day checklist (set alarm, carry extra pens, make sure of the time-table etc.). I have about five different types of travel checklists, depending on whether I am flying, driving, going camping, or to a conference etc. I have checklists for writing proposals, and papers and on and on.

You get the idea!

So it was with much delight that I read Atul Gawande's article in the New Yorker called The Checklist a few years ago, and his book called "The Checklist Manifesto" which fleshes out some of the major themes from his article, last month.

While it was "preaching to the choir", I thought some parts were very interesting. I particularly enjoyed the chapter called "Hero in the age of checklists".

It was a worthwhile read.

Here is a link to a recent talk at TED by Atul Gawande.

Monday, May 14, 2012

Commencement Speech: David Foster Wallace

A really nice, and somewhat offbeat, commencement speech by the late David Foster Wallace. I remember enjoying his book "Inifinite Jest".

YouTube also has the actual speech:




Wednesday, May 9, 2012

Shitty Economics

Juvenile but entertaining! (via Barry Ritzholtz)



Monday, May 7, 2012

Puzzle Links

1. Interesting "Martin Gardner" puzzle if you have not tried it before

2. Three puzzles in this post by Tanya Khovanova

3. My wife told me about a "lateral" puzzle she heard in her office. Complete the series:

|||, |CC, C, |C, ?


Thursday, May 3, 2012

Targetted Advertizing

I noticed a curious "coincidence" today.

I received an innocuous looking deal on vacuum cleaner bags from Amazon in my email. It was for the exact make and model of vacuum cleaner that I have owned for about two years now.

Nothing out of the ordinary here. Except for two things.

One, I was just about to run out of bags. I remember, because I made a mental note last weekend. Two, I did not buy the vacuum cleaner from Amazon!

I scoured the Trash folder of my email to see if there were prior deals like this, and if I just happened to notice it this time.  Nothing! Nada!

I recalled this story that appeared in the NYT article "How companies learn your secrets" earlier this year:
... a man walked into a Target outside Minneapolis and demanded to see the manager. He was clutching coupons that had been sent to his daughter, and he was angry, according to an employee who participated in the conversation. 
“My daughter got this in the mail!” he said. “She’s still in high school, and you’re sending her coupons for baby clothes and cribs? Are you trying to encourage her to get pregnant?” 
The manager didn’t have any idea what the man was talking about. He looked at the mailer. Sure enough, it was addressed to the man’s daughter and contained advertisements for maternity clothing, nursery furniture and pictures of smiling infants. The manager apologized and then called a few days later to apologize again. 
On the phone, though, the father was somewhat abashed. “I had a talk with my daughter,” he said. “It turns out there’s been some activities in my house I haven’t been completely aware of. She’s due in August. I owe you an apology.”
If you haven't read the article before (it got a lot of secondary press) you should!

Wednesday, May 2, 2012

User-defined Colors in Grace

Last week, I was working on a manuscript with a collaborators. As we approached the final draft, we decided to make the "colors" in our figures consistent. Unfortunately, we used different programs to make our graphs, and "red" in program X is not the same shade of "red" in program Y.

Grace (the program I use) offers a choice of 16 standard colors, from the usual pull down menu.

However, it allows for endless fine-tuning.

Once you save your graph (as graph.agr), open it up in a text editor. The color-definitions are clustered in lines that look like:
@map color 0 to (255, 255, 255), "white"
@map color 1 to (0, 0, 0), "black"
@map color 2 to (255, 0, 0), "red"
Let us say, we were not happy with the shade of grey available (it does seem a little too light). Look for the line:
@map color 7 to (220, 220, 220), "grey"
and change it to something like:
@map color 7 to (100, 100, 100), "grey"
You can tweak all the other color definitions to suit your taste or requirements.