1# -*- coding: utf-8 -*- 2""" 3 4Created on Wed Feb 19 12:39:49 2014 5 6Author: Josef Perktold 7""" 8 9import numpy as np 10from scipy import stats 11 12from statsmodels.sandbox.distributions.extras import (SkewNorm_gen, skewnorm, 13 ACSkewT_gen, 14 NormExpan_gen, pdf_moments, 15 ExpTransf_gen, LogTransf_gen) 16from statsmodels.stats.moment_helpers import mc2mvsk, mnc2mc, mvsk2mnc 17 18 19def example_n(): 20 21 print(skewnorm.pdf(1,0), stats.norm.pdf(1), skewnorm.pdf(1,0) - stats.norm.pdf(1)) 22 print(skewnorm.pdf(1,1000), stats.chi.pdf(1,1), skewnorm.pdf(1,1000) - stats.chi.pdf(1,1)) 23 print(skewnorm.pdf(-1,-1000), stats.chi.pdf(1,1), skewnorm.pdf(-1,-1000) - stats.chi.pdf(1,1)) 24 rvs = skewnorm.rvs(0,size=500) 25 print('sample mean var: ', rvs.mean(), rvs.var()) 26 print('theoretical mean var', skewnorm.stats(0)) 27 rvs = skewnorm.rvs(5,size=500) 28 print('sample mean var: ', rvs.mean(), rvs.var()) 29 print('theoretical mean var', skewnorm.stats(5)) 30 print(skewnorm.cdf(1,0), stats.norm.cdf(1), skewnorm.cdf(1,0) - stats.norm.cdf(1)) 31 print(skewnorm.cdf(1,1000), stats.chi.cdf(1,1), skewnorm.cdf(1,1000) - stats.chi.cdf(1,1)) 32 print(skewnorm.sf(0.05,1000), stats.chi.sf(0.05,1), skewnorm.sf(0.05,1000) - stats.chi.sf(0.05,1)) 33 34 35def example_T(): 36 skewt = ACSkewT_gen() 37 rvs = skewt.rvs(10,0,size=500) 38 print('sample mean var: ', rvs.mean(), rvs.var()) 39 print('theoretical mean var', skewt.stats(10,0)) 40 print('t mean var', stats.t.stats(10)) 41 print(skewt.stats(10,1000)) # -> folded t distribution, as alpha -> inf 42 rvs = np.abs(stats.t.rvs(10,size=1000)) 43 print(rvs.mean(), rvs.var()) 44 45 46 47def examples_normexpand(): 48 skewnorm = SkewNorm_gen() 49 rvs = skewnorm.rvs(5,size=100) 50 normexpan = NormExpan_gen(rvs, mode='sample') 51 52 smvsk = stats.describe(rvs)[2:] 53 print('sample: mu,sig,sk,kur') 54 print(smvsk) 55 56 dmvsk = normexpan.stats(moments='mvsk') 57 print('normexpan: mu,sig,sk,kur') 58 print(dmvsk) 59 print('mvsk diff distribution - sample') 60 print(np.array(dmvsk) - np.array(smvsk)) 61 print('normexpan attributes mvsk') 62 print(mc2mvsk(normexpan.cnt)) 63 print(normexpan.mvsk) 64 65 mnc = mvsk2mnc(dmvsk) 66 mc = mnc2mc(mnc) 67 print('central moments') 68 print(mc) 69 print('non-central moments') 70 print(mnc) 71 72 73 pdffn = pdf_moments(mc) 74 print('\npdf approximation from moments') 75 print('pdf at', mc[0]-1,mc[0]+1) 76 print(pdffn([mc[0]-1,mc[0]+1])) 77 print(normexpan.pdf([mc[0]-1,mc[0]+1])) 78 79 80def examples_transf(): 81 ##lognormal = ExpTransf(a=0.0, xa=-10.0, name = 'Log transformed normal') 82 ##print(lognormal.cdf(1)) 83 ##print(stats.lognorm.cdf(1,1)) 84 ##print(lognormal.stats()) 85 ##print(stats.lognorm.stats(1)) 86 ##print(lognormal.rvs(size=10)) 87 88 print('Results for lognormal') 89 lognormalg = ExpTransf_gen(stats.norm, a=0, name = 'Log transformed normal general') 90 print(lognormalg.cdf(1)) 91 print(stats.lognorm.cdf(1,1)) 92 print(lognormalg.stats()) 93 print(stats.lognorm.stats(1)) 94 print(lognormalg.rvs(size=5)) 95 96 ##print('Results for loggamma') 97 ##loggammag = ExpTransf_gen(stats.gamma) 98 ##print(loggammag._cdf(1,10)) 99 ##print(stats.loggamma.cdf(1,10)) 100 101 print('Results for expgamma') 102 loggammaexpg = LogTransf_gen(stats.gamma) 103 print(loggammaexpg._cdf(1,10)) 104 print(stats.loggamma.cdf(1,10)) 105 print(loggammaexpg._cdf(2,15)) 106 print(stats.loggamma.cdf(2,15)) 107 108 109 # this requires change in scipy.stats.distribution 110 #print(loggammaexpg.cdf(1,10)) 111 112 print('Results for loglaplace') 113 loglaplaceg = LogTransf_gen(stats.laplace) 114 print(loglaplaceg._cdf(2)) 115 print(stats.loglaplace.cdf(2,1)) 116 loglaplaceexpg = ExpTransf_gen(stats.laplace) 117 print(loglaplaceexpg._cdf(2)) 118 stats.loglaplace.cdf(3,3) 119 #0.98148148148148151 120 loglaplaceexpg._cdf(3,0,1./3) 121 #0.98148148148148151 122 123if __name__ == '__main__': 124 example_n() 125 example_T() 126 examples_normexpand() 127 examples_transf() 128