In [2]:

    

%matplotlib inline
import sys
import os.path
sys.path.append(os.path.join(os.pardir,os.pardir))
import disaggregator.GreenButtonDatasetAdapter as gbda
import matplotlib.pyplot as plt
import pymc


    

In [12]:

    

xml_path = '../../../xml_data.xml'
with open(xml_path,'r') as f:
    xml_string = f.read()

trace = gbda.get_trace_from_xml(xml_string)


    

In [13]:

    

plt.plot(trace.series)
plt.show()
hist = plt.hist(trace.series.astype(float),bins=150)


    

In [14]:

    

a,a,a = plt.hist(trace.series.astype(float),bins=150)
plt.yscale('log')


    

In [15]:

    

a,a,a = plt.hist(trace.series.diff().abs().fillna(0).astype(float),bins=150)
plt.yscale('log')


    

In [28]:

    

data = trace.series.astype(float).values[:100]/1000
print data


    

    




[ 1.714  0.362  0.382  0.17   0.114  0.169  0.132  0.156  0.12   0.144
  0.117  0.148  0.121  0.161  0.113  0.155  0.099  0.315  0.233  0.53
  1.841  3.232  1.934  0.958  1.002  1.107  1.28   1.38   0.121  0.483
  0.79   1.     0.929  0.911  1.113  1.125  1.032  0.857  0.795  0.851
  0.791  0.492  0.669  0.375  0.382  0.317  0.132  2.35   0.129  0.081
  0.135  0.079  0.14   0.092  0.137  0.092  0.174  0.123  0.154  0.113
  0.113  0.113  0.103  0.116  0.096  0.119  0.091  0.127  0.093  0.417
  0.387  0.773  0.695  0.829  0.96   1.166  1.311  1.263  1.209  1.41
  1.356  1.27   1.511  1.173  1.205  1.281  1.073  0.926  0.812  0.651
  0.777  0.428  0.63   0.34   0.402  0.095  0.162  0.357  0.136  0.322]

In [41]:

    

n_states = 3
mean_lower_bounds = [0,0,0]
mean_upper_bounds = [1,1,1]
mean_bounds = zip(mean_lower_bounds,mean_upper_bounds)
beta_lower_bounds = [0,0,0]
beta_upper_bounds = [1,1,1]
beta_bounds = zip(beta_lower_bounds,beta_upper_bounds)
prob_priors = [.9,.5,.1]



alpha_params = [pymc.Uniform('mean_param_{}'.format(i),lower=lower, upper=upper) for i,(lower,upper) in enumerate(mean_bounds)]
beta_params = [pymc.Uniform('std_param_{}'.format(i),lower=lower, upper=upper) for i,(lower,upper) in enumerate(beta_bounds)]
appliance_ab = zip(alpha_params,beta_params)

appliances = [pymc.Gamma('appliance_params_{}'.format(i),alpha=alpha,beta=beta) for i,(alpha,beta) in enumerate(appliance_ab)]

prob_params = [pymc.Bernoulli('prob_param_{}'.format(i),p=p) for i,p in enumerate(prob_priors)]

@pymc.stochastic
def actual_output(appliances=appliances,probs=prob_params,value=data,observed=True):
    out = 0
    for appliance,prob in zip(appliances,probs):
        if prob:
            out = out + appliance
    print out
    return out

@pymc.stochastic
def pred_output(appliances=appliances,probs=prob_params,value=.2):
    out = 0
    for appliance,prob in zip(appliances,probs):
        if prob:
            out = out + appliance
    return out

M = pymc.MCMC(pymc.Model([pred_output,actual_output]))


    

    




0.547678669603

In [1]:

    

M.sample(1000,500,10)
print


    

In [40]:

    

pymc.Matplot.plot(M)


    

    




Plotting appliance_params_2
Plotting prob_param_0
Plotting appliance_params_0
Plotting appliance_params_1
Plotting prob_param_1
Plotting pred_output
Plotting prob_param_2






    













    













    













    













    













    













    

In [ ]: