diff --git a/g6pd/__init__.py b/g6pd/__init__.py index dee2d5b..2360059 100644 --- a/g6pd/__init__.py +++ b/g6pd/__init__.py @@ -27,25 +27,25 @@ def check_data(input): if np.any(input[pos]<0) or np.any(input[pos]>input[n]): raise ValueError, 'Some observations are negative.' -def male_def(sp_sub, ceiling): +def male_def(sp_sub): allele = sp_sub.copy('F') - allele = invlogit(allele)*ceiling + allele = invlogit(allele) return allele -def fem_def_conservative(sp_sub, ceiling): - hom = male_def(sp_sub, ceiling) +def fem_def_conservative(sp_sub): + hom = male_def(sp_sub) fast_inplace_mul(hom,hom) return hom -def hw_hetero(sp_sub, ceiling): - p = male_def(sp_sub, ceiling) +def hw_hetero(sp_sub): + p = male_def(sp_sub) q = fast_inplace_scalar_add(-p,1) fast_inplace_mul(p,q) return 2*p -def fem_def(sp_sub, ceiling, a, b): - homo = male_def(sp_sub, ceiling) - hetero = hw_hetero(sp_sub, ceiling) +def fem_def(sp_sub, a, b): + homo = male_def(sp_sub) + hetero = hw_hetero(sp_sub) het_def = pm.rbeta(a,b) hetero *= het_def return hetero+homo @@ -55,15 +55,15 @@ def fem_def(sp_sub, ceiling, a, b): def validate_male(data): obs = data.male_pos n = data.n_male - def f(sp_sub, ceiling, n=n): - return pm.rbinomial(n=n,p=pm.invlogit(sp_sub)*ceiling) + def f(sp_sub, n=n): + return pm.rbinomial(n=n,p=pm.invlogit(sp_sub)) return obs, n, f def validate_female(data): obs = data.fem_pos n = data.n_fem - def f(sp_sub, ceiling, a, b, n=n): - p = pm.invlogit(sp_sub)*ceiling + def f(sp_sub, a, b, n=n): + p = pm.invlogit(sp_sub) h = pm.rbeta(a,b,size=len(sp_sub)) p_def = g6pd.p_fem_def(p,h) return pm.rbinomial(n=n, p=p) @@ -81,8 +81,8 @@ def area_male(gc): def h(**region): return np.array(region.values()[0]) - def f(sp_sub, x, ceiling): - p = pm.invlogit(sp_sub(x))*ceiling + def f(sp_sub, x): + p = pm.invlogit(sp_sub(x)) return p g = {gc.keys()[0]: f} @@ -96,8 +96,8 @@ def area_fem_def_conservative(gc): def h(**region): return np.array(region.values()[0]) - def f(sp_sub, x, ceiling): - p = pm.invlogit(sp_sub(x))*ceiling + def f(sp_sub, x): + p = pm.invlogit(sp_sub(x)) return p**2 g = {gc.keys()[0]: f} @@ -111,8 +111,8 @@ def area_fem_def(gc): def h(**region): return np.array(region.values()[0]) - def f(sp_sub, x, ceiling, a, b): - p = pm.invlogit(sp_sub(x))*ceiling + def f(sp_sub, x, a, b): + p = pm.invlogit(sp_sub(x)) h = pm.rbeta(a,b,size=len(p)) return g6pd.p_fem_def(p,h)