Plotting Continuous Functions

Community matplotlib-users
1

> For parsimony, I think you're probably best off just using the Gaussian
> equation:
>
> def fwhm2k(fwhm):
> '''converts fwhm value to k (see above)'''
> return fwhm/(2 * n.sqrt( n.log( 2 ) ) )
>
> def gauss1d(r, fwhm, c):
> '''returns the 1d gaussian given by fwhm (full-width at half-max),
> and c (centre) at positions given by r
> '''
> return exp( -(r-c)**2 / fwhm2k( fwhm )**2 )
>
> (released to public domain)

Angus,

   I'm trying to find the context for the above so I know what to feed fwhm2k
as the fwhm value.

fwhm is the full-width at half the maximum height, i.e. it's the
difference between the two values of x when:

r - c| = 0.5

The fwhm is a shape parameter (like std dev) - it determines the width
of the curve. The combination of width and the range of values you
plot (r) determine how close the function gets to zero, and how much
of it is plotted. As Jeff said, it'll never actually reach zero, so
you have to decide how close is close enough.

You don't need to call fwhm2k yourself; it's called by the gauss1d
function. I just do it that way because the equation uses k, but I'm
always interested in fwhm.

   It's been decades since I needed to work with continuous distributions and
my insights and skills have rusted.

Perhaps the easiest thing is to shove it into some quick code and play
around with the values so you see how it works:

import pylab as p, numpy as n
x = n.arange(100) - 50.
fwhm = 25
centre = 0
y = gauss1d(x, fwhm, centre)
p.plot(x,y)

If you have other questions, you'll need to be a bit more specific so
we can address them directly.

Angus.

···

On 24/11/2007, Rich Shepard <rshepard@...695...> wrote:

On Fri, 23 Nov 2007, Angus McMorland wrote:

--
AJC McMorland, PhD Student
Physiology, University of Auckland

2

fwhm is the full-width at half the maximum height, i.e. it's the
difference between the two values of x when:

>r - c| = 0.5

Angus,

   The additional explanation helps a lot.

The fwhm is a shape parameter (like std dev) - it determines the width of
the curve. The combination of width and the range of values you plot (r)
determine how close the function gets to zero, and how much of it is
plotted. As Jeff said, it'll never actually reach zero, so you have to
decide how close is close enough.

You don't need to call fwhm2k yourself; it's called by the gauss1d
function. I just do it that way because the equation uses k, but I'm
always interested in fwhm.

   Ah, I missed seeing that.

Perhaps the easiest thing is to shove it into some quick code and play
around with the values so you see how it works:

   That's what I intend to do. I've been running ipython and also writing
small scripts to understand how the functions work ... and plot. Playing
with parameters clarifies everything.

Much appreciated,

Rich

···

On Sat, 24 Nov 2007, Angus McMorland wrote:

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

3

> fwhm is the full-width at half the maximum height, i.e. it's the
> difference between the two values of x when:
>
> >r - c| = 0.5

Looking at my reply, I realised this was rubbish - sorry about that.
The fwhm is the difference between the two values of x that give Y =
0.5.

   The additional explanation helps a lot.

Great. Hopefully this correction will make things even more clear.

Angus.

···

On 24/11/2007, Rich Shepard <rshepard@...695...> wrote:

On Sat, 24 Nov 2007, Angus McMorland wrote:

--
AJC McMorland, PhD Student
Physiology, University of Auckland

4

Looking at my reply, I realised this was rubbish - sorry about that. The
fwhm is the difference between the two values of x that give Y = 0.5.

   Now that makes much more sense. Having control over the x values for the
inflection point allows us to tune the curve shape to more precisely reflect
the underlying semantics.

   The additional explanation helps a lot.

Great. Hopefully this correction will make things even more clear.

   I had scanned your previous response without looking deeply at it because
I was not going to work on this until today in any case.

   For anyone curious about the context for my questions, it is an
