Showing posts with label math. Show all posts
Showing posts with label math. Show all posts

Jun 24, 2010

If Cantor Were A Programmer...

When you take math in university, one thing you'll probably come across is the concept of countability. This means that given a set, you can find some way of assigning a natural number to each element of that set - you can count the elements in it.

This guy Cantor discovered in the late 19th century that there are different types of infinities. He did this by analyzing various sets of numbers and discovered that certain infinite sets are in fact "bigger" than other infinite sets. Which is kinda weird, but it works.

I thought it might be fun (my definition of fun is apparently kinda weird) to see what would happen if instead of a mathematician Cantor were a programmer, and what his proof might look like.

Let's take a hypothetical language called C∞, which is kinda like C except that it runs on our magical computer that has an infinite amount of RAM, and thus the integer types in C∞ have infinite precision.

In C∞, we say a type T is countable if we can construct a function like this:
T foo(unsigned int n)
Where foo() satisfies:
1) For every possible unsigned int value n foo() will return a T - this means no exceptions, no crashes, etc.
2) For every possible T value t, there is exactly one unsigned int n1 such that:
foo(n1) == t
We call functions like this one-to-one.

What are some countable types? Well, signed int is countable. Here's a function that will work:
int foo(unsigned int n){
  if (n % 2 == 0){
    return n / 2;
  }else{
    return -(n + 1) / 2;
  }
}
If you start from 0 and continue, this function will produce 0, -1, 1, -2, 2, -3, ...
This turns out to be a one-to-one function, so therefore the signed ints are countable.

It also turns out that signed int pairs are also countable. Suppose we have a type that looks like this:
struct pair {
  int a, b;
};
This type is also countable. The proof is fairly simple when you use a picture, see here for a nice one that is already done. It uses rational numbers, which are a subset of pairs of integers.

What about the real numbers though? This one is a bit trickier. Let's represent a real number like this:
struct real {
  int integerPart;
  bit decimal[∞];
};
What we're doing here for the decimal part is representing it as an infinite series of 1s and 0s. Since we're all computer geeks here, we know that a binary string is just as capable of representing a number as a decimal string, so we still are able to represent any real number using the above type.

Now we show that the reals are not countable using a contradiction. First, we assume that the reals are in fact countable. Let's simplify it a little and only work with the real numbers between 0 and 1, so we'll have a real0 type that is the same as real but with the integerPart field set to 0.
So since real0 is countable (by assumption), it means that we can create a one-to-one function like this:
real0 foo_real(unsigned int n)
Since this function is one-to-one, for every possible real0 value x there is an unsigned int n that will produce x when you call foo_real(n). So let's create an array that has all those real numbers in it:
real0 all_the_reals[∞];
for (unsigned int i = 0; i < ∞; i++){
  all_the_reals[i] = foo_real(i);
}
Now let's cause a problem. Consider the following code:
real0 y;
for (unsigned int i = 0; i < ∞; i++){
  y.decimal[i] = !all_the_reals[i].decimal[i];
}
Why is this a problem? Well, we assumed that all possible values for real0 are in all_the_reals. But for any i, we know that y != all_the_reals[i] because the bit at the ith decimal place is inverted. Therefore we know that y is not in all_the_reals, which means that all_the_reals does not contain all possible real0 values. Therefore there is a real0 value that cannot be produced by foo_real() — a contradiction since we assumed foo_real() is one-to-one. So real0 is not countable, and since real0 is a sub-type of real, real is also not countable.

What does this tell us? Well, we know that with our infinite precision computer there are an infinite number of possible unsigned ints we could have. At the same time we have shown that there is not a one-to-one mapping between the reals and the unsigned ints, so therefore the infinity of the reals is effectively "bigger" than the infinity of the unsigned ints.

Why would a programmer care? Well, this has an implication for computer science. But first, let's ask another question - if there are real0 numbers missing from all_the_reals, how many are there? Is there just one? Or even a finite number? Well, no, and here's why. Suppose there are a finite number of missing real0s. We could then create an array of them:
real0 missing_reals[N];
However we could then still create a one-to-one function:
real0 foo_real_with_missing(unsigned int n){
  if (n < N){
    return missing_reals[n];
  }else{
    return foo_real(n - N);
  }
}
So there must be an infinite number of missing real0s. Are the missing ones here countable? Assume that they are, and create another one-to-one function foo_real2() which gets all these missing real0s. But then we could still create a one-to-one that will grab all the reals returned by foo_real and foo_real2:
real0 foo_real_combined(unsigned int n){
  if (n % 2 == 0){
    return foo_real(n / 2);
  }else{
    return foo_real2((n - 1) / 2);
  }
}
So even the real0s that are missing from foo_real() are uncountable. Crap! In fact based on this we know that even the reals missing from the combination of foo_real() and foo_real2() is uncountable. If we create another one foo_real3() which grabs the real0s outside of foo_real_combined(), it is still uncountable. No matter how many functions we create, they will still not capture all the real0s. What does this have to do with programming? Suppose that our magical computer not only has infinite RAM, but infinite processors. The ith processor will compute foo_real(i) in the first clock cycle, foo_real2(i) in the second, etc. But since the real0s are uncountable, even with infinite processors we wouldn't be able to enumerate them all in finite time. That's something a programmer might care about.

Jun 7, 2010

Discovering Fractals: The Discrete Logistic Map

I've discovered another fractal for you, although this one doesn't have a really pretty picture like the other ones do. Instead the beauty of this one is in how weird the behaviour is. It is a very simple equation, yet can still manage to exhibit chaotic behaviour.

Here's the equation in C:
x = a * x * (1.0 - x);
Very simple. The variable x is some variable in between 0 and 1, and a is a fixed constant between 0 and 4.

The interesting part of this is when we execute this equation over and over. For small values of a the result is kinda boring, it converges to some fixed point and stays there. However once you start increasing a, some interesting things pop up. Here is a graph with a = 2.9..4 on the x-axis, and x = 0..1 on the y-axis:
At a = 3, the system no longer has a single fixed point, and begins to oscillate between two points. And around 3.45 it forks again into 4 points, and so on. Eventually you move into regions where it is chaotic, and doesn't really settle down. However even more interesting is how the chaotic regions break up into non-chaotic regions temporarily, and then go back to being chaotic. If you zoom in on these non-chaotic intervals you end up seeing that the region looks the same as the larger picture - which is why this is a fractal, it shows self-similarity.

Why would a programmer care about this? Well it's an example of a simple, completely deterministic discrete system that does not always stabilize, depending on the values of the parameters. It's possible for us to create much more complicated software systems that maybe have the same property - if we run the system 1000 times to see how it behaves, the 1001st (th?) time may still have completely different results. Same for the 1 000 000th and 1 000 001st times. The picture above is with 1000 iterations, but when I up it to 1 000 000 iterations the picture still looks the same (just takes longer to render). Because parts of this system are chaotic, no matter how many times you iterate it will not converge.
In short, an awareness of chaotic systems might be helpful when testing your programs.

If you want to see the code that generated image I have put up a gist that you can look at.

Apr 26, 2010

Discovering Fractals: Brownian Trees

After I got back into generating fractals yesterday, I decided to do up another interesting one. This one is called the Brownian Tree. The difference between this one and all the other ones that I've done is that this one is stochastically generated, where the other ones are all deterministic. This means that each image is random, so I can't tell you how to generate each one exactly the same - although technically since computers use pseudorandom numbers and not real random numbers, I could just give you the seed and you're set. However I didn't record the seeds.

The way this one works is based on something in reality. You start with a world. Within this world, you fix a seed particle. Then you repeatedly add new particles to the world and have them float around. When the new particle bumps into the seed it becomes part of the seed.
In pseudo-C, it would look something like this:
bool world[SIZE][SIZE]; // assume it is pre-set to false
world[rand() % SIZE][rand() % SIZE] = true;

