matplotlib-users

[Matplotlib-users] [matplotlib-users] How to plot digamma function (psi)

[Fabien L. <laf...@gm...>]

Hello everyone,

I try to plot the digamma function of (1/2 + 1/x) but I'm not sure that I'm
plotting the good one.

I've tried:

special.polygamma(0, (1/2 + 1/x))

and

special.polygamma(1, (1/2 + 1/x))

but I don't have the same result as with mathcad.

I've tried to code it like that:

def F(x): return  mpmath.diff(lambda x: gamma(1/2 + 1/x),1)/gamma(1/2 + 1/x)

But It returns zero division error even when x is in ]0,1]

Any idea?

Re: [Matplotlib-users] [matplotlib-users] How to plot digamma function (psi)

[Damon M. <dam...@gm...>]

On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:
> Hello everyone,
> 
> I try to plot the digamma function of (1/2 + 1/x) but I'm not sure that I'm
> plotting the good one.
> 
> I've tried:
> 
> special.polygamma(0, (1/2 + 1/x))
> 
> and
> 
> special.polygamma(1, (1/2 + 1/x))

You want special.polygamma(0, (1/2 + 1/x)). See
http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.polygamma.html

The number specifies which derivative of the digamma function you want.
Surely you want the 0th derivative?

> But It returns zero division error even when x is in ]0,1]

I think it blows up at x = 0. What is the type of x in your usecase? Is
it an array? If x contains the element 0, you will get a zero
division error. You could try plotting the points explicitly:

from numpy import linspace
from pylab import *

x = linspace(0.5, 2, num=100, endpoint=True)
y = special.polygamma(0, (1/2 + 1/x))
plot(x, y)
show()

You can compare output against this:
http://www.wolframalpha.com/input/?i=digamma%281%2F2+%2B+1%2Fx%29+between+0.5+and+2

Hope this helps.

-- 
Damon McDougall
http://damon-is-a-geek.com
B2.39
Mathematics Institute
University of Warwick
Coventry
West Midlands
CV4 7AL
United Kingdom

Re: [Matplotlib-users] [matplotlib-users] How to plot digamma function (psi)

[Benjamin R. <ben...@ou...>]

On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall
<dam...@gm...>wrote:

> On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:
> > Hello everyone,
> >
> > I try to plot the digamma function of (1/2 + 1/x) but I'm not sure that
> I'm
> > plotting the good one.
> >
> > I've tried:
> >
> > special.polygamma(0, (1/2 + 1/x))
> >
> > and
> >
> > special.polygamma(1, (1/2 + 1/x))
>
> You want special.polygamma(0, (1/2 + 1/x)). See
>
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.polygamma.html
>
> The number specifies which derivative of the digamma function you want.
> Surely you want the 0th derivative?
>
> > But It returns zero division error even when x is in ]0,1]
>
> I think it blows up at x = 0. What is the type of x in your usecase? Is
> it an array? If x contains the element 0, you will get a zero
> division error. You could try plotting the points explicitly:
>
> from numpy import linspace
> from pylab import *
>
> x = linspace(0.5, 2, num=100, endpoint=True)
> y = special.polygamma(0, (1/2 + 1/x))
> plot(x, y)
> show()
>
> You can compare output against this:
>
> http://www.wolframalpha.com/input/?i=digamma%281%2F2+%2B+1%2Fx%29+between+0.5+and+2
>
> Hope this helps.
>
>
Another problem might be the "1/2" part, which in python2.x would yield 0
unless one does "from __future__ import division".

Ben Root

Re: [Matplotlib-users] [matplotlib-users] How to plot digamma function (psi)

[Damon M. <dam...@gm...>]

On Tue, Jul 10, 2012 at 08:57:24AM -0400, Benjamin Root wrote:
> On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall
> <dam...@gm...>wrote:
> 
> > On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:
> >
> > > But It returns zero division error even when x is in ]0,1]
> >
> > I think it blows up at x = 0. What is the type of x in your usecase? Is
> > it an array? If x contains the element 0, you will get a zero
> > division error. You could try plotting the points explicitly:
> >
> Another problem might be the "1/2" part, which in python2.x would yield 0
> unless one does "from __future__ import division".
> 
> Ben Root

Wow, I can't believe I didn't spot that. Nice one.

I will update my answer according to Ben's astute observation:

from scipy import special
from pylab import *

x = linspace(0.5, 2.0, num=100, endpoint=True)
y = special.polygamma(0, 0.5 + 1.0/x)
plot(x, y)
show()

Thanks Ben.

-- 
Damon McDougall
http://damon-is-a-geek.com
B2.39
Mathematics Institute
University of Warwick
Coventry
West Midlands
CV4 7AL
United Kingdom

Re: [Matplotlib-users] [matplotlib-users] How to plot digamma function (psi)

[Fabien L. <laf...@gm...>]

Thanks! The problem came from the 1/2 ! For convenience I've found the
function digamma on numpy

*http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.psi.html

But it's quite hard to find it! Maybe we can ask to add digamma in the
title between parenthesis?


Fabien
*
2012/7/10 Damon McDougall <dam...@gm...>

> On Tue, Jul 10, 2012 at 08:57:24AM -0400, Benjamin Root wrote:
> > On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall
> > <dam...@gm...>wrote:
> >
> > > On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:
> > >
> > > > But It returns zero division error even when x is in ]0,1]
> > >
> > > I think it blows up at x = 0. What is the type of x in your usecase? Is
> > > it an array? If x contains the element 0, you will get a zero
> > > division error. You could try plotting the points explicitly:
> > >
> > Another problem might be the "1/2" part, which in python2.x would yield 0
> > unless one does "from __future__ import division".
> >
> > Ben Root
>
> Wow, I can't believe I didn't spot that. Nice one.
>
> I will update my answer according to Ben's astute observation:
>
> from scipy import special
> from pylab import *
>
> x = linspace(0.5, 2.0, num=100, endpoint=True)
> y = special.polygamma(0, 0.5 + 1.0/x)
> plot(x, y)
> show()
>
> Thanks Ben.
>
> --
> Damon McDougall
> http://damon-is-a-geek.com
> B2.39
> Mathematics Institute
> University of Warwick
> Coventry
> West Midlands
> CV4 7AL
> United Kingdom
>