Python scipy.stats.vonmises() Examples

def test_all_distributions():
    for dist in dists:
        distfunc = getattr(stats, dist)
        nargs = distfunc.numargs
        alpha = 0.01
        if dist == 'fatiguelife':
            alpha = 0.001

        if dist == 'frechet':
            args = tuple(2*rand(1))+(0,)+tuple(2*rand(2))
        elif dist == 'triang':
            args = tuple(rand(nargs))
        elif dist == 'reciprocal':
            vals = rand(nargs)
            vals[1] = vals[0] + 1.0
            args = tuple(vals)
        elif dist == 'vonmises':
            yield check_distribution, dist, (10,), alpha
            yield check_distribution, dist, (101,), alpha
            args = tuple(1.0+rand(nargs))
        else:
            args = tuple(1.0+rand(nargs))

        yield check_distribution, dist, args, alpha 

def check_vonmises_pdf_periodic(k,l,s,x):
    vm = stats.vonmises(k,loc=l,scale=s)
    assert_almost_equal(vm.pdf(x),vm.pdf(x % (2*numpy.pi*s))) 

def check_vonmises_cdf_periodic(k,l,s,x):
    vm = stats.vonmises(k,loc=l,scale=s)
    assert_almost_equal(vm.cdf(x) % 1,vm.cdf(x % (2*numpy.pi*s)) % 1) 

def cases_test_all_distributions():
    np.random.seed(1234)

    for dist in dists:
        distfunc = getattr(stats, dist)
        nargs = distfunc.numargs
        alpha = 0.01
        if dist == 'fatiguelife':
            alpha = 0.001

        if dist == 'trapz':
            args = tuple(np.sort(np.random.random(nargs)))
        elif dist == 'triang':
            args = tuple(np.random.random(nargs))
        elif dist == 'reciprocal' or dist == 'truncnorm':
            vals = np.random.random(nargs)
            vals[1] = vals[0] + 1.0
            args = tuple(vals)
        elif dist == 'vonmises':
            yield dist, (10,), alpha
            yield dist, (101,), alpha
            args = tuple(1.0 + np.random.random(nargs))
        else:
            args = tuple(1.0 + np.random.random(nargs))

        yield dist, args, alpha 

def check_vonmises_pdf_periodic(k, l, s, x):
    vm = stats.vonmises(k, loc=l, scale=s)
    assert_almost_equal(vm.pdf(x), vm.pdf(x % (2*numpy.pi*s))) 

def check_vonmises_cdf_periodic(k, l, s, x):
    vm = stats.vonmises(k, loc=l, scale=s)
    assert_almost_equal(vm.cdf(x) % 1, vm.cdf(x % (2*numpy.pi*s)) % 1) 

def test_vonmises_numerical():
    vm = stats.vonmises(800)
    assert_almost_equal(vm.cdf(0), 0.5) 

def _fitstart_beta(self, x, fixed=None):
    '''method of moment estimator as starting values for beta distribution

    Parameters
    ----------
    x : array
        data for which the parameters are estimated
    fixed : None or array_like
        sequence of numbers and np.nan to indicate fixed parameters and parameters
        to estimate

    Returns
    -------
    est : tuple
        preliminary estimates used as starting value for fitting, not
        necessarily a consistent estimator

    Notes
    -----
    This needs to be written and attached to each individual distribution

    References
    ----------
    for method of moment estimator for known loc and scale
    http://en.wikipedia.org/wiki/Beta_distribution#Parameter_estimation
    http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm
    NIST reference also includes reference to MLE in
    Johnson, Kotz, and Balakrishan, Volume II, pages 221-235

    '''
    #todo: separate out this part to be used for other compact support distributions
    #      e.g. rdist, vonmises, and truncnorm
    #      but this might not work because it might still be distribution specific
    a, b = x.min(), x.max()
    eps = (a-b)*0.01
    if fixed is None:
        #this part not checked with books
        loc = a - eps
        scale = (a - b) * (1 + 2*eps)
    else:
        if np.isnan(fixed[-2]):
            #estimate loc
            loc = a - eps
        else:
            loc = fixed[-2]
        if np.isnan(fixed[-1]):
            #estimate scale
            scale = (b + eps) - loc
        else:
            scale = fixed[-1]

    #method of moment for known loc scale:
    scale = float(scale)
    xtrans = (x - loc)/scale
    xm = xtrans.mean()
    xv = xtrans.var()
    tmp = (xm*(1-xm)/xv - 1)
    p = xm * tmp
    q = (1 - xm) * tmp

    return (p, q, loc, scale)  #check return type and should fixed be returned ? 

def _fitstart_poisson(self, x, fixed=None):
    '''maximum likelihood estimator as starting values for Poisson distribution

    Parameters
    ----------
    x : array
        data for which the parameters are estimated
    fixed : None or array_like
        sequence of numbers and np.nan to indicate fixed parameters and parameters
        to estimate

    Returns
    -------
    est : tuple
        preliminary estimates used as starting value for fitting, not
        necessarily a consistent estimator

    Notes
    -----
    This needs to be written and attached to each individual distribution

    References
    ----------
    MLE :
    http://en.wikipedia.org/wiki/Poisson_distribution#Maximum_likelihood

    '''
    #todo: separate out this part to be used for other compact support distributions
    #      e.g. rdist, vonmises, and truncnorm
    #      but this might not work because it might still be distribution specific
    a = x.min()
    eps = 0 # is this robust ?
    if fixed is None:
        #this part not checked with books
        loc = a - eps
    else:
        if np.isnan(fixed[-1]):
            #estimate loc
            loc = a - eps
        else:
            loc = fixed[-1]

    #MLE for standard (unshifted, if loc=0) Poisson distribution

    xtrans = (x - loc)
    lambd = xtrans.mean()
    #second derivative d loglike/ dlambd Not used
    #dlldlambd = 1/lambd # check

    return (lambd, loc)  #check return type and should fixed be returned ? 

def _fitstart_beta(self, x, fixed=None):
    '''method of moment estimator as starting values for beta distribution

    Parameters
    ----------
    x : array
        data for which the parameters are estimated
    fixed : None or array_like
        sequence of numbers and np.nan to indicate fixed parameters and parameters
        to estimate

    Returns
    -------
    est : tuple
        preliminary estimates used as starting value for fitting, not
        necessarily a consistent estimator

    Notes
    -----
    This needs to be written and attached to each individual distribution

    References
    ----------
    for method of moment estimator for known loc and scale
    http://en.wikipedia.org/wiki/Beta_distribution#Parameter_estimation
    http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm
    NIST reference also includes reference to MLE in
    Johnson, Kotz, and Balakrishan, Volume II, pages 221-235

    '''
    #todo: separate out this part to be used for other compact support distributions
    #      e.g. rdist, vonmises, and truncnorm
    #      but this might not work because it might still be distribution specific
    a, b = x.min(), x.max()
    eps = (a-b)*0.01
    if fixed is None:
        #this part not checked with books
        loc = a - eps
        scale = (a - b) * (1 + 2*eps)
    else:
        if np.isnan(fixed[-2]):
            #estimate loc
            loc = a - eps
        else:
            loc = fixed[-2]
        if np.isnan(fixed[-1]):
            #estimate scale
            scale = (b + eps) - loc
        else:
            scale = fixed[-1]

    #method of moment for known loc scale:
    scale = float(scale)
    xtrans = (x - loc)/scale
    xm = xtrans.mean()
    xv = xtrans.var()
    tmp = (xm*(1-xm)/xv - 1)
    p = xm * tmp
    q = (1 - xm) * tmp

    return (p, q, loc, scale)  #check return type and should fixed be returned ? 

def _fitstart_poisson(self, x, fixed=None):
    '''maximum likelihood estimator as starting values for Poisson distribution

    Parameters
    ----------
    x : array
        data for which the parameters are estimated
    fixed : None or array_like
        sequence of numbers and np.nan to indicate fixed parameters and parameters
        to estimate

    Returns
    -------
    est : tuple
        preliminary estimates used as starting value for fitting, not
        necessarily a consistent estimator

    Notes
    -----
    This needs to be written and attached to each individual distribution

    References
    ----------
    MLE :
    http://en.wikipedia.org/wiki/Poisson_distribution#Maximum_likelihood

    '''
    #todo: separate out this part to be used for other compact support distributions
    #      e.g. rdist, vonmises, and truncnorm
    #      but this might not work because it might still be distribution specific
    a = x.min()
    eps = 0 # is this robust ?
    if fixed is None:
        #this part not checked with books
        loc = a - eps
    else:
        if np.isnan(fixed[-1]):
            #estimate loc
            loc = a - eps
        else:
            loc = fixed[-1]

    #MLE for standard (unshifted, if loc=0) Poisson distribution

    xtrans = (x - loc)
    lambd = xtrans.mean()
    #second derivative d loglike/ dlambd Not used
    #dlldlambd = 1/lambd # check

    return (lambd, loc)  #check return type and should fixed be returned ?