Skip to content
Snippets Groups Projects
Commit e6e68a8e authored by Jörg Martin's avatar Jörg Martin
Browse files

Expanded uncertainties added

parent 71bd85b2
No related branches found
No related tags found
No related merge requests found
...@@ -156,6 +156,19 @@ class poi_ssd_framework(generic_ssd_framework): ...@@ -156,6 +156,19 @@ class poi_ssd_framework(generic_ssd_framework):
post_std = np.sqrt(self.alpha + n * mean)/(self.beta + n) post_std = np.sqrt(self.alpha + n * mean)/(self.beta + n)
return post_mean, post_std return post_mean, post_std
def get_post_quantiles(self, mean, n):
"""
Returns the median, lower_expanded_uncertainty, and upper_expanded_uncertainty given the sufficient statistic `mean` and the sample size `n`.
"""
post_alpha = self.alpha + mean
post_beta = self.beta + n
ppf = lambda q: scipy.stats.gamma.ppf(q, a=post_alpha)/post_beta
median = ppf(0.5)
lower_expanded_uncertainty = median - ppf(0.025)
upper_expanded_uncertainty = ppf(0.975) - median
return median, lower_expanded_uncertainty, upper_expanded_uncertainty
class normal_ssd_framework(generic_ssd_framework): class normal_ssd_framework(generic_ssd_framework):
""" """
...@@ -236,3 +249,23 @@ class normal_ssd_framework(generic_ssd_framework): ...@@ -236,3 +249,23 @@ class normal_ssd_framework(generic_ssd_framework):
+ n/(lamb*n_lamb**2*n_alpha) * (mean - mu_0)**2 + n/(lamb*n_lamb**2*n_alpha) * (mean - mu_0)**2
) )
return post_mean, post_std return post_mean, post_std
def get_post_quantiles(self, mean, std, n):
"""
Returns the median, lower expanded uncertainty, and upper expanded uncertainty given the sufficient statistics `mean`, `std` and the sample size `n`.
**Note**: For compatability both, lower and upper expanded_uncertainty are specified although they are identical for the Student t marginal.
"""
# use post mean and std
post_mean, post_std = self.make_inference(mean, std, n)
# compute median
median = post_mean
# compute expanded uncertainty
post_alpha = self.alpha + n/2
df = 2 * post_alpha
std_t = np.sqrt(df/(df-2))
expansion_factor = scipy.stats.t.ppf(0.975, df=df)/std_t
expanded_uncertainty = post_std * expansion_factor
lower_expanded_uncertainty = expanded_uncertainty
upper_expanded_uncertainty = expanded_uncertainty
return median, lower_expanded_uncertainty, upper_expanded_uncertainty
...@@ -8,4 +8,3 @@ from vpvc_interface.gui import run_gui ...@@ -8,4 +8,3 @@ from vpvc_interface.gui import run_gui
run_gui(fontsize=11) run_gui(fontsize=11)
...@@ -153,28 +153,45 @@ class InferenceMenu(): ...@@ -153,28 +153,45 @@ class InferenceMenu():
self.stat_entries = (self.mean_entry, ) self.stat_entries = (self.mean_entry, )
self.computation_frame = tk.Frame(self.master) self.computation_frame = tk.Frame(self.master)
self.compute_button = tk.Button(self.computation_frame, self.compute_button = tk.Button(self.computation_frame,
text=button_text, command=lambda: self.compute_inference(ssd)) text=button_text, command=lambda: self.compute_inference(ssd, distribution_name))
self.compute_button.grid(row=0, column=0, sticky=tk.W) self.compute_button.grid(row=0, column=0, sticky=tk.W)
self.estimate_text = tk.Text(self.computation_frame, height=1, width=30, bd=0, font=self.font, state='disabled') self.estimate_text = tk.Text(self.computation_frame, height=1, width=50, bd=0, font=self.font, state='disabled')
self.uncertainty_text = tk.Text(self.computation_frame, height=1, bd=0, width=30, font=self.font, state='disabled') self.uncertainty_text = tk.Text(self.computation_frame, height=1, bd=0, width=50, font=self.font, state='disabled')
self.estimate_text.grid(row=1, column=0, sticky=tk.E) self.estimate_text.grid(row=1, column=0, sticky=tk.E)
self.uncertainty_text.grid(row=2, column=0, sticky=tk.E) self.uncertainty_text.grid(row=2, column=0, sticky=tk.E)
# quantiles of the posterior
if distribution_name == 'Normal':
self.expanded_uncertainty_text = tk.Text(self.computation_frame, height=1, bd=0, width=50, font=self.font, state='disabled')
self.expanded_uncertainty_text.grid(row=3, column=0, sticky=tk.E)
if distribution_name == 'Poisson':
self.median_text = tk.Text(self.computation_frame, height=1, width=50, bd=0, font=self.font, state='disabled')
self.upper_expanded_uncertainty_text = tk.Text(self.computation_frame, height=1, bd=0, width=50, font=self.font, state='disabled')
self.lower_expanded_uncertainty_text = tk.Text(self.computation_frame, height=1, bd=0, width=50, font=self.font, state='disabled')
self.median_text.grid(row=3, column=0, sticky=tk.E)
self.lower_expanded_uncertainty_text.grid(row=4, column=0, sticky=tk.E)
self.upper_expanded_uncertainty_text.grid(row=5, column=0, sticky=tk.E)
self.entry_frame.grid(row=0,column=0, sticky=tk.W) self.entry_frame.grid(row=0,column=0, sticky=tk.W)
self.computation_frame.grid(row=1, column=0, sticky=tk.W) self.computation_frame.grid(row=1, column=0, sticky=tk.W)
def compute_inference(self, ssd): def compute_inference(self, ssd, distribution_name):
# preliminiaries for pretty printing
## return the magnitude of a number
magnitude = lambda x: np.floor(np.log10(np.abs(x)))
## format string into pretty printing
def pretty_string(prefix, x):
mag_x = magnitude(x)
if mag_x > 4 or mag_x < -4:
return prefix + '%e' % (x,)
else:
return prefix + '%f' % (x,)
# make inference
try: try:
args = [float(e.get()) for e in self.stat_entries] args = [float(e.get()) for e in self.stat_entries]
n = int(self.n_entry.get()) n = int(self.n_entry.get())
## posterior mean and std
estimate, uncertainty = ssd.make_inference(*args, n) estimate, uncertainty = ssd.make_inference(*args, n)
if estimate != 0: ## median and expanded uncertainties
estimate_magnitude = np.floor(np.log10(np.abs(estimate))) median, lower_expanded_uncertainty, upper_expanded_uncertainty = ssd.get_post_quantiles(*args, n)
else:
estimate_magnitude = 0
if uncertainty != 0:
uncertainty_magnitude = np.floor(np.log10(uncertainty))
else:
uncertainty_magnitude = 0
except (ValueError, AssertionError): except (ValueError, AssertionError):
self.status_variable.set('Invalid values encountered') self.status_variable.set('Invalid values encountered')
raise ValueError('Invalid values encountered') raise ValueError('Invalid values encountered')
...@@ -183,16 +200,30 @@ class InferenceMenu(): ...@@ -183,16 +200,30 @@ class InferenceMenu():
self.status_variable.set('Bayesian inference') self.status_variable.set('Bayesian inference')
self.uncertainty_text.delete("1.0",tk.END) self.uncertainty_text.delete("1.0",tk.END)
self.estimate_text.delete("1.0",tk.END) self.estimate_text.delete("1.0",tk.END)
if estimate_magnitude > 4 or estimate_magnitude < -4: self.estimate_text.insert(tk.END, pretty_string('Estimate: ', estimate))
self.estimate_text.insert(tk.END, 'Estimate: %e' % (estimate, )) self.uncertainty_text.insert(tk.END, pretty_string('Uncertainty: ', uncertainty))
else:
self.estimate_text.insert(tk.END, 'Estimate: %f' % (estimate, ))
if uncertainty_magnitude > 4 or uncertainty_magnitude < -4:
self.uncertainty_text.insert(tk.END, 'Uncertainty: %e' % (uncertainty, ))
else:
self.uncertainty_text.insert(tk.END, 'Uncertainty: %f' % (uncertainty, ))
self.estimate_text.configure(state='disabled') self.estimate_text.configure(state='disabled')
self.uncertainty_text.configure(state='disabled') self.uncertainty_text.configure(state='disabled')
if distribution_name == 'Normal':
self.expanded_uncertainty_text.configure(state='normal')
self.expanded_uncertainty_text.delete("1.0",tk.END)
self.expanded_uncertainty_text.insert(tk.END, pretty_string('Expanded uncertainty: ', upper_expanded_uncertainty))
self.expanded_uncertainty_text.configure(state='disabled')
if distribution_name == 'Poisson':
self.median_text.configure(state='normal')
self.lower_expanded_uncertainty_text.configure(state='normal')
self.upper_expanded_uncertainty_text.configure(state='normal')
self.median_text.delete("1.0",tk.END)
self.lower_expanded_uncertainty_text.delete("1.0",tk.END)
self.upper_expanded_uncertainty_text.delete("1.0",tk.END)
self.median_text.insert(tk.END, pretty_string('Median: ', median))
self.lower_expanded_uncertainty_text.insert(tk.END, pretty_string('L. exp. uncertainty: ', lower_expanded_uncertainty))
self.upper_expanded_uncertainty_text.insert(tk.END, pretty_string('U. exp. uncertainty: ', upper_expanded_uncertainty))
self.median_text.configure(state='disabled')
self.lower_expanded_uncertainty_text.configure(state='disabled')
self.upper_expanded_uncertainty_text.configure(state='disabled')
class ResultMenu(): class ResultMenu():
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment