[Fractint] Probability distribution functions and the Triangle Map

24 Jan 2003

      Just spent a couple of hours with a statistics text. Looked at Kerry 
Mitchell's tent map formulae. Went hnnnnnnnnnn...... pdfs.

My initial idea was to iterate the usual nonlinear transformation with one 
hump that lead to Feigenbaum's discovery of his number. Traditionally, a 
parabola is used, with its height governed by a parameter lambda that runs 
from 0 to 4. But a piecewise linear curve - a triangle, in fact, could be used 
instead and (among other things) leads to exactly the same Feigenbaum 
constant. I thought, hnnnn, what about a Gaussian curve? Okay, it's not 
exactly confined to the interval (0,1), but why should that be any hindrance 
to me? I mean, I'm not even going to be confining myself to the reals, let 
alone a particular finite interval of them. I mean, some of these 
distributions are even supposed to be _discrete_!

Kerry Mitchell had chucked in a user-definable "c" parameter to adjust the 
size of the triangle, and to provide a parameter for selecting elements of its 
Julia space. I left it in for the latter reason, but also tossed in the 
various parameters that are used to describe the various distributions (I 
could have used one of these instead of leaving in "c", but which one?). I 
also put in an offset parameter "o", since many of these distributions are 
centred on the origin, and they're probably more useful at, say o=1. I haven't 
done any massaging of these parameters - they still have their statistical 
significance - so some hunting around may be needed to find something 
interesting (I found a few things during my initial explorations, including 
the ubiquitous almond bread man). The bailout is also user-definable, but if 
none is supplied it defaults to 100. Generally speaking, a full zoom-out may 
be a good start. Starting with values like 1 or 2 may be beneficial (note, of 
course, that for some values the resulting distribution may be undefined - a 
look at the formula should be enough in most cases to pick the bad choices).

There are Julia and Mandelbrot sets for each distribution. The Mandelbrot sets 
for some proved tricky - for those, like the exponential distribution, with a 
critical point at infinity, I've arbitrarily initialised z to 100.

Some I'm not yet convinced have anything to offer. The Laplace distribution 
may even be buggy - I used a hacked version of Sylvie Gallett's Graph formula 
to check that the formulae I'd written had the right shape, but despite this 
laplace_distro_m just doesn't seem to want to deliver.

Morgan L. Owens
"Now on to more important matters. Whoopdy-doo."

gaussian_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	o = imag(p1)
	c1 = 2*sqr(sigma)
	z = pixel, c = p3/sqrt(pi*c1):
	z = c*exp(-sqr(z-o)/c1)
	|z|<bailout
}

gaussian_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	o = imag(p1)
	c1 = 2*sqr(sigma)
	z = o, c = pixel/sqrt(pi*c1)
	:
	z = c*exp(-sqr(z-o)/c1)
	|z|<bailout
}

exponential_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	mu = real(p1)
	o = imag(p1)
	c1 = mu
	z = pixel, c = p3/c1:
	z = c*exp((o-z)/c1)
	|z|<bailout
}

exponential_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	mu = real(p1)
	o = imag(p1)
	c1 = mu
	z = 100, c = pixel/c1:
	z = c*exp((o-z)/c1)
	|z|<bailout
}

laplace_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z = pixel, c = p3*a/2:
	z = c*exp(-abs((z-o)/a))
	|z|<bailout
}

laplace_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z = o, c = pixel*a/2:
	z = c*exp(-abs((z-o)/a))
	|z|<bailout
}

cauchy_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z = pixel, c = p3/(a*pi):
	z = c/(1+sqr((z-o)/a))
	|z|<bailout
}

cauchy_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z = o, c = pixel/(a*pi):
	z = c/(1+sqr((z-o)/a))
	|z|<bailout
}

rayleigh_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	o = imag(p1)
	c1 = sqr(sigma)
	c2 = -2*c1
	z = pixel, c = p3/c1:
	z = z-o
	z = c*z*exp(sqr(z)/c2)
	|z|<bailout
}

