simd/statistics.py at master · fomy/simd · GitHub

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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#
# A Class that couples samples from a simulation
# with some statistics functions
#
# Kevin Greenan (kmgreen@cs.ucsc.edu)
#
# Improved by Min Fu (fumin@hust.edu.cn)
import random
import math

#
# A Class that incapsulates a set of samples with
# operations over those samples (i.e. statistics)
#
class Samples:

    #
    # Construct new instance with a list of samples
    #
    # @param samples: a set of samples, most observed in simulation
    #
    def __init__(self):

        self.value_sum = 0L
        self.value2_sum = 0L
        self.prob_sum = 0L

        self.num_samples = 0

        #
        # A static table used to estimate the confidence
        # interval around a sample mean
        #
        self.conf_lvl_lku = {}
        self.conf_lvl_lku["0.80"] = 1.281
        self.conf_lvl_lku["0.85"] = 1.440
        self.conf_lvl_lku["0.90"] = 1.645
        self.conf_lvl_lku["0.95"] = 1.960
        self.conf_lvl_lku["0.995"] = 2.801

        self.value_mean = None
        self.value_mean2 = None
        self.value_dev = None
        self.value_ci = None
        self.value_re = None

        # used to calculate the probability of data loss
        self.prob_mean = None
        self.prob_mean2 = None
        self.prob_dev = None
        self.prob_ci = None
        self.prob_re = None

    # only non-zeros in samples, num shows the actual number of samples */
    def addSamples(self, samples, num):
        for sample in samples:
            self.value_sum += sample
            self.value2_sum += pow(sample, 2)
            self.prob_sum += 1

        self.num_samples += num

    def addSample(self, sample):
        if sample > 0:
            self.value_sum += sample
            self.value2_sum += pow(sample, 2)
            self.prob_sum += 1

        self.num_samples += 1L

    def addZeros(self, num):
        self.num_samples += long(num)

    #
    # Calculate the sample mean based on the samples for this instance
    #
    def calcMean(self):
        self.value_mean = self.value_sum / float(self.num_samples)
        self.value2_mean = self.value2_sum / float(self.num_samples)

        self.prob_mean = self.prob_sum / float(self.num_samples)

    #
    # Calculate the standard deviation based on the samples for this instance
    # dev = E(X-EX)^2 = EX^2 - (EX)^2
    #
    def calcStdDev(self ):
        self.calcMean()
        self.value_dev = math.sqrt(self.value2_mean - pow(self.value_mean, 2))
        self.prob_dev = math.sqrt(self.prob_mean - pow(self.prob_mean, 2))

    #
    # Calculate the confidence interval around the sample mean
    #
    # @param conf_level: the probability that the mean falls within the interval
    #
    def calcConfInterval(self, conf_level):

        if conf_level not in self.conf_lvl_lku.keys():
            print "%s not a valid confidence level!" % conf_level
            return None

        self.calcStdDev()

        self.value_ci = abs(self.conf_lvl_lku[conf_level] * (self.value_dev / math.sqrt(self.num_samples)))
        self.prob_ci = abs(self.conf_lvl_lku[conf_level] * (self.prob_dev / math.sqrt(self.num_samples)))

    #
    # Calculate the relative error
    #
    # self.conf_lvl_lku[conf_level] * sqrt(Var)/sqrt(num_samples) / mean
    #
    def calcRE(self, conf_level):

        self.calcConfInterval(conf_level)

        if self.value_mean == 0:
            self.value_re = 0
            self.prob_re = 0
        else:
            self.value_re = self.value_ci / self.value_mean
            self.prob_re = self.prob_ci / self.prob_mean

    # zeros have been eliminated
    def calcResults(self, conf_level):

        self.calcRE(conf_level)

#
# Generate samples from a known distribution and verify the statistics
#
def test():

    num_samples = 1000
    samples = []
    mean = 0.5
    std_dev = 0.001
    for i in range(num_samples):
        samples.append(random.gauss(mean, std_dev))

    s = Samples(1000)
    s.calcResults("0.9", samples)

    print "Mean: %s (%s): " % (s.calcMean(), mean)
    print "Std Dev: %s (%s): " % (s.calcStdDev(), std_dev)
    print "Conf. Interval: (%s, %s)" % s.calcConfInterval("0.995")

    (a,b,c,d) = s.getResults()
    print "Mean: %s (%s): " % (a, mean)
    print "Conf. Interval: (%s, %s)" % (b,c)
    print "Relative Error: (%s)" % d

if __name__ == "__main__":
    test()