On Mon, 2005-03-07 at 12:04 -0800, Stephen Walton wrote:
Robert Kern wrote:
If it's univariate, and you can write out the pdf or cdf as a function, then I believe you can subclass scipy.stats.rv_continuous, and it's rvs() method will numerically invert the cdf to generate it's random numbers.
This is so cool! I had a desire to generate values on (0,1) where values near 0.5 were less probable than values at the endpoints. Here's the implementation:
#---------------begin------------------------------ import Numeric as Num from scipy.stats.distributions import rv_continuous
# # CDF for the sunspot generation function. If we make it a subclass # of rv_continuous we get sunspot.rvs for free :-) #
class sunspot_gen(rv_continuous): pmin=0.1 def _pdf(self,x): pmax=(Num.pi-2*self.pmin)/(Num.pi-2) return((pmax-self.pmin)*(1-Num.sin(Num.pi*x))+self.pmin) def _cdf(self,x): pmax=(Num.pi-2*self.pmin)/(Num.pi-2) return -(cos(pi*x)*pmin-cos(pi*x)*pmax-pmax*x*pi-pmin+pmax)/pi sunspot = sunspot_gen(a=0.,b=1.,name='sunspot', longname='A sunspot subdivision',extradoc="""
Sunspot distribution
This distribution represents the probability of a subdivision of a sunspot into two spots of size x and (1-x), where x is a value from sunspot.rvs(0). The PDF has a high probability of x near 0 or 1 and a low probability of x near 0.5.
""") #-----------------end--------------------------------
My only question: what should I replace the 'import Numeric as Num' with if I want to be able to work within the framework of using either Numeric or numarray? 'import numerix as Num' doesn't seem to work.
Try "import scipy_base.numerix as Num" but understand that it probably won't work with numarray until scipy as a whole is ported (right now all we've got is scipy_base). It should work with Numeric however. Regards, Todd