# -*- coding: utf-8 -*-
"""
Created on Sat Mar 13 07:56:22 2010
Author: josef-pktd
"""
from __future__ import print_function
import sympy as sy
[docs]def pdf(x, mu, sigma):
"""Return the probability density function as an expression in x"""
#x = sy.sympify(x)
return 1/(sigma*sy.sqrt(2*sy.pi)) * sy.exp(-(x-mu)**2 / (2*sigma**2))
[docs]def cdf(x, mu, sigma):
"""Return the cumulative density function as an expression in x"""
#x = sy.sympify(x)
return (1+sy.erf((x-mu)/(sigma*sy.sqrt(2))))/2
mu = sy.Symbol('mu')
sigma = sy.Symbol('sigma')
sigma2 = sy.Symbol('sigma2')
x = sy.Symbol('x')
y = sy.Symbol('y')
df = sy.Symbol('df')
s = sy.Symbol('s')
dldxnorm = sy.log(pdf(x, mu,sigma)).diff(x)
print(sy.simplify(dldxnorm))
print(sy.diff(sy.log(sy.gamma((s+1)/2)),s))
print(sy.diff((df+1)/2. * sy.log(1+df/(df-2)), df))
#standard t distribution, not verified
tllf1 = sy.log(sy.gamma((df+1)/2.)) - sy.log(sy.gamma(df/2.)) - 0.5*sy.log((df)*sy.pi)
tllf2 = (df+1.)/2. * sy.log(1. + (y-mu)**2/(df)/sigma2) + 0.5 * sy.log(sigma2)
tllf2std = (df+1.)/2. * sy.log(1. + y**2/df) + 0.5
tllf = tllf1 - tllf2
print(tllf1.diff(df))
print(tllf2.diff(y))
dlddf = (tllf1-tllf2).diff(df)
print(dlddf)
print(sy.cse(dlddf))
print('\n derivative of loglike of t distribution wrt df')
for k,v in sy.cse(dlddf)[0]: print(k,'=',v)
print(sy.cse(dlddf)[1][0])
print('\nstandard t distribution, dll_df, dll_dy')
tllfstd = tllf1 - tllf2std
print(tllfstd.diff(df))
print(tllfstd.diff(y))
print('\n')
print(dlddf.subs(dict(y=1,mu=1,sigma2=1.5,df=10.0001)))
print(dlddf.subs(dict(y=1,mu=1,sigma2=1.5,df=10.0001)).evalf())
# Note: derivatives of nested function doesn't work in sympy
# at least not higher order derivatives (second or larger)
# looks like print(failure