jcbozonier

tiny_numbers = np.array([-19285006, -19275006, -19275003, -19275002, -19275001])
scaled = tiny_numbers - np.max(tiny_numbers)
print(scaled)

jcbozonier

import numpy as np

aggregate_product = np.log(1)
iterations = 0

while aggregate_product != float('-inf'):
    aggregate_product += np.log(0.5)
    iterations += 1
    
# break program execution and then print out

jcbozonier

aggregate_product = 1
iterations = 0

while aggregate_product != 0:
    aggregate_product *= 0.51
    iterations += 1
   
print(iterations)

jcbozonier

mle = {}

for var_name in ['a','h','k','sigma']:
    var_mass, var_bins = np.histogram(trace[var_name], bins=np.arange(-50000,50000, 0.1), density=True)
    var_averaged_bins = 0.5*(var_bins[1:] + var_bins[:-1])
    bin_index = np.argmax(var_mass)
    mle[var_name] = var_averaged_bins[bin_index]

print(mle);

jcbozonier

%matplotlib inline

import pymc3 as pm

with pm.Model() as model:
    _a = pm.Uniform('a', lower=-10, upper=0)
    _h = pm.Uniform('h', lower=4, upper=12)
    _k = pm.Uniform('k', lower=18000, upper=30000)
    _sigma = pm.Uniform('sigma', lower=5, upper=25)

jcbozonier

max(evaluated_hypotheses, key=lambda x: x[-1])

jcbozonier

def normalize(log_p):
  shifted_p = np.exp(log_p - np.max(log_p))
  normalized_log_p = shifted_p/shifted_p.sum()
  return normalized_log_p

jcbozonier

# We'll be working in log space as an optimization
# So that our numbers don't get so small that the computer 
# rounds them to zero.
# Start with an uninformative prior
# Every hypothesis gets an equal weighting.
# Not important that they sum to one, just that
# they're equally weighted.
hypothesis_likelihoods = np.log(np.array([1]))

for i, hypotheses_tuple in enumerate(hypotheses):

jcbozonier

# Come up with a hypothesis for each parameter
a_hypotheses = np.array([-4,-3,-2,-1])
h_hypotheses = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
k_hypotheses = np.array([18000, 20000, 23000, 25000, 26000, 28000, 30000])
sigma_hypotheses = np.array([5, 10, 15, 20, 25])

# Create list of every possible combination of them.
hypotheses = list(itertools.product(a_hypotheses, h_hypotheses, k_hypotheses, sigma_hypotheses))

jcbozonier

def reverse_aggregated_partition_list(partitions):
    # We're going to now merge partitions working from right
    # to left to try and simplify the partitioned data
    reversed_partition_list = list(reversed(partitions))
    aggregated_partition_list = []

    current_partition = reversed_partition_list[0]
    for i in range(len(reversed_partition_list)):
      if i+1 >= len(reversed_partition_list):
        # If we don't have a "next" partition to compare to