1# -*- coding: utf-8 -*- 2""" 3 4Created on Sun Apr 21 07:59:26 2013 5 6Author: Josef Perktold 7""" 8 9from statsmodels.compat.python import lmap 10import numpy as np 11import matplotlib.pyplot as plt 12 13import statsmodels.stats.proportion as sms 14import statsmodels.stats.weightstats as smw 15 16from numpy.testing import assert_almost_equal 17 18 19# Region, Eyes, Hair, Count 20ss = '''\ 211 blue fair 23 1 blue red 7 1 blue medium 24 221 blue dark 11 1 green fair 19 1 green red 7 231 green medium 18 1 green dark 14 1 brown fair 34 241 brown red 5 1 brown medium 41 1 brown dark 40 251 brown black 3 2 blue fair 46 2 blue red 21 262 blue medium 44 2 blue dark 40 2 blue black 6 272 green fair 50 2 green red 31 2 green medium 37 282 green dark 23 2 brown fair 56 2 brown red 42 292 brown medium 53 2 brown dark 54 2 brown black 13''' 30 31dta0 = np.array(ss.split()).reshape(-1,4) 32dta = np.array(lmap(tuple, dta0.tolist()), dtype=[('Region', int), ('Eyes', 'S6'), ('Hair', 'S6'), ('Count', int)]) 33 34xfair = np.repeat([1,0], [228, 762-228]) 35 36# comparing to SAS last output at 37# http://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/viewer.htm#procstat_freq_sect028.htm 38# confidence interval for tost 39ci01 = smw.confint_ztest(xfair, alpha=0.1) 40assert_almost_equal(ci01, [0.2719, 0.3265], 4) 41res = smw.ztost(xfair, 0.18, 0.38) 42 43assert_almost_equal(res[1][0], 7.1865, 4) 44assert_almost_equal(res[2][0], -4.8701, 4) 45 46nn = np.arange(200, 351) 47pow_z = sms.power_ztost_prop(0.5, 0.72, nn, 0.6, alpha=0.05) 48pow_bin = sms.power_ztost_prop(0.5, 0.72, nn, 0.6, alpha=0.05, dist='binom') 49plt.plot(nn, pow_z[0], label='normal') 50plt.plot(nn, pow_bin[0], label='binomial') 51plt.legend(loc='lower right') 52plt.title('Proportion Equivalence Test: Power as function of sample size') 53plt.xlabel('Number of Observations') 54plt.ylabel('Power') 55 56plt.show() 57