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