File: stats_isf.py

package info (click to toggle)
pyferret 7.0.0-2
  • links: PTS
  • area: main
  • in suites: stretch
  • size: 117,852 kB
  • ctags: 21,106
  • sloc: fortran: 214,303; ansic: 25,025; python: 23,717; java: 1,755; sh: 1,375; makefile: 978; pascal: 569; csh: 285; awk: 18
file content (88 lines) | stat: -rw-r--r-- 3,558 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
"""
Returns the array of inverse survival function values for
a probability distribution and set of quantile values.
"""
from __future__ import print_function
from past.builtins import xrange
import numpy
import scipy.stats
import pyferret
import pyferret.stats


def ferret_init(id):
    """
    Initialization for the stats_isf python-backed ferret external function
    """
    axes_values = [ pyferret.AXIS_IMPLIED_BY_ARGS ] * pyferret.MAX_FERRET_NDIM
    true_influences = [ True ] * pyferret.MAX_FERRET_NDIM
    false_influences = [ False ] * pyferret.MAX_FERRET_NDIM
    retdict = { "numargs": 3,
                "descript": "Returns inverse survival function values for a probability distribution",
                "axes": axes_values,
                "argnames": ("PROBS", "PDNAME", "PDPARAMS"),
                "argdescripts": ("Probabilities (0-1) at which to calculate the inverse survival function values",
                                 "Name of a probability distribution",
                                 "Parameters for this probability distribution"),
                "argtypes": (pyferret.FLOAT_ARRAY, pyferret.STRING_ONEVAL, pyferret.FLOAT_ARRAY),
                "influences": (true_influences, false_influences, false_influences),
              }
    return retdict


def ferret_compute(id, result, resbdf, inputs, inpbdfs):
    """
    Assigns result with the inverse survival function values for the probability
    distribution indicated by inputs[1] (a string) using the parameters given in
    inputs[2] at the quantile values given by inputs[0].  For undefined quantile
    values, the result value will be undefined.
    """
    distribname = inputs[1]
    distribparams = inputs[2].reshape(-1)
    distrib = pyferret.stats.getdistrib(distribname, distribparams)
    badmask = ( numpy.fabs(inputs[0] - inpbdfs[0]) < 1.0E-5 )
    badmask = numpy.logical_or(badmask, numpy.isnan(inputs[0]))
    goodmask = numpy.logical_not(badmask)
    result[badmask] = resbdf
    # array[goodmask] is a flattened array
    result[goodmask] = distrib.isf(inputs[0][goodmask])


#
# The rest of this is just for testing this module at the command line
#
if __name__ == "__main__":
    # make sure ferret_init does not have problems
    info = ferret_init(0)

    # Normal distribution along the Y axis
    dimen = 25
    mu = 5.0
    sigma = 2.0
    distf = scipy.stats.norm(mu, sigma)
    qvals = numpy.linspace(0.05, 0.95, dimen)
    isfvals = distf.isf(qvals)

    pfname = "norm"
    pfparams = numpy.array([mu, sigma], dtype=numpy.float64)
    inpbdfs = numpy.array([-1.0, 0.0, 0.0], dtype=numpy.float64)
    resbdf = numpy.array([-2.0], dtype=numpy.float64)
    quantile = numpy.empty((1, dimen, 1, 1, 1, 1), dtype=numpy.float64, order='F')
    expected = numpy.empty((1, dimen, 1, 1, 1, 1), dtype=numpy.float64, order='F')
    for j in xrange(dimen):
        if (j % 7) == 3:
            quantile[0, j, 0, 0, 0, 0] = inpbdfs[0]
            expected[0, j, 0, 0, 0, 0] = resbdf[0]
        else:
            quantile[0, j, 0, 0, 0, 0] = qvals[j]
            expected[0, j, 0, 0, 0, 0] = isfvals[j]
    result = -888.0 * numpy.ones((1, dimen, 1, 1, 1, 1), dtype=numpy.float64, order='F')
    ferret_compute(0, result, resbdf, (quantile, pfname, pfparams), inpbdfs)
    if not numpy.allclose(result, expected):
        print ("Expected (flattened) = %s" % str(expected.reshape(-1)))
        print ("Result (flattened) = %s" % str(result.reshape(-1)))
        raise ValueError("Unexpected result")

    # All successful
    print ("Success")