Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.graphics.apps.gnuplot > #2030 > unrolled thread

plotting a nice random mountain on a plain

Back to article view | Back to comp.graphics.apps.gnuplot

Contents

  plotting a nice random mountain on a plain Karl <mail.kfr@gmx.net> - 2013-09-12 14:38 +0200
    Re: plotting a nice random mountain on a plain sfeam <sfeam@users.sourceforge.net> - 2013-09-12 20:38 -0700

#2030 — plotting a nice random mountain on a plain

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?

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

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))

[toc] | [next] | [standalone]

#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))

[toc] | [prev] | [standalone]

Back to top | Article view | comp.graphics.apps.gnuplot

csiph-web