approximate reasoning (i.e., expert system) model using fuzzy sets and fuzzy
logic. You can read all about it in my book, "Quantifying Environmental
Impact Assessments Using Fuzzy Logic," published in 2005 by Springer-Verlag.

Rich

···

On Sat, 24 Nov 2007, Angus McMorland wrote:

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

5

While the functions and equations are now clear, I get an error that was
present in matplotlib-0.87, but which should be fixed in -0.90.1:

Traceback (most recent call last):
   File "normal-curve.py", line 19, in ?
     G = gauss1d(x, fwhm, center)
   File "normal-curve.py", line 16, in gauss1d
     return exp(-(r-center)**2 / fwhm2k(fwhm)**2)
TypeError: only length-1 arrays can be converted to Python scalars

   According to /usr/lib/python2.4/site-packages/numpy/version.py, the
installed version is 1.0b5.

   Do I have a version incompatibility here?

Rich

···

On Sat, 24 Nov 2007, Angus McMorland wrote:

Great. Hopefully this correction will make things even more clear.

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

6

I'm not completely sure, but I suspect that this is an implementation
bug, rather than a version bug, particularly because the line in
question isn't involving matplotlib at all. If you post the relevant
code (normal-curve.py, by the looks of things), it might be easy to
spot the problem.

I've found it easiest to solve these sorts of bugs by running the code
in an ipython shell, with automatic pdb calling. That way you can
inspect the values of the parameters in question - one of which is, I
think, the problem here.

Angus.

···

On 25/11/2007, Rich Shepard <rshepard@...695...> wrote:

On Sat, 24 Nov 2007, Angus McMorland wrote:

> Great. Hopefully this correction will make things even more clear.

   While the functions and equations are now clear, I get an error that was
present in matplotlib-0.87, but which should be fixed in -0.90.1:

Traceback (most recent call last):
   File "normal-curve.py", line 19, in ?
     G = gauss1d(x, fwhm, center)
   File "normal-curve.py", line 16, in gauss1d
     return exp(-(r-center)**2 / fwhm2k(fwhm)**2)
TypeError: only length-1 arrays can be converted to Python scalars

   According to /usr/lib/python2.4/site-packages/numpy/version.py, the
installed version is 1.0b5.

   Do I have a version incompatibility here?

--
AJC McMorland, PhD Student
Physiology, University of Auckland

7

I'm not completely sure, but I suspect that this is an implementation bug,
rather than a version bug, particularly because the line in question isn't
involving matplotlib at all. If you post the relevant code
(normal-curve.py, by the looks of things), it might be easy to spot the
problem.

Angus,

   I've seen the same error trying to plot other curves in the past couple of
days, but not those using the Boltzmann distribution. Here's the file:

import matplotlib.numerix as nx
import pylab as p
from math import *

center = 50.0
fwhm = 50.0

def fwhm2k(fwhm):
   '''converts fwhm value to k (see above)'''
   return fwhm/(2 * nx.sqrt(nx.log(2)))

def gauss1d(r, fwhm, center):
   '''returns the 1d gaussian given by fwhm (full-width at half-max),
   and c (centre) at positions given by r
   '''
   return exp(-(r-center)**2 / fwhm2k(fwhm)**2)

x = nx.arange(0, 100, 0.1)
G = gauss1d(x, fwhm, center)
p.plot(x, G, color='red', lw=2)
p.axis([0, 100, 0.0, 1.0])
p.xlabel('Universe of Discourse')
p.ylabel('Membership Grade')
p.show()

I've found it easiest to solve these sorts of bugs by running the code in
an ipython shell, with automatic pdb calling. That way you can inspect the
values of the parameters in question - one of which is, I think, the
problem here.

   I've not run ipython with pdb; I'll look at the docs to learn how. I do
use winpdb on the application.

Thanks,

Rich

···

On Sun, 25 Nov 2007, Angus McMorland wrote:

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

8

If you type %pdb *before* running your scripts, then any exception
that fires will automatically activate pdb.

