undid earlier changes in PDFStatistics

This commit is contained in:
sebastianp 2016-07-13 11:00:30 +02:00
parent 8cef8ff2db
commit 38d225fbcf
2 changed files with 30 additions and 13 deletions

View File

@ -339,6 +339,12 @@ class PDFDictionary(object):
class PDFstatistics(object): class PDFstatistics(object):
'''
To do:
plots for std, quantiles,
calc: quantiles
encountering error when trying to get quantile... float.. utcdatetime error
'''
def __init__(self, directory): def __init__(self, directory):
self.directory = directory self.directory = directory
self.stations = {} self.stations = {}
@ -347,19 +353,20 @@ class PDFstatistics(object):
self.makeDirlist() self.makeDirlist()
self.getData() self.getData()
self.getStatistics() self.getStatistics()
self.arraylen = 0
self.Theta015 = []
self.Theta25 = []
self.Theta485 = []
#self.showData() #self.showData()
def makeDirlist(self, fn_pattern='*'): def makeDirlist(self, fn_pattern='*'):
self.dirlist = glob.glob1(self.directory, fn_pattern) self.dirlist = glob.glob1(self.directory, fn_pattern)
#print self.dirlist
self.evtdict = {} self.evtdict = {}
for rd in self.dirlist: for rd in self.dirlist:
#print rd
self.evtdict[rd] = glob.glob1((os.path.join(self.directory, rd)), '*.xml') self.evtdict[rd] = glob.glob1((os.path.join(self.directory, rd)), '*.xml')
#print rd, self.evtdict[rd]
def getData(self): def getData(self):
arraylen = 0
for dir in self.dirlist: for dir in self.dirlist:
for evt in self.evtdict[dir]: for evt in self.evtdict[dir]:
self.stations[evt] = [] self.stations[evt] = []
@ -371,27 +378,30 @@ class PDFstatistics(object):
try: try:
p_std = pdfs['P'].standard_deviation() p_std = pdfs['P'].standard_deviation()
except KeyError: except KeyError:
p_std = 0 p_std = np.nan
try: try:
s_std = pdfs['S'].standard_deviation() s_std = pdfs['S'].standard_deviation()
except KeyError: except KeyError:
s_std = 0 s_std = np.nan
self.stations[evt].append(station) self.stations[evt].append(station)
self.p_std[evt].append(p_std) self.p_std[evt].append(p_std)
self.s_std[evt].append(s_std) self.s_std[evt].append(s_std)
arraylen += 1 self.arraylen += 1
self.makeArray(arraylen)
def makeArray(self, arraylen):
def makeArray(self):
index = 0 index = 0
self.p_stdarray = np.zeros(arraylen) self.p_stdarray = np.zeros(self.arraylen)
self.s_stdarray = np.zeros(arraylen) self.s_stdarray = np.zeros(self.arraylen)
for evt in self.p_std.keys(): for evt in self.p_std.keys():
for n in range(len(self.p_std[evt])): for n in range(len(self.p_std[evt])):
self.p_stdarray[index] = self.p_std[evt][n] self.p_stdarray[index] = self.p_std[evt][n]
self.s_stdarray[index] = self.s_std[evt][n] self.s_stdarray[index] = self.s_std[evt][n]
index += 1 index += 1
def showData(self): def showData(self):
#figure(p, s) = plt.pyplot.subplots(2) #figure(p, s) = plt.pyplot.subplots(2)
#p.hist(self.p_stdarray, Bins=100, range=[min(self.p_stdarray),max(self.p_stdarray)/256]) #p.hist(self.p_stdarray, Bins=100, range=[min(self.p_stdarray),max(self.p_stdarray)/256])
@ -405,6 +415,7 @@ class PDFstatistics(object):
def getPDFDict(self, month, evt): def getPDFDict(self, month, evt):
self.pdfdict = PDFDictionary.from_quakeml(os.path.join(self.directory,month,evt)) self.pdfdict = PDFDictionary.from_quakeml(os.path.join(self.directory,month,evt))
def getStatistics(self): def getStatistics(self):
self.p_mean = self.p_stdarray.mean() self.p_mean = self.p_stdarray.mean()
self.p_std_std = self.p_stdarray.std() self.p_std_std = self.p_stdarray.std()
@ -412,3 +423,8 @@ class PDFstatistics(object):
self.s_mean = self.s_stdarray.mean() self.s_mean = self.s_stdarray.mean()
self.s_std_std = self.s_stdarray.std() self.s_std_std = self.s_stdarray.std()
self.s_median = np.median(self.s_stdarray) self.s_median = np.median(self.s_stdarray)
if __name__ == "__main__":
rootdir = '/home/sebastianp/Data/Reassessment/Insheim/'
Insheim = PDFstatistics(rootdir)

View File

@ -105,6 +105,7 @@ def exp_branches(k, (mu, sig1, sig2, a)):
''' '''
def _func(k, mu, sig1, sig2, a): def _func(k, mu, sig1, sig2, a):
mu = float(mu)
if k < mu: if k < mu:
rval = a * np.exp(sig1 * (k - mu)) rval = a * np.exp(sig1 * (k - mu))
else: else:
@ -311,7 +312,7 @@ class ProbabilityDensityFunction(object):
mu = self.mu mu = self.mu
rval = 0 rval = 0
for x in self.axis: for x in self.axis:
rval += (x - mu) ** 2 * self.data(x) rval += (x - float(mu)) ** 2 * self.data(x)
return rval * self.incr return rval * self.incr
def prob_lt_val(self, value): def prob_lt_val(self, value):