#!/usr/bin/env python2 # -*- coding: utf-8 -*- """ Created on Wed Oct 10 14:06:21 2018 @author: bob Two way BANOVA """ from __future__ import division import numpy as np import pymc3 as pm import pandas as pd import matplotlib.pyplot as plt plt.style.use('seaborn-darkgrid') from scipy.stats import norm from sklearn.preprocessing import LabelEncoder from theano import tensor as tt from sys import exit # =============== load data ============================ # single measurements dataF = pd.read_csv("../Data/FLT1-7_50x.csv") dataQ = pd.read_csv("../Data/QTFU2-10_50x.csv") dataR = pd.read_csv("../Data/RS0.4u_50x_4x1tiles.csv") # unite them data = pd.concat([dataF,dataQ], ignore_index=True,sort=True) # measures measuresAll = ["Sq","Ssk","Sku","Sp","Sv","Sz","Sa","Smr","Smc","Sxp","Sal","Str","Std","Sdq","Sdr","Vm","Vv","Vmp","Vmc","Vvc","Vvv","Sk","Spk","Svk","Smr1","Smr2","Svq","Smq","Spq"] measuresSuggested = ["Sa", "Sdr", "Str", "Std", "Smr", "Vmc"] measuresSingle = ['Sa'] measuresBig = ['Sdr','Smr'] # get only the needed rows and columns for idMeasure in measuresAll: print idMeasure filtered = data[['Objective','Location','Sample',idMeasure]] # transform them to standard form melt = pd.DataFrame() melt["Factor1"] = filtered['Objective'] melt["Factor2"] = filtered["Location"] + filtered["Sample"] melt['Value'] = filtered[idMeasure] print melt #estimate magnitude of variance due to change of objective and change of other factor estimatorDF_1 = melt.groupby('Factor2').std() print estimatorDF_1 stdGuessFactor1 = estimatorDF_1['Value'].max() print "Std guess b1 hyper:",stdGuessFactor1 estimatorDF_2 = melt.groupby('Factor1').std() print estimatorDF_2 stdGuessFactor2 = estimatorDF_2['Value'].max() print "Std guess b2 hyper:",stdGuessFactor2 combinedStd = np.sqrt(np.power(stdGuessFactor1,2)+np.power(stdGuessFactor2,2)) # extract from data frame y = melt['Value'].values x1 = pd.Categorical(melt['Factor1']).codes x1names = melt['Factor1'].unique() x2 = pd.Categorical(melt['Factor2']).codes x2names = melt['Factor2'].unique() # ============ prepare for contrasts Ntotal = len(y) Nx1Lvl = len(set(x1)) Nx2Lvl = len(set(x2)) x1contrast_dict = {'Objective_0_vs_1': [1, -1]} # # ============ compute mean and standard devs meanData = np.mean(y) stdData = np.std(y) # ================== specify model y ~ N(mu,std) with mu = b0 + b1*x1+ b2*x2 + x1*M*x2 ================================ with pm.Model() as model: # define the priors b0 = pm.Deterministic('b0', meanData +stdData * pm.Normal('b0-Tilde',mu=0,sd=1)) b1 = pm.Normal('b1', mu=0, sd=stdGuessFactor1, shape=Nx1Lvl) b2 = pm.Normal('b2', mu=0, sd=stdGuessFactor2, shape=Nx2Lvl) M = pm.Normal('M' , mu=0, sd=0.05*combinedStd ,shape=[Nx1Lvl, Nx2Lvl]) mu = pm.Deterministic('mu',b0 + b1[x1]+ b2[x2] + M[x1,x2]) sigma = pm.Uniform('Like_std', upper=0.2*np.min([stdGuessFactor1,stdGuessFactor2])) # define likelihood like = pm.Normal('like', mu=mu, sd=sigma, observed=y) # Generate a MCMC chain trace = pm.sample(6000,cores=6,tune=4000,target_accept=0.99, init='auto') # =================== results ================================= interestingVars = ['b0','b1','b2','M'] # Print summary for each trace print pm.summary(trace).round(2) # Check for mixing and autocorrelation pm.energyplot(trace) # vars throws exception plt.savefig('./MeasuredQuantities/'+idMeasure + '_Auto.pdf') ## Plot KDE and sampled values for each parameter. #pm.traceplot(trace,varnames=interestingVars) pm.traceplot(trace) plt.savefig('./MeasuredQuantities/'+idMeasure + '_Trace.pdf') if False: pm.pairplot(trace, sub_varnames=['b0','Hyper_b1_std'], divergences=True, color='C3', figsize=(10, 5), kwargs_divergence={'color':'C2'}) plt.title('scatter plot between log(tau) and theta[0]') plt.savefig('./MeasuredQuantities/'+idMeasure + '_Diagnostic.pdf') # Extract values of 'a' b0_sample = trace['b0'] b1_sample = trace['b1'] b2_sample = trace['b2'] b1b2_sample = trace['M'] # =========== plotting ======================= # compute figure size figBaseLength = 3 figSizeX = len(x1names)+1 figSizeY = len(x2names)+1 # create figure plt.figure(figsize=(figSizeX*figBaseLength,figSizeY*figBaseLength)) # (width, height) for i in range(figSizeX*figSizeY): ax = plt.subplot(figSizeY, figSizeX, i+1) ax.set_title(r'${}$'.format(i)) # plot baseline ax = plt.subplot(figSizeY,figSizeX,1) # nrows, ncolumns,index pm.plot_posterior(b0_sample, bins=50, ax=ax) ax.set_title(r'$\beta0$') # plot first predictor coefficients for j in range(len(x1names)): ax = plt.subplot(figSizeY, figSizeX, j+2) pm.plot_posterior(b1_sample[:,j], ax=ax) ax.set_title(r'$\beta1_{}$'.format(j)) # plot second predictor coefficients for i in range(len(x2names)): ax = plt.subplot(figSizeY, figSizeX, (i+1)*figSizeX+1) pm.plot_posterior(b2_sample[:,i], bins=50, ax=ax) ax.set_title(r'$\beta2_{}$'.format(i)) # plot combinated predictor coefficients for j in range(len(x1names)): for i in range(len(x2names)): ax = plt.subplot(figSizeY, figSizeX, (i+1)*figSizeX+2+j) pm.plot_posterior(b1b2_sample[:,j,i], bins=50, ax=ax) ax.set_title(r'$M$ at $i={} j={}$'.format(i, j)) plt.tight_layout() plt.savefig('./MeasuredQuantities/'+idMeasure + '_HistoMatrix.pdf') ## Display contrast analyses figBaseLengthContrasts = 5 numberContrasts = len(x1contrast_dict.items()) plt.figure(figsize=(numberContrasts*figBaseLengthContrasts,figBaseLengthContrasts)) count = 1 for key, value in x1contrast_dict.items(): contrast = np.dot(b1_sample, value) ax = plt.subplot(1,numberContrasts, count) pm.plot_posterior(contrast, ref_val=0.0, bins=50, ax=ax) ax.set_title('Contrast {}'.format(key)) count += 1 plt.tight_layout() plt.savefig('./MeasuredQuantities/'+idMeasure + '_Contrasts.pdf')