In [44]:

    

import itertools

import matplotlib.pyplot as plt
import matplotlib as mpl
from pymc3 import Model, Normal, Slice
from pymc3 import sample
from pymc3 import traceplot
from pymc3.distributions import Interpolated
from theano import as_op
import theano.tensor as tt
import numpy as np
from scipy import stats
import tqdm

%matplotlib inline

%load_ext version_information

%version_information pymc3, scipy


    

    




The version_information extension is already loaded. To reload it, use:
  %reload_ext version_information






    
Out[44]:




Software
Version
Python
3.6.2 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]
IPython
6.1.0
OS
Darwin 15.6.0 x86_64 i386 64bit
pymc3
3.1
scipy
0.19.1
Mon Sep 18 12:13:55 2017 MDT

In [3]:

    

stats.kappa3?


    

In [14]:

    

k3_rvs = stats.kappa3.rvs(a=1, loc=90, scale=45, size=10000)
plt.hist(k3_rvs, 50);


    

In [111]:

    

# stats.kappa4?
h=100
k=0
k4_pdf = stats.kappa4.rvs(h=k, k=k, loc=90, scale=20, size=10000)
plt.hist(k4_pdf, 50);
mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
print(mean, var, skew, kurt)


    

    




4.621500393020947 0.023720538535741298 17.440372233406986 431.65646930084574






    

In [24]:

    

import scipy.optimize


    

In [92]:

    

def moments(x):
    h, k = x
    mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
    return skew

minme = lambda x: np.abs(moments(x))

x0 = [.2, 4]

minme(x0)


    

    
Out[92]:





7.5529025101903082

In [97]:

    

# optimization seems slow, lets try a brute force to get a feel
N=7
h = np.linspace(0, 5, N)
k = np.linspace(0, 5, N)

H, K = np.meshgrid(h, k)

ans = np.zeros((N,N), dtype=float)


    

In [98]:

    

for i, j in tqdm.tqdm(itertools.product(range(N), range(N)), total=N*N):
    ans[i,j] = minme((H[i,j], K[i,j]))


    

    




 86%|████████▌ | 42/49 [00:52<00:02,  2.66it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.
  warnings.warn(msg, IntegrationWarning)
100%|██████████| 49/49 [00:58<00:00,  2.84it/s]

In [100]:

    

plt.contourf(H,K,ans)
plt.colorbar()
plt.contour(H,K,ans, levels=[0,.1,.2], colors='w')

print(np.unravel_index(np.argmin(ans), ans.shape))
print(ans.min())
min = np.unravel_index(np.argmin(ans), ans.shape)
plt.scatter(H[min], K[min], marker='o', c='r')
print(H[min], K[min])


    

    




(1, 1)
0.00946777071186
0.833333333333 0.833333333333






    

In [91]:

    

niter = 0
def callbackF(x):
    global niter
    print("{} : {}, {}".format(niter, x[0], x[1]))
    niter += 1

  
scipy.optimize.minimize(minme, x0, method='TNC', callback=callbackF, bounds=[[0,1], [0,1]])


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-91-6a675168301b> in <module>()
      6 
      7 
----> 8 scipy.optimize.minimize(minme, x0, method='TNC', callback=callbackF, bounds=[[0,1], [0,1]])

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    451     elif meth == 'tnc':
    452         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
--> 453                              **options)
    454     elif meth == 'cobyla':
    455         return _minimize_cobyla(fun, x0, args, constraints, **options)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/tnc.py in _minimize_tnc(fun, x0, args, jac, bounds, eps, scale, offset, mesg_num, maxCGit, maxiter, eta, stepmx, accuracy, minfev, ftol, xtol, gtol, rescale, disp, callback, **unknown_options)
    407                                         offset, messages, maxCGit, maxfun,
    408                                         eta, stepmx, accuracy, fmin, ftol,
--> 409                                         xtol, pgtol, rescale, callback)
    410 
    411     funv, jacv = func_and_grad(x)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/tnc.py in func_and_grad(x)
    365         def func_and_grad(x):
    366             f = fun(x, *args)