for (int i = 0; i < NUM_PARTICLES; i++){
particle = [rand() % SIZE, rand() % SIZE];

while (true){
projection = particle + random direction
if (projection out of bounds){
// do something
}if (world[projection] == true){
world[particle] = true;
break;
}
particle = projection;
}
}
plot(world);
What this ends up doing is generating structures that look really organic. Here's an example:It looks like a shrubbery!

So two things to note about the algorithm. When the particle goes out of bounds, I just put it in some random other spot in the world. Also when the particle collides with the structure, I keep track of how long it took to get there (it resets to 0 when it goes out of bounds). Based on how long it takes, I give it a different colour. This leads to the nice layering effects that you see there.

You can change this algorithm in quite a number of ways. Here's a modification where when it bumps into the side of the window, it just sticks there:Or one with multiple seeds:Or one where instead of using a point for the seed, you use a collection of points in a ring (I stole this idea from the Wikipedia page):One thing I read briefly in my searching is that you can also generate music using fractals. This would be quite interesting, and would force me to learn how audio files work! Maybe I will post something on it at some point.

Discovering Fractals: The Lorenz Attractor

It's been a long time since I posted anything about fractals. I've been looking at them a bit more recently and have decided to put up more pretty pictures for you all. Today we're looking at the Lorenz Attractor, which is one of the earliest discovered chaotic systems - technically many systems had been discovered before that, but this is one of the ones that was discovered after someone (that someone being Lorenz) had figured out what chaos was.

It's a fairly simple system, although less simple than the previous ones I've written about. It follows a set of differential equations(I just grabbed these from the Wikipedia page):
dx/dt = σ * (y - x)
dy/dt = x * (ρ - z) - y
dz/dt = x * y - β * z
If we translate these to C (I do all my fractal stuff in C, since it seems to be hopelessly slow in anything else) we get:
for (t = 0; t < NUM_ITERATIONS; t++){
x1 = x + dt * sigma * (y - z);
y1 = y + dt * (x * (rho - z) - y);
z1 = z + dt * (x * y - beta * z);

x = x1;
y = y1;
z = z1;

// plot point
}
You'll want some initial conditions, I used x = 0.1, y = 0.0, z = 0.0. You also have a time-step dt, which I set to 0.001. This determines how "smooth" your picture will look.

Here's a picture of what it looks like:This is with NUM_ITERATIONS = 100 000. I just used a mostly orthographic projection - the x just maps to x, and y maps to y when you plot it on the screen, and you ignore z. I did a bit of scaling and translating to the x and y so that it fit nicely into an 800x800 window.

The interesting thing is the colours (it may look ugly, but there is some science behind it!). This picture is effectively a trace of a particle moving through time. The more red it is in this picture, the earlier it is in the particle's path, and the more blue, the later (obviously purple ones are in the middle). It's interesting because except for the little red swirls at the beginning there, the red, purple and blue paths are all fairly mixed up (it might look like the blue ones are more clustered, but keep in mind that blue pixels are plotted later, and therefore will overwrite any pixels that might have been there before). This is evidence of something called topological mixing, which as I understand from my brief digging through Wikipedia basically means that the paths taken at different time ranges will inevitably overlap.

Anyway, that's enough about the actual math behind chaos theory. You can probably even take on Jeff Goldblum now! I'll work on other fractals sometime and give you guys more neat pictures.

Jan 26, 2010

Max/Min Problems in LaTeX

Suppose you want to write a maximization or minimization problem in LaTeX, and you want the variables that you are maximizing with respect to under the "max" keyword, kinda like this:

Unfortunately out of the box LaTeX can't do this, but there's a handy little trick. Right before the \begin{document} put:
\usepackage{amsmath}
\DeclareMathOperator*{\Max}{Max}
\begin{document}
Then to write the max problem above you can put:
V(p, y) = \Max_x U(x) : p \cdot x \le y

Nov 17, 2009

Learning Scala: Euler's Method

I've been messing around more with Scala to see some of its more interesting features, two of which are pattern matching and currying. These are really handy in certain situations, and I decided to share them with you today.

The example I'm using is the Euler method for approximating a solution to a differential equation. This is probably not the most intuitive way of doing the Euler method, but it uses a lot of Scala's interesting features and I think it is a good way of illustrating them.

