In [22]:
#you have to make sure that you have all of these installed
import cProfile
import re
import math
import numpy as np
import scipy as sp
from scipy import stats
from scipy import optimize as opt
import pandas as pd
import random as rnd
from matplotlib import pyplot as plt
import time
import numpy.random
import warnings
warnings.filterwarnings('ignore')
import multiprocessing as mp
In [2]:
import chen_utils as ch
import user_simulation_utils as sim
In [3]:
rate_input = 0.3
share_observed = 0.7
observations_total = 7
print(rate_input, share_observed, observations_total)
In [4]:
type(sim.simulate_data()[0])
Out[4]:
In [5]:
x = sim.simulate_data()
p = 0.5
beta = 0.5
#runtime estimate
x
Out[5]:
In [6]:
#check function call
ch._gamma_calc(x[1],2 * (np.array(range(101))+1),0.5).shape
Out[6]:
In [9]:
np.random.seed(1234)
true_param = [0.5,0.5]
x = sim.simulate_data(true_param[0], true_param[1], 100)
ch.maximum_likelihood_estimate(x)
Out[9]:
In [10]:
ch.total_pipeline(x, true_param)
In [11]:
p = 0.6
beta = 0.4
x = sim.simulate_data(beta,p, 60)
traj = ch.metropolis(x, starting_point = [0.5,0.5],chain_length = 5000, burn_in = 0.2)
print(traj)
In [10]:
x
Out[10]:
In [11]:
print(np.mean(traj[:,0]), np.mean(traj[:,1]))
In [12]:
traj_new = ch.metropolis_new(x, starting_point = [0.5,0.5],chain_length = 5000, burn_in = 0.2)
In [13]:
np.mean(traj_new, axis = 0)
Out[13]:
In [14]:
x
Out[14]:
In [29]:
cProfile.run('ch.metropolis_new(x, starting_point = [0.5,0.5],chain_length = 5000, burn_in = 0.1)')
In [30]:
%timeit ch.metropolis_new(x, starting_point = [0.5,0.5],chain_length = 5000, burn_in = 0.1)
In [15]:
np.max(traj_new[:,1])
Out[15]:
In [16]:
ch.plot_trajectory(traj)
In [17]:
ch.plot_trajectory(traj_new)
In [18]:
ch.plot_region(traj, x, 80.0)
In [19]:
p = 0.6
beta = 0.3
x = sim.simulate_data(beta,p, 60)
time0 = time.time()
traj = ch.metropolis(x, starting_point = [0.5,0.5],chain_length = 100000, burn_in = 0.2)
print(time.time()-time0)
print(traj)
In [20]:
ch.plot_trajectory(traj)
ch.plot_region(traj, x, 90.0)
In [56]:
acceptedTraj = np.column_stack((sim.logit(traj[:,0]), np.log(traj[:,1])))
# Display the sampled points
"""Consider square shaped plot"""
# par( pty="s" ) # makes plots in square axes.
XY = np.transpose(acceptedTraj)
plt.plot(XY[0], XY[1], 'b')
plt.xlim([-3.0, 3.0])
plt.ylim([-3.0, 3.0])
plt.xlabel('P')
plt.ylabel(r'$\beta$')
# Display means and rejected/accepted ratio in plot.
if (meanTraj[0] > .5):
xpos, xadj = 0.01, 0.01
else:
xpos, xadj = 0.95, 0.95
if (meanTraj[1] > .5):
ypos, yadj = 0.01, 0.01
else:
ypos, yadj = 0.95, 0.95
plt.show()