From: Rich Shepard <rshepard@ap...>  20071122 23:16:31

I see that I've been immortalized on the SciPy MatPlotLib Cookbook web page for my enquiry on plotting S and Zcurves. The Boltzman function serves very well for that purpose, and I've tweaked the example code to allow me to pass in the two endpoints and the midpoint for each of these curves. Now I need to plot normal curves (a.k.a. Gaussian or bell curves, depending on the background of the speaker/writer). I see that SciPy has a class for the normal curve in its stats package, and that the curve shape is defined by the mean and standard deviation. My need is to draw these curves based on the midpoint (== mean) and tail endpoints (which are not the same as the s.d.). I can copy  or call  the SciPy class from my application code, but I don't know if this is the most parsimonious approach. Your thoughts are appreciated. These plots will be used in two different media, and it two forms. Within the wxPython application I want to display each curve on one panel, and the set of curves related to one linguistic variable on a second panel. Then I need to have the set incorporated into a ReportLab report in .pdf. Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Angus McMorland <amcmorl@gm...>  20071123 00:08:21

On 23/11/2007, Rich Shepard <rshepard@...> wrote: > Now I need to plot normal curves (a.k.a. Gaussian or bell curves, > depending on the background of the speaker/writer). I see that SciPy has a > class for the normal curve in its stats package, and that the curve shape is > defined by the mean and standard deviation. 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 (fullwidth at halfmax), and c (centre) at positions given by r ''' return exp( (rc)**2 / fwhm2k( fwhm )**2 ) (released to public domain) > My need is to draw these curves based on the midpoint (== mean) and tail > endpoints (which are not the same as the s.d.). The midpoint here is c. It's not clear what you mean by endpoints  if you mean you want to be able to specify the y value at a given x deltax away from c, then it should be relatively simple to solve the equation to find the required fullwidth at halfmax to achieve these endpoints. After a very quick look (i.e. definitely needs verification), I think k = sqrt( (Rc)**2/log(Y) ) where Y is the desired value at distance Rc from the centre. >Your thoughts are appreciated. I hope that's what you're after. Angus,  AJC McMorland, PhD Student Physiology, University of Auckland 
From: Rich Shepard <rshepard@ap...>  20071123 04:16:54

On Fri, 23 Nov 2007, Angus McMorland wrote: > 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 (fullwidth at halfmax), > and c (centre) at positions given by r > ''' > return exp( (rc)**2 / fwhm2k( fwhm )**2 ) Thank you, Angus. I'll look at the Gaussian explanation to understand the input values. > The midpoint here is c. OK. > It's not clear what you mean by endpoints  if you mean you want to be > able to specify the y value at a given x deltax away from c, then it > should be relatively simple to solve the equation to find the required > fullwidth at halfmax to achieve these endpoints. After a very quick > look (i.e. definitely needs verification), I think What I mean is the x value where the tails of the curve have y == 0.0. These curves are defined by the range of x over which they are valid, and assume the midpoint is where y == 1.0 (the maximum value). The inflection points are at y = 0.5; in rare situations that may change. > I hope that's what you're after. I'll look at it in detail tomorrow (my time) and the weekend. I, too, hope that it's what I need. Much appreciated, Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Rich Shepard <rshepard@ap...>  20071124 00:01:02

On Thu, 22 Nov 2007, Rich Shepard wrote: >> 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 (fullwidth at halfmax), >> and c (centre) at positions given by r >> ''' >> return exp( (rc)**2 / fwhm2k( fwhm )**2 ) > > Thank you, Angus. I'll look at the Gaussian explanation to understand the > input values. It's been decades since I've needed to write any statistical or distributional function code; I used SysStat in DOS for quite some time and R for the past decade with linux. There's a great difference between using working code and writing one's one code. :) Looking more carefully at the SciPy stats module, I see the normal distribution is included as a class. However, I don't understand all the functions in that class  other than normal.pdf(x)  and do not see where loc and scale (representing mean and std. deviation) are used. Toward understanding how to use this code to draw the curves I need, I extracted the one method into a standalone file and tried to plot a Normal/Bell curve from 0 to 100: import matplotlib.numerix as nx import pylab as p from math import * def normal(x): return exp(x**2.0/2.0)/sqrt(2.0*pi) x = nx.arange(0, 100, 0.1) N = normal(x) p.plot(x, N, color='red', lw=2) p.axis([0, 100, 0.0, 1.0]) p.xlabel('Universe of Discourse') p.ylabel('Membership Grade') p.show() The error returned is: Traceback (most recent call last): File "normalcurve.py", line 17, in ? N = normal(x) File "normalcurve.py", line 14, in normal return exp(x**2.0/2.0)/sqrt(2.0*pi) TypeError: only length1 arrays can be converted to Python scalars A clue stick to the meaning of this error message, and how to fix the problem is needed. Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Jeff Whitaker <jswhit@fa...>  20071123 13:53:17

Rich Shepard wrote: > On Fri, 23 Nov 2007, Angus McMorland wrote: > > >> 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 (fullwidth at halfmax), >> and c (centre) at positions given by r >> ''' >> return exp( (rc)**2 / fwhm2k( fwhm )**2 ) >> > > Thank you, Angus. I'll look at the Gaussian explanation to understand the > input values. > > >> The midpoint here is c. >> > > OK. > > >> It's not clear what you mean by endpoints  if you mean you want to be >> able to specify the y value at a given x deltax away from c, then it >> should be relatively simple to solve the equation to find the required >> fullwidth at halfmax to achieve these endpoints. After a very quick >> look (i.e. definitely needs verification), I think >> > > What I mean is the x value where the tails of the curve have y == 0.0. > These curves are defined by the range of x over which they are valid, and > assume the midpoint is where y == 1.0 (the maximum value). The inflection > points are at y = 0.5; in rare situations that may change. > Rich: The tails of a Gaussian never reach zero  they just asymptote to zero for large x. Jeff  Jeffrey S. Whitaker Phone : (303)4976313 NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 325 Broadway Boulder, CO, USA 803053328 
From: Rich Shepard <rshepard@ap...>  20071123 15:55:53

