diff --git a/solution1.py b/solution1.py new file mode 100644 index 0000000..3c1a09f --- /dev/null +++ b/solution1.py @@ -0,0 +1,75 @@ +import numpy +import pandas +from scipy.optimize import minimize +from scipy.stats import norm +from scipy.stats import chi2 +from plotnine import * + +def linear(p, obs): + B0 = p[0] + B1 = p[1] + sigma = p[2] + + expected = B0 + (B1 * obs.mutation) + nll = -1 * norm(expected, sigma).logpdf(obs.ponzr1Counts).sum() + return nll + +def nullH(p, obs): + B0 = p[0] + sigma = p[1] + + + expected = B0 + nll = -1 * norm(expected, sigma).logpdf(obs.ponzr1Counts).sum() + return nll + + +data = pandas.read_csv('ponzr1.csv', header = 0, sep = ',') +data['mutation'] = data["mutation"].map({'WT' : 0, 'M124K' : 1, 'V456D' : 2, 'I213N' : 3}) + +subset=data.loc[data.mutation.isin(['0','1']),:] + +initialLinGuess = numpy.array([1,1,1]) +linfit = minimize(linear, initialLinGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +initialNullGuess = numpy.array([1,1]) +nullfit = minimize(nullH, initialNullGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +linfit.fun +nullfit.fun + +D=2*(nullfit.fun - linfit.fun) +print("M124K P value:") +print(1-chi2.cdf(x=D, df=1)) + +subset=data.loc[data.mutation.isin(['0','2']),:] + +initialLinGuess = numpy.array([1,1,1]) +linfit = minimize(linear, initialLinGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +initialNullGuess = numpy.array([1,1]) +nullfit = minimize(nullH, initialNullGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +linfit.fun +nullfit.fun + +D=2*(nullfit.fun - linfit.fun) +print("V456D P value:") +print(1-chi2.cdf(x=D, df=1)) + +subset=data.loc[data.mutation.isin(['0','3']),:] + +initialLinGuess = numpy.array([1,1,1]) +linfit = minimize(linear, initialLinGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +initialNullGuess = numpy.array([1,1]) +nullfit = minimize(nullH, initialNullGuess, method="Nelder-Mead", options={'disp': True}, args=subset) + +linfit.fun +nullfit.fun + +D=2*(nullfit.fun - linfit.fun) +print("I213N P value:") +print(1-chi2.cdf(x=D, df=1)) + + diff --git a/solution2.py b/solution2.py new file mode 100644 index 0000000..6601d28 --- /dev/null +++ b/solution2.py @@ -0,0 +1,22 @@ +import numpy +import pandas +from scipy.optimize import minimize +from scipy.stats import norm +from scipy.stats import chi2 +from plotnine import * + +def monod(p, obs): + UMAX = p[0] + KS = p[1] + sigma = p[2] + + expected = UMAX*(obs.S / (obs.S + KS)) + nll = -1 * norm(expected, sigma).logpdf(obs.u).sum() + return nll + +data = pandas.read_csv('MmarinumGrowth.csv', header = 0, sep = ',') + +initialGuess = numpy.array([1,1,1]) +fit = minimize(monod, initialGuess, method="Nelder-Mead", options={'disp': True}, args=data) + +print(fit.x) \ No newline at end of file diff --git a/solution3.py b/solution3.py new file mode 100644 index 0000000..d6507d5 --- /dev/null +++ b/solution3.py @@ -0,0 +1,53 @@ +import numpy +import pandas +from scipy.optimize import minimize +from scipy.stats import norm +from scipy.stats import chi2 +from plotnine import * + +def const(p, obs): + a = p[0] + sigma = p[1] + + expected = a + nll = -1 * norm(expected, sigma).logpdf(obs.decomp).sum() + return nll + +def linear(p, obs): + a = p[0] + b = p[1] + sigma = p[2] + + expected = a + b * obs.Ms + nll = -1 * norm(expected, sigma).logpdf(obs.decomp).sum() + return nll + +def hump(p, obs): + a = p[0] + b = p[1] + c = p[2] + sigma = p[3] + + expected = a + b*obs.Ms + c*(obs.Ms) * (obs.Ms) + nll = -1 * norm(expected, sigma).logpdf(obs.decomp).sum() + return nll + +data = pandas.read_csv('leafDecomp.csv', header = 0, sep = ',') + +initialGuess = numpy.array([1,1,1,1]) +constfit = minimize(const, initialGuess, method="Nelder-Mead", options={'disp': True}, args=data) +linfit = minimize(linear, initialGuess, method="Nelder-Mead", options={'disp': True}, args=data) +humpfit = minimize(hump, initialGuess, method="Nelder-Mead", options={'disp': True}, args=data) + + +D=2*(constfit.fun - linfit.fun) +print("linear over const fit") +print(1-chi2.cdf(x=D, df=1)) + +D=2*(humpfit.fun - linfit.fun) +print("linear over hump fit") +print(1-chi2.cdf(x=D, df=1)) + +D=2*(constfit.fun - humpfit.fun) +print("hump over const fit") +print(1-chi2.cdf(x=D, df=2))