--> 367             g = approx_fprime(x, fun, epsilon, *args)
    368             return f, g
    369     else:

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/optimize.py in approx_fprime(xk, f, epsilon, *args)
    686 
    687     """
--> 688     return _approx_fprime_helper(xk, f, epsilon, args=args)
    689 
    690 

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
    626         ei[k] = 1.0
    627         d = epsilon * ei
--> 628         grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
    629         ei[k] = 0.0
    630     return grad

<ipython-input-31-62019fb4b4c9> in <lambda>(x)
      4     return skew
      5 
----> 6 minme = lambda x: np.abs(moments(x))
      7 
      8 x0 = [.2, 4]

<ipython-input-31-62019fb4b4c9> in moments(x)
      1 def moments(x):
      2     h, k = x
----> 3     mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
      4     return skew
      5 

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in stats(self, *args, **kwds)
   1025             if 'm' in moments:
   1026                 if mu is None:
-> 1027                     mu = self._munp(1, *goodargs)
   1028                 out0 = default.copy()
   1029                 place(out0, cond, mu * scale + loc)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _munp(self, n, *args)
    783         # Silence floating point warnings from integration.
    784         olderr = np.seterr(all='ignore')
--> 785         vals = self.generic_moment(n, *args)
    786         np.seterr(**olderr)
    787         return vals

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
   2732             vargs.extend([kwargs[_n] for _n in names])
   2733 
-> 2734         return self._vectorize_call(func=func, args=vargs)
   2735 
   2736     def _get_ufunc_and_otypes(self, func, args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
   2808                       for a in args]
   2809 
-> 2810             outputs = ufunc(*inputs)
   2811 
   2812             if ufunc.nout == 1:

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _mom1_sc(self, m, *args)
   1603 
   1604     def _mom1_sc(self, m, *args):
-> 1605         return integrate.quad(self._mom_integ1, 0, 1, args=(m,)+args)[0]
   1606 
   1607     def _pdf(self, x, *args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    321     if (weight is None):
    322         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323                        points)
    324     else:
    325         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    386     if points is None:
    387         if infbounds == 0:
--> 388             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    389         else:
    390             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _mom_integ1(self, q, m, *args)
   1600     # moment calculated using ppf
   1601     def _mom_integ1(self, q, m, *args):
-> 1602         return (self.ppf(q, *args))**m
   1603 
   1604     def _mom1_sc(self, m, *args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in ppf(self, q, *args, **kwds)
   1898         q, loc, scale = map(asarray, (q, loc, scale))
   1899         args = tuple(map(asarray, args))
-> 1900         cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
   1901         cond1 = (0 < q) & (q < 1)
   1902         cond2 = cond0 & (q == 0)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py in _argcheck(self, h, k)
   3496                              [f0, f1, f1, f0, f1, f1],
   3497                              [h, k],
-> 3498                              default=np.nan)
   3499         return h == h
   3500 

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/_lib/_util.py in _lazyselect(condlist, choicelist, arrays, default)
     95             continue
     96         cond, _ = np.broadcast_arrays(cond, arrays[0])
---> 97         temp = tuple(np.extract(cond, arr) for arr in arrays)
     98         np.place(out, cond, func(*temp))
     99     return out

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/_lib/_util.py in <genexpr>(.0)
     95             continue
     96         cond, _ = np.broadcast_arrays(cond, arrays[0])
---> 97         temp = tuple(np.extract(cond, arr) for arr in arrays)
     98         np.place(out, cond, func(*temp))
     99     return out

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in extract(condition, arr)
   2324 
   2325     """
-> 2326     return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
   2327 
   2328 

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in take(a, indices, axis, out, mode)
    132            [5, 7]])
    133     """
--> 134     return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)
    135 
    136 

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in _wrapfunc(obj, method, *args, **kwds)
     55 def _wrapfunc(obj, method, *args, **kwds):
     56     try:
---> 57         return getattr(obj, method)(*args, **kwds)
     58 
     59     # An AttributeError occurs if the object does not have

KeyboardInterrupt: 

In [112]:

    

kappa = 3.99390425811
mean, var, skew, kurt = stats.vonmises.stats(kappa, moments='mvsk')
print(mean, var, skew, kurt)


    

    




4.995631231732412e-17 0.2988042418167408 4.572598569626176e-16 0.6905715325848036

In [122]:

    

x = np.linspace(stats.vonmises.ppf(0.01, kappa), stats.vonmises.ppf(0.99, kappa), 100)

plt.plot(x, stats.vonmises.pdf(x, kappa), 'r-', lw=5, alpha=0.6, label='vonmises pdf')

x = np.linspace(stats.norm.ppf(0.01, loc=0, scale=.4, ), stats.norm.ppf(0.99, loc=0, scale=.4, ), 100)

plt.plot(x, stats.norm.pdf(x, loc=0, scale=.4), 'k', lw=2, alpha=0.6, label='normal pdf')


    

    
Out[122]:





[<matplotlib.lines.Line2D at 0x1228fb240>]






    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: