A probability distribution for IT project planning and the question that has answer 42

Reading The Black Swan by Nassim Nicholas Taleb (page 159) I came across the idea of looking at the conditional expectation of a random variable \(X\), given that \(X>k\) for some \(k\). If you life expectancy is 80 years, but you are already 87 years old, the conditional life expectancy is now, say 92 years. In math notation \(E(X|X>87)=92\). Look here for real data. The point here is that the more extreme the condition, the less the added additional expectation. So for example, \(E(X|X>119) \le 120\).

Now think about an IT project. A similar project took 18 weeks to complete. Now already 40 weeks have passed. We are interested in \(E[X|X>40)\). Anyone who has worked in IT knows that the new estimate for this project is unlikely to be 41 weeks. To get any further, we need to make assumptions about the probability distribution of \(X\). If \(X\) has a Poisson distribution with mean 18, it turns out \(E[X|X>40)=40.7\) so we are obviously looking for a wilder distribution. The wildest distribution I know is the Pareto distribution (also called power law).

The Pareto distribution takes two parameters, the minimum value and a shape parameter \(\alpha\). Let's say the shortest time a project can take is 1 (week or any unit of time). Say, from experience, that 1 estimated week in a few historical projects actually turned out to be 1, 3, 7, 2, 1, 2, 1, 4 weeks. These number are not outrageous - I expect you to sigh out of boredom: This is just the nature of IT projects. Ok, then we can use the maximum likelihood method to fit the shape parameter of the Pareto distribution. That's a oneliner in Mathematica. It comes out at \(\alpha=1.375\).

Now we can generate some random numbers out of this and see what we can expect in future projects! I get 1, 1, 1, 1, 2, 2, 3, 3, 7, 70. That's not exactly matching my intuition. We get something close to the expected value most of the time, with an occasional wild outlier. I'd rather see something a bit more in between the expected value and wild outliers like 70. Before moving on, let's summarize in Mathematica what we did so far:

I want a distribution somewhat like the Poisson distribution below, but that is fatter in the right-end tail, something that falls off less quickly.

The Gamma distribution comes to mind. It's density function has a polynomial factor like \(x^2\) that dominates when \(x\) is small and an exponential factor \(e^{-\beta x}\) that dominates for large \(x\), ensuring rapid fall-off of probability in the right tail. To attenuate the fall-off, let's take the square root of \(x\) in the exponential. We get our candidate distribution for project planning

$$f_0(x):=\int_0^\infty x^2 e^{-\sqrt{x}}$$

Well, the total probability has to be 1 so we must normalize it. I suppose you can make the substitution \(t=\sqrt{x}\) and grind it with integration by parts. But Mathematica can do it symbolically and it comes out as 240. So we actually have

$$f(x):=\frac{1}{240}\int_0^\infty x^2 e^{-\sqrt{x}}$$

Let's look at a plot of it!

Ok, that is somewhat promising. What is the expected value? Later on, we can add location and shape parameters, but let's just consider this "canonical" version first. Again trusting Mathematica on the details, lets just enter the definition of expected value \(\int_0^\infty x f(x)\) and press enter. It comes out as

That's right - the answer to the Ultimate Question of Life, the Universe, and Everything! Wow. After the initial shock, let's try to generate some random variates from this to see if we get results that matches our experience from IT projects. For that, we can pass in uniform \([0,\,1]\) random numbers into the inverse of the cumulative distribution function, that's a well known result in statistics. I get 10, 10, 11, 22, 23, 28, 43, 60, 61, 170. The most probable outcome is around 20 at the peak of the density function. But the expected value is 42 since occasional tail outcomes raise the average with nothing compensating from below. This is maybe not a great distribution to model IT projects, but it has the nice property that you can legitimately hope for getting the 20 outcome but should expect 42 weeks until completion...

This blog post is already too long. I will leave you with the last piece of Mathematica code.

What distribution do YOU use to model your IT projects? Comment below!

Optimal bus size with self-driving vehicles

Ok, back after a seven year break from the blog. Today when driving my car behind a bus after having read The Principles of Product Development Flow by Donald Reinhertsen I realized that we will maybe not have large buses in the future of driverless cars.

The reason is that the optimum batch size decreases when the transaction cost decreases. That is, the optimum number of items (=persons) per batch (=vehicle) is a function of the transaction cost. The transaction cost is cost of gas, maintenance and not the least -- the driver.

In figure 5-4 in the book, we can see the deduction for the optimum batch size Q is found at

   Q = sqrt(2F*N/H)

where F is the fixed cost per batch, N is the number of items per year and H is the holding cost per year. The holding cost is the cost of having people waiting for reaching their destination. It's hard to estimate so we just let it be. The transaction cost F is not so hard to estimate, let's say the cost of gas and maintenance for driving half an hour is $8 and the cost of the driver (including HR, average sick leaves etc) is $80. Then we can cut F to a tenth by having a driverless vehicle. So the new optimum batch size becomes

  Q' = sqrt(2*(F/10)*N/H) = sqrt(1/10) * Q  ~ 32 % of Q

So if a bus takes 60 people, the new optimum vehicle will take around 20 people. The closest one gets to that is a rather small bus.

Hopefully I will check in my blog in 20 years or so to see how my prediction turned out :)

