Probability distribution functions and the Triangle Map
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 }
Dear reader, hallo Morgan I took a short look to you frm's about distributions. well that's to say only to the frm of the normal distribution. I hope you don't mind I changed it a bit. gaussian_distm (XAxis) {;standard distribution take real(p2)=1, the other variables 0 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 } and not gaussian_distm2 (XAxis) {;standard distribution take real(p2)=1, the other variables 0.Real (p2) must never be 0 if(real(p4)==0) bailout=100 else bailout=real(p4) endif sigma =p2 mu =p1 z = 0, c =1/(sqr(2*pi)*sigma),scalingfactor=50, c1=scalingfactor*pixel: z = c*exp(-.5*((z-mu)/sigma)^2)*c1 |z|<bailout } notes: 1.The 50 is only there because then you get when using p1=(0,0) and sigma=(1,0), that's to say the standard normal distribution, directly something on the screen. If you don't put it in then you have only a part of the basin of attraction. 2.I don't know why you omit the fasctor .5 in the exonent of e. 3. I don't see any reason not to allow complex numbers for mean and variance. 4. I don't know ehy first squaring sigma and then take the squareroot of sigma Cheers, Jos Hendriks
jos hendriks wrote:
Dear reader, hallo Morgan
notes: 1.The 50 is only there because then you get when using p1=(0,0) and sigma=(1,0), that's to say the standard normal distribution, directly something on the screen. If you don't put it in then you have only a part of the basin of attraction. 2.I don't know why you omit the fasctor .5 in the exonent of e. 3. I don't see any reason not to allow complex numbers for mean and variance. 4. I don't know ehy first squaring sigma and then take the squareroot of sigma
1 and 3 are pretty much answered by the original post. I didn't mess around with the definition much: 1 (-(x-mu)^2/(2 sigma^2)) p(x) = ------------------ e sqrt(2 pi sigma^2) Where sigma is the standard deviation, and mu the mean (which I called o for consistency with the others). My biggest fiddle was to chuck in a c parameter: c (-(x-o)^2/(2 sigma^2)) p(x) = ------------------ e sqrt(2 pi sigma^2) which in the julia sets is an arbitrary complex constant, and varies over the plane for the Mandelbrot. Short answer to 2 and 4: For 2, the factor is there if you look for it. For 4, the answer is because it's cheaper to compute that way. Longer answer to 2 and 4: Since c ------------------ sqrt(2 pi sigma^2) is constant, and the only place where c appears, I set c equal to that. c = {p3|pixel}/sqrt(2*pi*sqr(sigma)): z = c*exp(-sqr(z-o)/(2*sqr(sigma))) 2*sqr(sigma) appears twice, and it's constant, so it's (a) a waste of time computing it twice, and (b) a waste of time computing it inside the loop. So, pull it out into another constant, c1: c1 =2*sqr(sigma) c = {p3|pixel}/sqrt(pi*c1): z = c*exp(-sqr(z-o)/c1) The difference between sqrt(sqr(sigma)) and sigma is of course that for sigma real, sqrt(sqr(sigma)) is nonnegative. Using just sigma would (a) change the behaviour, and (b) mean an extra multplication in the initialiser, viz: c1 =2*sqr(sigma) c = {p3|pixel}/(sigma*sqrt(pi*2)): z = c*exp(-sqr(z-o)/c1) Morgan L. Owens "The bit that doesn't go away when you stop believing in it."
Hallo Morgan Thank you for the clear reaction Sorry I missed your factor 1/2. I prefer to use names like mu and sigma, because they are widely agreed upon. I took a look to the log normal distribution and found it interesting. Maybe you appreciate this: Log_Normal_-1_-1 { ; Jos Hendriks ,2002 reset=2000 type=formula formulafile=fractint.frm formulaname=lognormal_distm passes=1 center-mag=0.169752/-0.0312433/0.1751883 params=-1/-1/0/0 float=y colors=000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z0\ 0e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz\ 0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0z\ zzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00\ z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL000\ 0<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zz\ LzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e0\ 0L00zzLzL0000<2>00z<2>0zzzzezz0<2>z00e00L00zzLzL0zL0ccc<29>ccc } Much greetings
jos hendriks wrote:
Hallo Morgan
Thank you for the clear reaction Sorry I missed your factor 1/2.
I suffer from an embarrassing psychosomatic condition called "premature optimisation" ... that's why I showed an interest in Gerald K. Dobiasovsky's technique (A fellow sufferer!).
I prefer to use names like mu and sigma, because they are widely agreed upon.
So do I; I only changed "mu" to "o" in the gaussian distribution because there the two are equivalent, but the same isn't true in general. The main reason is because the effects of 'mu', 'zeta' and 'sigma' are already familiar to those familiar with the distribution in question. Maybe I should take your suggestion on - be even more relaxed than I already am and allow them to become complex... There are always all sorts of directions one can take, and one can take only a finite number of them ...
I took a look to the log normal distribution and found it interesting. Maybe you appreciate this:
See? Told ya! Midgets! I _wasn't_ imagining them! :) Had a look around the image; saw, among other things, a vive voce in progress. Didn't record that one, though; recorded this one instead: Taken_in_Lieu { ; A Critical Point in Negotiations ; ('snot your job) ; After Jos Hendriks, 2002 ; Version 2001 Patchlevel 8 reset=2001 type=formula formulafile=distro.frm formulaname=lognormal_distm passes=1 center-mag=0.1353/0/1.5 params=-1/-1/0/0/0/0/0/0 float=y outside=fmod colors=000_XR<43>_hj_hj_hk<2>_il_il`jk<60>tnYunYunX<2>vnXvnXumW<73>ZZJYZ\ JYZJ<3>XYIFN9WYI<52>FN9 } Seems to give the engine a bit of grief here and there - as tends to happen when dividing by (really small) z. Morgan L. Owens "_'Allow'?_"
participants (2)
-
jos hendriks -
Morgan L. Owens