Working version to fit data

three_parameter_screening.py

@@ -173,26 +173,73 @@ def likelihood_function(D, I, n, s, r, alpha, beta):
     float
         The likelihood value.
     """
-    L = 1.0
+    L = []
     for i in range(len(D)):
-        s=numpy.sum(s[i])
-        term1 = D[i] ** s
+        qs=np.sum(s[i])
+        term1 = D[i] ** qs
         term2 = I[i] ** r[i]
-        term3 = (1 - D[i] - I[i]) ** (n[i] - s - r[i])
-        product_term = np.prod([(alpha[j] / beta) ** s[i][j] for j in range(len(alpha))])
+        term3 = (1 - D[i] - I[i]) ** (n[i] - qs - r[i])
+        product_term = np.prod((alpha / beta) ** s[i] )
         L_i = term1 * term2 * term3 * product_term
-        L *= L_i
+        L.append( L_i )
     return L
 
-def likelihood_function(pars,n,s,r):
+
+def mylog(x):
+   if x==0:
+      return - np.inf
+   return np.log(x)
+
+
+def loglikelihood_function(D, I, n, s, r, alpha, beta):
+# Same as likelihood above
+    L = []
+    for i in range(len(D)):
+        qs=np.sum(s[i])
+        term1 = qs*mylog(D[i])
+        term2=r[i]*mylog(I[i])
+        term3 = (n[i]-qs-r[i])*mylog(1 - D[i] - I[i]) 
+        product_term=0
+        if beta!=0:
+           product_term = np.sum(s[i]*mylog(alpha / beta))
+        L_i = term1 + term2 + term3 + product_term
+        L.append( L_i )
+    return L
+
+
+
+
+def likelihood_function1(pars,n,s,r):
+   w=pars[0]
+   mu=pars[1]
+   beta=pars[2]
    def Q(x):
-      return numpy.exp(-x/mu)
-   mu=pars[0]
-   beta=pars[1]
-   w=pars[2]
-   t=range(len(n))
-   I=probablity_of_clinical_incidence(beta,w,mu,t,Q)
+      return np.exp(-x/mu)
+   t=25+np.linspace(0,len(n),len(n)+1)
+   I=probability_of_clinical_incidence(beta,w,mu,t,Q)
+   D=probability_of_preclinical_diagnosis(beta,w,mu,t,Q)
+   qs=np.expand_dims(s,1)
+
+   alpha=np.array([beta])
+   L=loglikelihood_function(D,I,n,s,r,alpha,beta)
+   return np.sum(L)
    
+def fI(pars,t):
+   w=pars[0]
+   mu=pars[1]
+   beta=pars[2]
+   def Q(x):
+      return np.exp(-x/mu)
+   return probability_of_clinical_incidence(beta,w,mu,t,Q)
+
+def fD(pars,t):
+   w=pars[0]
+   mu=pars[1]
+   beta=pars[2]
+   def Q(x):
+      return np.exp(-x/mu)
+   return probability_of_preclinical_diagnosis(beta,w,mu,t,Q)
+
 
 def probability_of_clinical_incidence(beta, w, mu, t, Q):
 #Calculate the probability of clinical incidence.