On Fri, 23 Nov 2007, Jeff Whitaker wrote: > Rich: The tails of a Gaussian never reach zero  they just asymptote to zero > for large x. Jeff, For all practical purposes, that's fine. Usually any y value > 0.20 (the default) is considered functionally equivalent to zero. If the display looks close enough to the x axis, that fulfills that purpose. Because these curves are used to determine degree of set membership the tails are of no practical value. Thanks, Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Angus McMorland <amcmorl@gm...>  20071124 01:26:01

On 24/11/2007, Rich Shepard <rshepard@...> wrote: > On Fri, 23 Nov 2007, Angus McMorland wrote: > > > 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 (fullwidth at halfmax), > > and c (centre) at positions given by r > > ''' > > return exp( (rc)**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 fullwidth 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.  AJC McMorland, PhD Student Physiology, University of Auckland 
From: Rich Shepard <rshepard@ap...>  20071124 04:17:12

On Sat, 24 Nov 2007, Angus McMorland wrote: > fwhm is the fullwidth 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  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Angus McMorland <amcmorl@gm...>  20071124 04:29:58

On 24/11/2007, Rich Shepard <rshepard@...> wrote: > On Sat, 24 Nov 2007, Angus McMorland wrote: > > > fwhm is the fullwidth 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.  AJC McMorland, PhD Student Physiology, University of Auckland 
From: Rich Shepard <rshepard@ap...>  20071124 15:01:49

On Sat, 24 Nov 2007, Angus McMorland wrote: > 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 SpringerVerlag. Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Rich Shepard <rshepard@ap...>  20071124 15:41:17

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 matplotlib0.87, but which should be fixed in 0.90.1: Traceback (most recent call last): File "normalcurve.py", line 19, in ? G = gauss1d(x, fwhm, center) File "normalcurve.py", line 16, in gauss1d return exp((rcenter)**2 / fwhm2k(fwhm)**2) TypeError: only length1 arrays can be converted to Python scalars According to /usr/lib/python2.4/sitepackages/numpy/version.py, the installed version is 1.0b5. Do I have a version incompatibility here? Rich  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Angus McMorland <amcmorl@gm...>  20071124 23:08:32

On 25/11/2007, Rich Shepard <rshepard@...> 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 matplotlib0.87, but which should be fixed in 0.90.1: > > Traceback (most recent call last): > File "normalcurve.py", line 19, in ? > G = gauss1d(x, fwhm, center) > File "normalcurve.py", line 16, in gauss1d > return exp((rcenter)**2 / fwhm2k(fwhm)**2) > TypeError: only length1 arrays can be converted to Python scalars > > According to /usr/lib/python2.4/sitepackages/numpy/version.py, the > installed version is 1.0b5. > > Do I have a version incompatibility here? 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 (normalcurve.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.  AJC McMorland, PhD Student Physiology, University of Auckland 
From: Rich Shepard <rshepard@ap...>  20071124 23:18:05

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 > (normalcurve.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 (fullwidth at halfmax), and c (centre) at positions given by r ''' return exp((rcenter)**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  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Fernando Perez <fperez.net@gm...>  20071125 00:25:27

On Nov 24, 2007 4:17 PM, Rich Shepard <rshepard@...> 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. 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 
From: Angus McMorland <amcmorl@gm...>  20071125 01:12:34

On 25/11/2007, Rich Shepard <rshepard@...> 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 > > (normalcurve.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 (fullwidth at halfmax), > and c (centre) at positions given by r > ''' > return exp((rcenter)**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() > 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) arrayfriendly versions. If you change exp > nx.exp in your definition of gauss1d, all works okay. A.  AJC McMorland, PhD Student Physiology, University of Auckland 
From: Rich Shepard <rshepard@ap...>  20071125 02:04:19

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) arrayfriendly versions. If you change exp > > nx.exp in your definition of gauss1d, all works okay. 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  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 
From: Rich Shepard <rshepard@ap...>  20071125 02:17:21

On Sun, 25 Nov 2007, Angus McMorland wrote: > If you change exp > nx.exp in your definition of gauss1d, all works okay. 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  Richard B. Shepard, Ph.D.  Integrity Credibility Applied Ecosystem Services, Inc.  Innovation <http://www.applecosys.com>; Voice: 5036674517 Fax: 5036678863 