The Euler method works like this. You have a differential equation in the form y' = f(x, y) and an initial value (x0, y0). You calculate the slope at the initial value, which is f(x0, y0). You then take a step of size h in the x-direction following the slope, which brings you to (x0 + h, y0 + h * f(x0, y0) ). This becomes your new point, and you repeat until you are satisfied with the results. A smaller h will lead to more accurate results, but with more computation.

Here is the not-so-intuitive code:
object Euler{
def main(args: Array[String]){
// Use the differential equation y' = x^2 + y^2
val my_de = solve_de((x, y) => Math.pow(x, 2) + Math.pow(y, 2)) _

// Use a step size of 0.1, and do 10 steps
val generator = my_de(0.1, 10)

// try with different initial conditions
println(generator(0, 1))
println(generator(1, 1))

// create another generator with a smaller step size
val precise_generator = my_de(0.05, 20)
println(precise_generator(0, 1))
println(precise_generator(1, 1))
}

def solve_de(func: (Double, Double) => Double)(step : Double, iterations: Int)(x_0: Double, y_0: Double) =
(List(List(), List(x_0, y_0)) /: (0 until iterations))((s, i) => s match {
case List(res, List(x, y)) => {
val fxy = func(x, y)
List(res ::: List(fxy), List(x + step, y + step*fxy))
}
}).head

}
Look first at the definition for solve_de. It has three argument lists! One of them takes a function of type (Double, Double) => Double, also known as a function that takes two Doubles and returns a Double. The second takes a Double and an Int, and the third takes two Doubles.
The first parameter list takes a function representing the differential equation we are approximating, the second takes the parameters for the approximation which are the step size and the number of steps we take, the third takes the initial values.

If you look back at the main function, you'll see that we call solve_de with only one of its argument lists, and a _ at the end! This is called currying - it returns a new function with the first parameter list filled out. We save this as our variable my_de (note we use val, so this variable is immutable) which is a solution generator for our specific differential equation described above. You use the _ at the end to tell Scala that you are only partially applying the solve_de function.

Next, we call my_de with the values 0.1 and 10. This creates another function which solves our differential equation using a step size of 0.1 and 10 steps. We can then call this function with different initial conditions to get different solutions to the differential equation. Each time we call it, it returns a list of 10 points that lie along our solution curve. If we wanted to, we could then plot this curve using some graphing library.
Note: for some reason here you don't have to call my_de with the _ at the end, it is probably for some reason that I do not yet understand.

After that, we create a more precise generator with half the step size. I double the number of steps so that if we were to graph this alongside the first list, they would have the same range of x-values.