A nice way to visualize a solution of the heat equation on a three dimensional cube

Having nothing to do on the train today, I fired up Mathematica with the intention to solve a differential equation on a three dimensional domain just to see how well it worked. It turned out to be easy as π. Just type in the equations. And better yet... With the slightest amount of Mathematica magic, I was able to visualize the solution for a fixed time.

Since we are not four dimensional (or more) beings, it is hard four us to look inside three dimensional objects. Unless, of course, they have some kind of openings. An easy way to open up a cube is to make it into lots of small cubes, like this.

It is just a lot of Cuboid's. The color function which I called heatColor is essentially a wrapper around RGBColor. It has two parameters min and max, for temperatures T close to min we choose litte red and a lot of blue. Vice versa for T close to max.

In order to solve the partial differential equation I use NDSolve. The boundary conditions are zero on all sides of the cube except the bottom (z=0). I just chose something that would be consistent with the other boundary conditions as well as with the inital condition. The inital condition is also chosen as anything simple that will be zero when x, y or z get close to zero or one. The result from NDSolve comes in the form of an interpolated function, so you can evaluate the numerical solution in any point that you wish!

You can play around with the number of cubes and the opacity to find your personal optimum.

The code:

solution =
NDSolve[{D[u[x, y, z, t], t] ==
D[u[x, y, z, t], x, x] + D[u[x, y, z, t], y, y] +
D[u[x, y, z, t], z, z],
u[0, y, z, t] == 0,
u[1, y, z, t] == 0,
u[x, 0, z, t] == 0,
u[x, 1, z, t] == 0,
u[x, y, 0, t] == t x y (1 - x) (1 - y),
u[x, y, 1, t] == 0,
u[x, y, z, 0] == x y z (1 - x) (1 - y) (1 - z)},
u,
{x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {t, 0, 2}]


(* Let v be the first - and only - solution *)
v = (u /. solution)[[1]]


heatColor[T_, min_, max_] :=
With[{c = (T - min)/(max - min)}, RGBColor[{c, 0, 1 - c}]]


Graphics3D[Table[
{heatColor[v[i/5, j/5, k/5, 1], 0, 0.03],
Opacity[.50],
Cuboid[{i, j, k}/5, {i, j, k}/5 + .10]},
{i, 0, 5}, {j, 0, 5}, {k, 0, 5}], Axes -> True,
AxesLabel -> {x, y, z}]

Programming for multi-core processors

CPUs come with an increasingly number of cores. As a programmer, you just have to relate to this in some way. One way would be to make sure you really understand parallel programming to utilize the higher responsiveness and speedup made possible with multi-core CPUs. But parallel programming is complicated and it's quite common not see the speedup you hoped for. Here's an insight I just had:

If you parallelize your code yourself, you will also have to maintain and improve it yourself. Conversely, if you rely on automatic programming constructs, the provider will continue to improve the parallelization for you!

For instance, by using the Cilk extension of the C programming language to expose parallelism in your program, you can recompile it with an improved version of Cilk to get a better speedup. If you would have used for instance MPI, you would have to reprogram your code by yourself to get a faster program. We both know it; it's not gonna happen, it's way too boring to deal with old code :)

With this argument in mind, I think the only reason to learn about parallel programming is because it's fun. It's not necessary in order to remain competitive in a multi-core future. I guess this calls for a blog post on different ways to exploit multi-cores CPUs without having to do the actual parallelization yourself...

Bézier Curves in Mathematica

Picking up Applied Geometry for Computer Graphics and CAD by Duncan March, I read a chapter on Bézier curves. They are polynomial curves defined by control points that in a very intuitive way let a designer shape the curve. Bézier curves are really a killer app for the Manipulate and Locator commands in Mathematica. A cubic polynomial need four Locators (control points) in order to be defined:

You need a cubic polynomial in order to be able to create a loop or a cusp:

Is there any other interesting phenomena that becomes possible if we allow higer order polynomials? I don't know but here's a heart or the contour of a fox head. The LocatorAutoCreate option let us add Locators just by holding the alt key and click in the graphics area.

To do this in Mathematica (version 6 or more), define the Bernstein polynomials and the Bézier curve as on page 131 in March's book. Then in the Manipulate construction, you can make a list of Locators called b. Here's the code:

bernstein[i_, n_, t_] := n! (1 - t)^(n - i) t^i/((n - i)! i!)

bezier[t_, b_List] := 

 Module[{n = Length[b] - 1}, 

  Sum[b[[i + 1]] bernstein[i, n, t], {i, 0, n}]]

Manipulate[

 Module[

  {bez},

  bez[t_] = bezier[t, b];

  ListLinePlot[Table[bez[t], {t, 0, 1, 0.01}], 

   PlotRange -> {{0, 1}, {0, 5}}]

  ],

 {

  {b, Table[{i, 0.0}, {i, 0, 1, 1/3.0}]},

  Locator,

  LocatorAutoCreate -> True

  }

 ]

It is crucial that the bezier function takes an arbitrary length list b as parameter. Otherwise it wouldn't work to use LocatorAutoCreate. Consider doing this in Java, how many cups of coffee you would need. To finish, in version 7 you can use an in-built function called BezierCurve instead of defining your own. But it's not necessary from a performance perspective, the response is immediate with my version on a single-core 1.7 MHz Pentium M.