rayleigh_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	o = imag(p1)
	c1 = sqr(sigma)
	c2 = -2*c1
	z = o, c = pixel/c1:
	z = z-o
	z = c*z*exp(sqr(z)/c2)
	|z|<bailout
}

rayleightail_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	a = imag(p1)
	o = real(p2)
	c1 = sqr(sigma)
	c2 = 2*c1
	c3 = sqr(a)
	z = pixel, c = p3/c1:
	z = z - o
	z = c*z*exp((c3-sqr(z))/c2)
	|z|<bailout
}

rayleightail_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	a = imag(p1)
	o = real(p2)
	c1 = sqr(sigma)
	c2 = 2*c1
	c3 = sqr(a)
	z = o, c = pixel/c1:
	z = z-o
	z = c*z*exp((c3-sqr(z))/c2)
	|z|<bailout
}

lognormal_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	zeta = imag(p1)
	o = real(p2)
	c1 = 2*sqr(sigma)
	z = pixel, c = p3/sqrt(pi*c1):
	z = z-o
	z = c*exp(-sqr(log(z) - zeta)/c1)/z
	|z|<bailout
}

lognormal_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	sigma = real(p1)
	zeta = imag(p1)
	o = real(p2)
	c1 = sqr(sigma)
	c2 = 2*c1
	z = exp(zeta-c1)+o, c = pixel/sqrt(pi*c2):
	z = z-o
	z = c*exp(-sqr(log(z) - zeta)/c2)/z
	|z|<bailout
}

logistic_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z=pixel, c=p3:
	z = z-o
	z = c*exp(-z/a)/(z*sqr(1+exp(-z/a)))
	|z|<bailout
}

logistic_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	o = imag(p1)
	z=100, c=pixel:
	z = z-o
	z = c*exp(-z/a)/(z*sqr(1+exp(-z/a)))
	|z|<bailout
}

pareto_jul {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = a+1
	z = pixel, c = p3*a/b:
	z = z-o
	z = c/((z/b)^c1)
	|z|<bailout
}

pareto_man (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = a+1
	z = 100, c = pixel*a/b:
	z = z-o
	z = c/((z/b)^c1)
	|z|<bailout
}

weibull_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = b-1
	z = pixel, c = p3*b/a^b:
	z = z-o
	z = c * z^c1 * exp(-(z/a)^b)
	|z|<bailout
}

weibull_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = b-1
	z = o, c = pixel*b/a^b:
	z = z-o
	z = c * z^c1 * exp(-(z/a)^b)
	|z|<bailout
}

gumbel_1_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	z = pixel, c = p3*a*b:
	z = a*(z-o)
	z = c/exp(b/exp(z)+z)
	|z|<bailout
}

gumbel_1_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	z = 100, c = pixel*a*b:
	z = a*(z-o)
	z = c/exp(b/exp(z)+z)
	|z|<bailout
}

gumbel_2_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = a+1
	z = pixel, c = p3*a*b:
	z = z-o
	z = c/(z^c1*exp(b/z^a))
	|z|<bailout
}

gumbel_2_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	a = real(p1)
	b = imag(p1)
	o = real(p2)
	c1 = a+1
	z = 100, c = pixel*a*b:
	z = z-o
	z = c/(z^c1*exp(b/z^a))
	|z|<bailout
}

geometric_distj {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	p = real(p1)
	o = imag(p1)
	c1 = 1-p
	z = pixel, c = p3*p:
	z = c*c1^(z-o-1)
	|z|<bailout
}

geometric_distm (XAxis) {
	if(real(p4)==0)
		bailout=100
	else
		bailout=real(p4)
	endif
	p = real(p1)
	o = imag(p1)
	c1 = 1-p
	z = o+1, c = pixel*p:
	z = c*c1^(z-o-1)
	|z|<bailout
}