Path: csiph.com!newsfeed.hal-mli.net!feeder3.hal-mli.net!newsfeed.hal-mli.net!feeder1.hal-mli.net!feeder.erje.net!us.feeder.erje.net!newsfeed.fsmpi.rwth-aachen.de!eternal-september.org!feeder.eternal-september.org!news.eternal-september.org!.POSTED!not-for-mail
From: sfeam <sfeam@users.sourceforge.net>
Newsgroups: comp.graphics.apps.gnuplot
Subject: Re: plotting a nice random mountain on a plain
Followup-To: comp.graphics.apps.gnuplot
Date: Thu, 12 Sep 2013 20:38:20 -0700
Organization: gnuplot development team
Lines: 77
Message-ID: <l0u1bf$9r4$1@dont-email.me>
References: <l0sckn$6dq$1@news.rz.uni-karlsruhe.de>
Reply-To: sfeam@users.sourceforge.net
Mime-Version: 1.0
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 7Bit
Injection-Date: Fri, 13 Sep 2013 03:38:23 +0000 (UTC)
Injection-Info: mx05.eternal-september.org; posting-host="e947cb29d2533e27aa00183e80187dbb"; logging-data="10084"; mail-complaints-to="abuse@eternal-september.org";	posting-account="U2FsdGVkX1/MUXz5ZwTOmsdxCnM85czB"
User-Agent: KNode/4.4.9
Cancel-Lock: sha1:+1W6n7Ifmp1LS1VQ3Fjw1syjiNg=
Xref: csiph.com comp.graphics.apps.gnuplot:2031

Karl wrote:

> Hi,
> wrote the thing below yesterday, could be something for the demo set?
> Otherwise, just enjoy!
> 
> On thing puzzles me, though: the first plot just draws a straight line
> at y=1 in logscale y mode (commented out). Anyone know why?

The "smooth" routines don't handle log-scale.  See for example
https://sourceforge.net/p/gnuplot/bugs/1263/
It's not a simple coding error, but rather a total failure to consider
the possibility of log scaling when the code was written.


> And how is the "invnorm()" function done in gp? I thought there was no
> analytical form for it.

It uses a 2-piece polynomial approximation.  The terms of the polynomial
are in the source file specfun.c.   I didn't see a corresponding citation
(quick look only) but many of the polynomial approximations are taken from
Abramowitz & Stegun so I would look there first to see if the terms match.

The development version of gnuplot also offers the option to link 
to the external math library libcephes.  That has a different
implementation of inverse normal, but I think its accuracy is about
the same as the one in specfun.  I'd be interested in a more serious
investigation into their respective accuracy if someone wants to make
the effort to test.

	Ethan

> 
> Karl
> 
> -----------
> 
> #!/gnuplot
> # 2D plot of a mountain constructed over the density of normally
> # distributed random points on a plain
> 
> # Box-Muller method to create (pairs of) normally distributed
> # random numbers from a pair of [0:1] evenly distributed
> # random numbers
> # This should give the same as "invnorm(rand(0))"
> 
> normd(x) = (u1=rand(0),u2=rand(0), sqrt(-2*log(u2)) * cos(2*pi*u1))
> # or
> #normd(x) = (u1=rand(0),u2=rand(0), sqrt(-2*log(u2)) * sin(2*pi*u1))
> 
> #set logsc y
> set samp 1000
> plot "+" us (normd(0)):(1) smooth kdens title "Box-Muller", \
>      "+" us (invnorm(rand(0))):(1) smooth kdens title "gp46 invnorm"
> #unset logsc y
> pause -1
> 
> # normdc(x) returns normally distributed random numbers
> # if x is complex, the random number is complex too,
> # with the real and imaginary part being the pair of
> # numbers the method of Box and Muller returns
> 
> normdc(x) = (u1=rand(0),u2=rand(0), sqrt(-2*log(u2)) * \
>   (cos(2*pi*u1) + (imag(x)==0?0:(sqrt(-1)*sin(2*pi*u1)))))
> 
> 
> set xrange [-4:4]
> set yrange [-4:4]
> # less samples make the mountain more interesting ;-)
> set samp 1000
> # so does changing the parameters of kdensity2d
> set dgrid3d 80,80 gauss kdensity2d 0.3,0.3
> set hidden
> #set logsc z
> splot "+" us (real(p=normdc(sqrt(-1)))):(imag(p)):(1) w l
> 
> #stats "+" us (real(p=normdc(sqrt(-1)))):(imag(p))