The next bit is the solve_de function, which illustrates some of the more interesting features of Scala. First one (which probably isn't that interesting) is that there are no curly brackets around the body of solve_de. If you have a function in Scala that is only one line, you can just write:
def foo(x) = ...
You don't need to include curly brackets.
We have a fold using the /: operator. If you've used inject() in Ruby then you'll know what I'm talking about, otherwise take a look here for a description of what fold (aka reduce1) is. In Scala you can write this:
(0 /: myList)(some function f)
This does a left fold of myList with the initial value 0, using the function f - aka f(f(0, myList[0]), myList[1])....
The initial value is a list that looks like this:
[ [], [x0, y0] ]
The first element of this list is where we will be sticking the approximated values, the second element is our current point.
We fold this list into (0 until iterations), which is a range equivalent to 0..(iterations - 1) in Ruby. This is a very interesting piece of Scala, because it shows some of the fancier features. Scala is a pure object-oriented language so the 0 there is actually an object of class Int. Effectively what we are doing is calling 0.until(iterations), which returns a Range object that we can use fold on. In Scala for a method which only takes one argument, you don't need to put the . or the brackets.
However there is no until() method for Int. Where does this until() come from? Scala has a feature called implicit functions, which are used for implicitly converting one type into another - like the auto-boxing between int and Integer in Java. Little do you know, there is actually a class called RichInt which supplies the until() method, and a bunch of other handy things (you could write 0 to 5 if you like). When you call until() on 0, Scala first looks to see if Int has an until() method. Since it doesn't, it checks to see if there is an implicit conversion for Ints into a class that does have an until() method. Since there is only one such class (RichInt), it automatically replaces your statement with something like toRichInt(0).until(iterations). If there were more than one implicit conversions however, then Scala would give you a compile error and you would have to explicitly provide your conversion. The main difference between this and auto-boxing in Java is that you can provide your own implicit conversions between any classes you like, provided they don't result in ambiguities.

The next step is to provide a function to the fold operator to use for folding. After the => we see
s match {
This matches the variable s (the "accumulator") against a set of patterns. This is another feature of Scala called pattern matching. This example doesn't really do it justice since we only have one pattern here, and it is just so that we can have a nice way of extracting the variables out of s without using head() and tail(). I think I might post something more detailed on pattern matching in the future. Anyway, we use the expression List(res, List(x, y)) to match s, and this extracts out our current accumulated values as res, and the current position into x,y. We can then compute f(x, y) and put it in fxy (this is to save some time in computation) and then return:
[ new res, [ x + step, y + step * f(x, y) ] ]
The new res value is just res with f(x, y) stuck on the end (that's what the ::: operator does, it concatenates two lists).

Two small syntactic things to note:
- There is not a single semicolon in this program. Scala doesn't need the semicolons at the end of lines, although you can include them if you like.
- There are no return keywords in this program, even though we have functions. Scala doesn't require the return keyword, it will insert it where it thinks you are trying to return something.

So I'm not sure if I'd recommend you actually write Euler's method like this, instead you would probably write it something like this in Scala:
def solve_de(func: (Double, Double) => Double)(step : Double, iterations: Int)(x_0: Double, y_0: Double) = {
var x = x_0
var y = y_0
var res = List[Double]()

for (i <- 0 until iterations){
val fxy = func(x, y)
x = x + step
y = y + step * fxy
res = res ::: List(fxy)
}
res
}
However in this case, you wouldn't be able to use all those fun little toys that Scala gives you, so I did it in a different way.

1In Scala fold and reduce are two different things: fold takes an initial element, where reduce uses the first element of the list as the initial element. Reduce will throw an exception if used on an empty list, where fold will just return the initial element. In non-Scala languages, reduce and fold are the same thing.

Jul 18, 2009

Fractal Time Analysis Results

In my last post about fractals, I wrote up about how I wanted to analyze which seeds would take a long time to render and which would be fine.

The process goes like this:
real, im = -1, -1
times[IMG_SIZE][IMG_SIZE]
x, y = 0, 0

while x < IMG_SIZE
y = 0
while y < IMG_SIZE
render

times[x][y] = processing time
y++
x++

normalize times // min will be 0, max will be 255
output to bmp
So what you end up with is a greyscale image. The brighter it is, the longer it took to render.

The result is...well, amazing. Like really. See for yourself:



Does that look familiar to anyone?

I used a fairly small image size here (it was even smaller before). Problem is that rendering time is O(n2), so it takes 4 times as long when I double the size - originally I wanted to have an 800x800 image - turns out that'd take like a week for the thing to process! I might try tweaking the code to be multithreaded, that would probably half the time on my machine.

Jun 20, 2009

Julia Sets - Moving the seed

I've been at the fractals again. This time I made it morph:


This is a Julia set again, but instead I am tracing a path through the complex space with the seed value. The path taken in this video is a near-circle centred at 0.12 + 0.74i with an real radius of 0.11 and an imaginary radius of 0.10.

There were a number of other really cool ones, but the problem is that the path they take tends to go into areas which have lots of points that do not fly off into infinity, which means they take a long time to process. That means that I can't have a nice framerate like the video has.

One thing that was kinda neat was that if you flip the sign on the centre of the ellipse for the imaginary component, it flips the patterns displayed on the x-axis (or maybe it's the y-axis). So centred at 0.12 - 0.74i, those swirls are in the top-right and bottom-left instead of the top-left and bottom-right.

A thing I would like to try to make better animations would be to see where the framerates get low. I can probably do this by iterating between -1 and 1 in the real and imaginary components and spit out a time for each computation. My plan will be to spit that out into a greyscale image so that it is easy to see where the slow points are.

Jun 5, 2009

Newton's Iterative Method

Once in a while when you're doing math processing, it is useful to be able to estimate the zeroes of a function (if f(x) = 0, then x is a zero of f). An example is square roots. A square root of n is simply the zero of f(x) = x2 - n. This is a really useful operation.

I'm going to talk about Newton's Iterative Method, which is a way of estimating a square root. I say estimating because since most square roots are irrational numbers, it is impossible to represent them accurately using floating point arithmetic. So we just get as close as we can to it.

There are other ways of calculating square roots. The square root of n can also be expressed as eln(n)/2. However this is slower to calculate.

How does this method work? Well to understand it fully you need to understand calculus, however with square roots the calculus is really simple so you can just take my word for it if you don't know calculus.

Let's try to find the square root of 2. Our function is therefore:
f(x) = x2 - 2
The derivative1 f′ of this function is 2x.

We start off with a guess. Let's say 1. If you want you can start off with any other number except zero (I'll explain why you can't use zero in a bit), but I'm going to start with 1.

What we want to do now is find the linear approximation (first order Taylor series if you want to be more pedantic) of f around the point x = 1. We'll call this point x0 The function, we'll call it T, for a linear approximation is:
T(x) = f(x0) + f′(x0) * (x - x0)
For our function, it looks like this:
T(x) = x02 - 2 + 2x0(x - x0) = 2x0x - x02 - 2
At x0 = 1, we have:
T(x) = 2x - 3

Now for the more interesting part. Let's look at a picture of these two functions (the blue line is the x-axis):

We want to find out where the linear approximation hits the x-axis. By setting T to zero and solving for x we get:
x = x0 - (x02 - 2) / 2x0 = 1.5
We now use this 1.5 as our new x0, except that we will call it x1. Now we plug that into our formula:
x = x1 - (x12 - 2) / 2x1 = 1.416667
Call that x2, plug that in again, we get:
x = x2 - (x22 - 2) / 2x2 = 1.414216
We're getting pretty close eh?

When do we stop looping? Well, since we can't get absolute precision, we can stop looping when the xi is within a certain threshold of some error:
while |f(xi)| > ε
do iterative method
You can set this ε to whatever you want. A higher value will mean less precision but more speed - however the iterative method moves pretty damn quickly toward the zero so it is not a huge deal.

There's a few catches. What would happen if we used a guess of zero? Well, we'd end up with a divide-by-zero error. Basically the linear approximation to f would be flat, and never touch the x-axis.
What would happen if we used a negative guess? Well, that leads to a more interesting discussion. Newton's Iterative Method finds local zeroes. The closest zero to a negative number (say -1) is not around 1.41421. It is close to -1.41421, which is the other square root of 2. When there are multiple zeroes to a function, the iterative method will tend toward which ever one it is sloping toward. So if you're trying to find the zero of a function that crosses the x-axis several times, choose your initial guess wisely!

Finally, here's an interesting problem. Try finding the square root of -1 using this method. Try it with different guesses.

1If you don't know what a derivative is, it is the function which tells you the slope of f at any given value of x.

Mar 30, 2009

Discovering Fractals: The Julia Set

After my post yesterday about the Mandelbrot set I started digging around a bit more. I discovered this other thing called a Julia set, which is something similar.

It still uses this formula:
zi+1 = zi + c
However where Mandelbrot set started with z0 = 0 and c as the coordinates of the pixel, the Julia set has z0 as the coordinates of the pixel and c as some number. Depending on the c that you pick, you end up with a different picture.

I also looked into colouring it a bit better than black and blue, it is now based on HSV colours instead of RGB. Then if the formula escapes your radius at iteration i, you set the hue of the pixel to (i+b) mod 360 where b is some offset you can choose to change your colour scheme. With b = 0 you get a bright orange background, but I didn't really like it so I picked b = 180 which gives a nice lightish blue background. Unfortunately I'm not completely liking the bright pink/purple of the swirls, anyone know a few things about colour theory and have suggestions?

Here it is centred at 0.35 - 0.33i, with 4x zoom and c = (φ - 2) + (φ - 1)i, φ being the golden ratio:

Mar 29, 2009

Discovering Fractals: The Mandelbrot Set

I read a blog post a while back about content generation in video games, one section of which is terrain generation. One interesting way to generate realistic terrain is using fractals. One of the most well-known and relatively basic fractals is the Mandelbrot set. After reading through this article I decided to write a simple program that generates pictures based on this set.

Here's how the set works. You take a complex number c. You then take a sequence described by the following formula:
zi+1 = zi2 + c
Where z0 = 0. If this sequence is bounded (it does not escape to infinity), then c is in the set, otherwise it is not.

You can display this using a program by taking each pixel in your window and converting it to a complex number. You then test to see if the number is in the set by iterating over it a certain number of times to see what happens. Of course you won't be 100% precise, but with enough iterations you can get a pretty good approximation.

On top of this you can make it look even prettier by tweaking the colour depending on how quickly each point escapes to infinity. If it shoots off into infinity immediately, you make it darker. If it takes a while, then you make it lighter. I used the following formula:
colour = (0, 0, i / n * 255)
Where i is the iteration that the number escaped a certain threshold where the norm of the complex number is greater than this threshold (I used 2 and 3, didn't notice a huge difference between them) and n is the number of iterations (I used 500, 1000 and 2000, 1000 made the prettiest picture). This means that the colour of each pixel ranges between black and bright blue.

Here's what I got when I calculated it centred at c = 0.38 + 0.1i, with 64x zoom (where 1x zoom is in (1 + i to -1 - i)), 1000 iterations and a threshold of 2 - unfortunately when viewing the picture in the browser (at least for me on Ubuntu with Firefox) it doesn't look as nice as if you actually download the thing so I recommend you do that to appreciate it more:


That's pretty eh?

One thing that I find neat about this set is that it is self-similar. Each one of those little blue blobs is how the picture looks if you centre at 0 + 0i and have 1x zoom. You can then zoom into each one of those and see little bumps on the side that also look exactly the same as the whole. And each of those bumps has smaller bumps on it that also look the same. Pretty cool eh?

UPDATE: I wrote a follow-up about this here.

Mar 23, 2009

OORegress: Another Stats Package

A while back I complained about how statistics functionality in OpenOffice is sadly lacking. Then I discovered that JRuby can fairly easily tie into OpenOffice, so I started thinking I might be able to tie the two together.

I mucked around a bit and rolled up something that can help out. It is a stats program that ties into OpenOffice from the command line, allowing you to enter commands to do statistical things. For example, the following code will run a regression:
regress(Y = X)
It will evaluate the model:
Yi = α + β * X + ei
Note that you don't specify the intercept, it is implied. I'll be adding a feature eventually where you can force a zero intercept.
You can also do some fancier things:
regress(ln(Y) = X1^2 + D*X2)
Which will obviously regress
ln Yi = α + β1*X1i2 + Di*X2i + e1
Where X1, X2 and D are independent variables.

The interface follows some simple conventions. The columns of the spreadsheet contain the data, and the very first row contains the name of each column. In your regression equation you address the variables by those names. When you run the regression, the program will open up a new sheet in Calc with the regression output with a bunch of info about the coefficients, their significance, some properties of the variance of the regression, R2, etc. I pretty much just copied the stuff that Excel prints out when you run a regression because this is what my stats classes want. However I make it easier here since you can use an actual regression formula instead of having to copy-paste columns and apply formulas in the spreadsheet itself.

I'm working on documenting how to use the program, and also working on some new functionality like lagged variables and having ΔYi instead of just some function of Yi. However the regression itself doesn't completely work at the moment, so new stuff will have to wait.

You can check out the code if you like from here: http://code.google.com/p/ooregress/. I should have a more functional version coming out soon.