But for a while we've had a more convenient way to access pdb, which
is the new %debug command. At any time if you simply type %debug, the
pdb debugger will activate into the last unhandled exception. So as
long as you don't wait too long after seeing an exception (since the
system only works with the *last* one, if you get a new exception from
a typo at the command line you lose the chance to inspect your
program), you can use it in a more fluid way than letting %pdb
forcefully activate every single time.

Cheers,

f

···

On Nov 24, 2007 4:17 PM, Rich Shepard <rshepard@...695...> wrote:

On Sun, 25 Nov 2007, Angus McMorland wrote:

> I've found it easiest to solve these sorts of bugs by running the code in
> an ipython shell, with automatic pdb calling. That way you can inspect the
> values of the parameters in question - one of which is, I think, the
> problem here.

   I've not run ipython with pdb; I'll look at the docs to learn how. I do
use winpdb on the application.

9

As I suspected, this is a parameter issue- in this case caused by your
use of the ath module routines which require scalar input, rather than
numpy's (or matplotlib's numerix's) array-friendly versions. If you
change exp -> nx.exp in your definition of gauss1d, all works okay.

A.

···

On 25/11/2007, Rich Shepard <rshepard@...695...> wrote:

On Sun, 25 Nov 2007, Angus McMorland wrote:

> I'm not completely sure, but I suspect that this is an implementation bug,
> rather than a version bug, particularly because the line in question isn't
> involving matplotlib at all. If you post the relevant code
> (normal-curve.py, by the looks of things), it might be easy to spot the
> problem.

Angus,

   I've seen the same error trying to plot other curves in the past couple of
days, but not those using the Boltzmann distribution. Here's the file:

import matplotlib.numerix as nx
import pylab as p
from math import *

center = 50.0
fwhm = 50.0

def fwhm2k(fwhm):
   '''converts fwhm value to k (see above)'''
   return fwhm/(2 * nx.sqrt(nx.log(2)))

def gauss1d(r, fwhm, center):
   '''returns the 1d gaussian given by fwhm (full-width at half-max),
   and c (centre) at positions given by r
   '''
   return exp(-(r-center)**2 / fwhm2k(fwhm)**2)

x = nx.arange(0, 100, 0.1)
G = gauss1d(x, fwhm, center)
p.plot(x, G, color='red', lw=2)
p.axis([0, 100, 0.0, 1.0])
p.xlabel('Universe of Discourse')
p.ylabel('Membership Grade')
p.show()

--
AJC McMorland, PhD Student
Physiology, University of Auckland

10

Oh, d'oh! That's what I get for copying without closely examining the
code. With more experience I'll learn the differences between the Python
math functions and the equivalent ones from NumPy.

Many thanks,

Rich

···

On Sun, 25 Nov 2007, Angus McMorland wrote:

As I suspected, this is a parameter issue- in this case caused by your use
of the ath module routines which require scalar input, rather than numpy's
(or matplotlib's numerix's) array-friendly versions. If you change exp ->
nx.exp in your definition of gauss1d, all works okay.

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863

11

Angus,

    Yes, it works just fine. By adjusting the value of the fwhm parameter I
can produce the curves we need for both display and printing.

    Now I can spend some time tomorrow working out the functions for
trapezoidal and shouldered curves (comparatively quite simple), learning how
to embed them in a wxPython panel, and putting two, three, five, or seven
curves on the same axes. The embedded wiki page and User Guide should
provide what I need.

    My thanks to you and everyone else who responded. While I don't yet have
sufficient knowledge to help others on this mail list, I do provide help on
those lists where I've used the software for months or years and can answer
questions posed by others. Perhaps I'll be doing that here, too, soon.

Carpe weekend -- what's left of it,

Rich

···

On Sun, 25 Nov 2007, Angus McMorland wrote:

If you change exp -> nx.exp in your definition of gauss1d, all works okay.

--
Richard B. Shepard, Ph.D. | Integrity Credibility
Applied Ecosystem Services, Inc. | Innovation
<http://www.appl-ecosys.com> Voice: 503-667-4517 Fax: 503-667-8863