Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Improved setup options

  • Loading branch information...
commit 8bbb777ca5a768227e8ea6c2954a44d1d6966920 1 parent 017668b
@twiecki twiecki authored
View
4 builddpkg
@@ -0,0 +1,4 @@
+#!/bin/bash
+python setup.py --command-packages=stdeb.command bdist_deb
+cd deb_dist/HDDM-0.1
+dpkg-buildpackage -rfakeroot -uc -us
View
17 examples/RatcliffFrank.conf
@@ -0,0 +1,17 @@
+[data]
+load = RatcliffFrank.csv
+save = RatcliffFrank.txt
+
+[model]
+type = lba
+is_subj_model = False
+
+[depends]
+z = da,conf
+v = conf
+
+[mcmc]
+samples = 10000
+burn = 5000
+thin = 5
+dbname = RatcliffFrank.db
View
19,995 examples/RatcliffFrank.csv
19,995 additions, 0 deletions not shown
View
513 examples/RatcliffFrank.db
513 additions, 0 deletions not shown
View
25 examples/RatcliffFrank.txt
@@ -0,0 +1,25 @@
+Mean group estimates:
+a: 1.049030
+z_('med', 'nogo'): 0.272175
+z_('med', 'go'): 1.048957
+z_('lo', 'go'): 1.048987
+v0: 0.999877
+v1: 0.000123
+t: 0.479866
+V: 1.998987
+z_('lo', 'nogo'): 0.774927
+
+Standard deviations of group parameters:
+a: 0.000305
+z_('med', 'nogo'): 0.012303
+z_('med', 'go'): 0.000297
+z_('lo', 'go'): 0.000303
+v0: 0.000149
+v1: 0.000149
+t: 0.000068
+V: 0.000992
+z_('lo', 'nogo'): 0.005602
+
+General model stats:
+logp: -1843.656935
+dic: 3693.770048
View
19 examples/RatcliffFrankTer.txt
@@ -0,0 +1,19 @@
+Mean group estimates:
+a: 1.135953
+t_('lo', 'go'): 0.439632
+t_('lo', 'nogo'): 0.550229
+t_('med', 'go'): 0.443847
+t_('med', 'nogo'): 0.615013
+v: 2.495221
+
+Standard deviations of group parameters:
+a: 0.004984
+t_('lo', 'go'): 0.000836
+t_('lo', 'nogo'): 0.000741
+t_('med', 'go'): 0.000483
+t_('med', 'nogo'): 0.000329
+v: 0.018246
+
+General model stats:
+logp: 12126.898257
+dic: -24249.478863
View
BIN  examples/V.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/a.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/a_('lo',).png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/a_('med',).png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/a_lo.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/a_med.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
14 examples/example.conf
@@ -0,0 +1,14 @@
+[data]
+load = data/example_subj.csv
+save = data/example_subj_out.txt
+
+[model]
+type = simple
+is_subj_model = True
+debug = True
+
+[depends]
+
+[mcmc]
+samples=500
+burn=100
View
0  examples/example.csv
No changes.
View
BIN  examples/t_('lo', 'go').png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/t_('lo', 'nogo').png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/t_('med', 'go').png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/t_('med', 'nogo').png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/t_go.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/t_nogo.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/v.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/v0.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/v1.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
BIN  examples/z.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
View
19 setup.py
@@ -0,0 +1,19 @@
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Distutils import build_ext
+
+
+setup(
+ name="HDDM",
+ version="0.1",
+ author="Thomas V. Wiecki",
+ author_email="thomas_wiecki@brown.edu",
+ url="http://code.google.com/p/hddm",
+ packages=["hddm"],
+ description="HDDM is a python module that implements Hierarchical Bayesian estimation of Drift Diffusion Models.",
+ requires=['NumPy (>=1.3.0)', 'PyMC (>=2.0)', 'Cython'],
+ cmdclass = {'build_ext': build_ext},
+ ext_modules = [Extension("wfpt", ["src/wfpt.pyx"]),
+ Extension("lba", ["src/lba.pyx"])]
+)
+
View
7 setupegg.py
@@ -0,0 +1,7 @@
+#!/usr/bin/env python
+"""
+A setup.py script to use setuptools, which gives egg goodness, etc.
+"""
+
+from setuptools import setup
+execfile('setup.py')
View
151 src/lba64.pyx
@@ -1,151 +0,0 @@
-#!/usr/bin/python
-#
-# Cython version of the Navarro & Fuss, 2009 DDM PDF. Based directly
-# on the following code by Navarro & Fuss:
-# http://www.psychocmath.logy.adelaide.edu.au/personalpages/staff/danielnavarro/resources/wfpt.m
-#
-# This implementation is about 170 times faster than the matlab
-# reference version.
-#
-# Copyleft Thomas Wiecki (thomas_wiecki[at]brown.edu), 2010
-# GPLv3
-from __future__ import division
-from copy import copy
-import numpy as np
-cimport numpy as cnp
-
-import scipy as sp
-import scipy.stats
-from scipy.stats.distributions import norm
-
-cimport cython
-
-# Define data type
-DTYPE = np.double
-ctypedef cnp.double_t DTYPE_t
-cdef double PI = 3.1415926535897
-
-cdef extern from "math.h":
- double sin(double)
- double log(double)
- double exp(double)
- double sqrt(double)
- double fmax(double, double)
- double pow(double, double)
- double ceil(double)
- double floor(double)
- double fabs(double)
- double erf(double)
-
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def fptcdf(cnp.ndarray[DTYPE_t, ndim=1] t, double z, double a, double driftrate, double sddrift):
- cdef cnp.ndarray[DTYPE_t, ndim=1] zs=t*sddrift
- cdef cnp.ndarray[DTYPE_t, ndim=1] zu=t*driftrate
- cdef cnp.ndarray[DTYPE_t, ndim=1] aminuszu=a-zu
- cdef cnp.ndarray[DTYPE_t, ndim=1] xx=aminuszu-z
- cdef cnp.ndarray[DTYPE_t, ndim=1] azu=aminuszu/zs
- cdef cnp.ndarray[DTYPE_t, ndim=1] azumax=xx/zs
- cdef cnp.ndarray[DTYPE_t, ndim=1] tmp1=zs*(norm.pdf(azumax)-norm.pdf(azu))
- cdef cnp.ndarray[DTYPE_t, ndim=1] tmp2=xx*norm.cdf(azumax)-aminuszu*norm.cdf(azu)
- return 1+(tmp1+tmp2)/z
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def fptpdf(cnp.ndarray[DTYPE_t, ndim=1] t, double z, double a, double driftrate, double sddrift):
- cdef cnp.ndarray[DTYPE_t, ndim=1] zs=t*sddrift
- cdef cnp.ndarray[DTYPE_t, ndim=1] zu=t*driftrate
- cdef cnp.ndarray[DTYPE_t, ndim=1] aminuszu=a-zu
- cdef cnp.ndarray[DTYPE_t, ndim=1] azu=aminuszu/zs
- cdef cnp.ndarray[DTYPE_t, ndim=1] azumax=(aminuszu-z)/zs
-
- return (driftrate*(norm.cdf(azu)-norm.cdf(azumax)) + sddrift*(norm.pdf(azumax)-norm.pdf(azu)))/z
-
-def fptpdf2(cnp.ndarray[DTYPE_t, ndim=1] t, double z, double a, double driftrate, double sddrift):
- return driftrate*(norm.cdf((a-(t*driftrate))/(t*sddrift))-norm.cdf(a-(t*driftrate)-z/(t*sddrift))) + \
- sddrift*(norm.pdf(a-(t*driftrate)-z/(t*sddrift))-norm.pdf(a-((t*driftrate)/(t*sddrift))))/z
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def lba(cnp.ndarray[DTYPE_t, ndim=1] t, double z, double a, cnp.ndarray[DTYPE_t, ndim=1] drift, double sv):
- # Generates defective PDF for responses on node #1.
- cdef int i=0
- cdef int N=drift.shape[0] # Number of responses.
- #cdef np.ndarray[DTYPE_t, ndim=1] G = np.empty_like(t)
- #cdef np.ndarray[DTYPE_t, ndim=2] tmp = np.empty((t.shape[0], N-1), dtype=DTYPE)
- cdef cnp.ndarray[DTYPE_t, ndim=1] G=1-fptcdf(t=t,z=z,a=a,driftrate=drift[1],sddrift=sv)
- #if (N>2):
- #for i from 1 <= i <= N:
- #tmp[:,i]=fptcdf(t=t,z=z,a=a,driftrate=drift[i],sddrift=sv)
- #G=apply(1-tmp,1,prod)
- #np.prod(1-tmp, axis=0, out=G)
- # print ''
- #else:
- # cdef np.ndarray[DTYPE_t, ndim=1] G=1-fptcdf(t=t,z=z,a=a,driftrate=drift[1],sddrift=sv)
-
- return G*fptpdf(t=t,z=z,a=a,driftrate=drift[0],sddrift=sv)
-
-# Second approach
-cdef inline double norm_pdf(double x): return 0.3989422804014327 * exp(-(x**2/2))
-cdef inline double norm_cdf(double x): return .5 * (1 + erf((x)/1.4142135623730951))
-
-cdef DTYPE_t fptcdf_single(DTYPE_t t, double z, double a, double driftrate, double sddrift):
- cdef double zs=t*sddrift
- cdef double zu=t*driftrate
- cdef double aminuszu=a-zu
- cdef double xx=aminuszu-z
- cdef double azu=aminuszu/zs
- cdef double azumax=xx/zs
- cdef double tmp1 = zs * (norm_pdf(azumax) - norm_pdf(azu))
- cdef double tmp2 = (xx * norm_cdf(azumax)) - (aminuszu * norm_cdf(azu))
- return 1+(tmp1+tmp2)/z
-
-cdef DTYPE_t fptpdf_single(DTYPE_t t, double z, double a, double driftrate, double sddrift):
- cdef double zs=t*sddrift
- cdef double zu=t*driftrate
- cdef double aminuszu=a-zu
- cdef double azu=aminuszu/zs
- cdef double azumax=(aminuszu-z)/zs
-
- return (driftrate*(norm_cdf(azu)-norm_cdf(azumax)) + sddrift*(norm_pdf(azumax)-norm_pdf(azu)))/z
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def lba_single(cnp.ndarray[DTYPE_t, ndim=1] t, double z, double a, cnp.ndarray[DTYPE_t, ndim=1] drift, double sv, cnp.ndarray[DTYPE_t, ndim=1] out=None):
- # Generates defective PDF for responses on node #1.
- cdef unsigned int i=0
- cdef unsigned int max = <unsigned int> t.shape[0]
- if out is None:
- out = np.empty(max, dtype=DTYPE)
- elif out.shape[0] != t.shape[0]:
- raise ValueError('out array must have the same shape as input array t')
-
-
- for i from 0 <= i < max:
- out[i] = (1-fptcdf_single(t[i], z, a, drift[<unsigned int> 1], sv)) * \
- fptpdf_single(t[i], z, a, drift[<unsigned int> 0], sv)
-
- return out
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def lba_like(cnp.ndarray[DTYPE_t, ndim=1] value, cnp.ndarray[DTYPE_t, ndim=1] resps, double z, double a, double ter, cnp.ndarray[DTYPE_t, ndim=1] drift, double sv, unsigned int logp=0):
- # Rescale parameters so as to fit in the same range as the DDM
- cdef cnp.ndarray[DTYPE_t, ndim=1] rt = (np.abs(value) - ter)
- cdef unsigned int nresp = value.shape[0]
- cdef unsigned int i
- cdef cnp.ndarray[DTYPE_t, ndim=1] probs = np.empty(nresp, dtype=DTYPE)
-
- if a <= z:
- #print "Starting point larger than threshold!"
- return -np.Inf
-
- #print "z: %f, a: %f, ter: %i, v: %f, sv: %f" % (z, a, ter, v[0], sv)
- for i from 0 <= i < nresp:
- if resps[i] == 1:
- probs[i] = (1-fptcdf_single(rt[i], z, a, drift[1], sv)) * \
- fptpdf_single(rt[i], z, a, drift[0], sv)
- else:
- probs[i] = (1-fptcdf_single(rt[i], z, a, drift[0], sv)) * \
- fptpdf_single(rt[i], z, a, drift[1], sv)
-
- if logp == 1:
- return np.sum(np.log(probs))
- else:
- return probs
View
158 src/wfpt64.pyx
@@ -1,158 +0,0 @@
-#!/usr/bin/python
-#
-# Cython version of the Navarro & Fuss, 2009 DDM PDF. Based directly
-# on the following code by Navarro & Fuss:
-# http://www.psychocmath.logy.adelaide.edu.au/personalpages/staff/danielnavarro/resources/wfpt.m
-#
-# This implementation is about 170 times faster than the matlab
-# reference version.
-#
-# Copyleft Thomas Wiecki (thomas_wiecki[at]brown.edu), 2010
-# GPLv3
-from __future__ import division
-from copy import copy
-import numpy as np
-cimport numpy as np
-
-cimport cython
-
-cdef extern from "math.h":
- double sin(double)
- double log(double)
- double exp(double)
- double sqrt(double)
- double fmax(double, double)
- double pow(double, double)
- double ceil(double)
- double floor(double)
- double fabs(double)
-
-# Define data type
-DTYPE = np.double
-ctypedef double DTYPE_t
-
-cpdef double pdf(double t, double v, double a, double z, double err, unsigned int logp=0):
- """Compute the likelihood of the drift diffusion model using the method
- and implementation of Navarro & Fuss, 2009.
- """
- if t <= 0:
- if logp == 0:
- return 0
- else:
- return -np.Inf
- cdef double tt = t/(pow(a,2)) # use normalized time
- cdef double w = z/a # convert to relative start point
-
- cdef double kl, ks, p
- cdef double PI = 3.1415926535897
- cdef double PIs = 9.869604401089358 # PI^2
- cdef int k, K, lower, upper
-
- # calculate number of terms needed for large t
- if PI*tt*err<1: # if error threshold is set low enough
- kl=sqrt(-2*log(PI*tt*err)/(PIs*tt)) # bound
- kl=fmax(kl,1./(PI*sqrt(tt))) # ensure boundary conditions met
- else: # if error threshold set too high
- kl=1./(PI*sqrt(tt)) # set to boundary condition
-
-
- # calculate number of terms needed for small t
- if 2*sqrt(2*PI*tt)*err<1: # if error threshold is set low enough
- ks=2+sqrt(-2*tt*log(2*sqrt(2*PI*tt)*err)) # bound
- ks=fmax(ks,sqrt(tt)+1) # ensure boundary conditions are met
- else: # if error threshold was set too high
- ks=2 # minimal kappa for that case
-
- # compute f(tt|0,1,w)
- p=0 #initialize density
- if ks<kl: # if small t is better (i.e., lambda<0)
- K=<int>(ceil(ks)) # round to smallest integer meeting error
- lower = <int>(-floor((K-1)/2.))
- upper = <int>(ceil((K-1)/2.))
- for k from lower <= k <= upper: # loop over k
- p=p+(w+2*k)*exp(-(pow((w+2*k),2))/2/tt) # increment sum
- p=p/sqrt(2*PI*pow(tt,3)) # add constant term
-
- else: # if large t is better...
- K=<int>(ceil(kl)) # round to smallest integer meeting error
- for k from 1 <= k <= K:
- p=p+k*exp(-(pow(k,2))*(PIs)*tt/2)*sin(k*PI*w) # increment sum
- p=p*PI # add constant term
-
- # convert to f(t|v,a,w)
- if logp == 0:
- return p*exp(-v*a*w -(pow(v,2))*t/2.)/(pow(a,2))
- else:
- return log(p) + (-v*a*w -(pow(v,2))*t/2.) - 2*log(a)
-
-cpdef double pdf_sign(double t, double v, double a, double z, double ter, double err, int logp=0):
- """Wiener likelihood function for two response types. Lower bound
- responses have negative t, upper boundary response have positive t"""
- if a<z or z<=0 or z<0 or a<0:
- return -np.Inf
-
- if t<0:
- return pdf(fabs(t)-ter, v, a, z, err, logp)
- else:
- return pdf(t-ter, -v, a, a-z, err, logp)
-
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def pdf_array(np.ndarray[DTYPE_t, ndim=1] x, double v, double a, double z, double ter, double err, int logp=0):
- cdef unsigned int size = x.shape[0]
- cdef unsigned int i
- cdef np.ndarray[DTYPE_t, ndim=1] y = np.empty(size, dtype=DTYPE)
- for i from 0 <= i < size:
- y[i] = pdf_sign(x[i], v, a, z, ter, err, logp)
- return y
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def pdf_array_multi(np.ndarray[DTYPE_t, ndim=1] x, v, a, z, ter, double err, int logp=0, multi=None):
- cdef unsigned int size = x.shape[0]
- cdef unsigned int i
- cdef np.ndarray[DTYPE_t, ndim=1] y = np.empty(size, dtype=DTYPE)
-
- if multi is None:
- return pdf_array(x, v=v, a=a, z=z, ter=ter, err=err, logp=logp)
- else:
- params = {'v':v, 'z':z, 'ter':ter, 'a':a}
- params_iter = copy(params)
- for i from 0 <= i < size:
- for param in multi:
- params_iter[param] = params[param][i]
-
- y[i] = pdf_sign(x[i], params_iter['v'], params_iter['a'], params_iter['z'], params_iter['ter'], err=err, logp=logp)
-
- return y
-
-@cython.boundscheck(False) # turn of bounds-checking for entire function
-def wiener_like_full_avg(np.ndarray[DTYPE_t, ndim=1] t, double v, double sv, z, double sz, double ter, double ster, a, double err=.0001, int logp=0, unsigned int reps=10, a_is_multi=False):
- cdef unsigned int num_resps = t.shape[0]
- cdef unsigned int rep, i
-
- if logp == 1:
- zero_prob = -np.Inf
- else:
- zero_prob = 0
-
- # Create samples
- cdef np.ndarray[DTYPE_t, ndim=2] ter_samples = np.random.uniform(size=(reps, num_resps), low=ter-ster/2., high=ter+ster/2.)
- cdef np.ndarray[DTYPE_t, ndim=2] z_samples = np.random.uniform(size=(reps, num_resps), low=z-sz/2., high=z+sz/2.)
- cdef np.ndarray[DTYPE_t, ndim=2] v_samples = np.random.normal(size=(reps, num_resps), loc=v, scale=sv)
- cdef np.ndarray[DTYPE_t, ndim=2] probs = np.empty((reps, num_resps), dtype=DTYPE)
-
- for rep from 0 <= rep < reps:
- for i from 0 <= i < num_resps:
- if a_is_multi:
- a_val = a[i]
- else:
- a_val = a
- if (fabs(t[i])-ter_samples[rep,i]) < 0:
- probs[rep,i] = zero_prob
- elif a_val <= z_samples[rep,i]:
- probs[rep,i] = zero_prob
- else:
- probs[rep,i] = pdf_sign(t[i], v_samples[rep,i], a_val, z_samples[rep,i], ter_samples[rep,i], err=err, logp=logp)
-
- return np.mean(probs, axis=0)
-
View
3  stdeb.cfg
@@ -0,0 +1,3 @@
+[HDDM]
+Depends: python (>= 2.5.0), python-numpy (>= 1.1.1), python-scipy, python-matplotlib, python-traits, python-traitsgui
+Build-Depends: cython
Please sign in to comment.
Something went wrong with that request. Please try again.