Correcting sum terms in probablity of clinical incidence and preclinical diagnosis

three_parameter_screening.py

@@ -175,42 +175,54 @@ def likelihood_function(D, I, n, s, r, alpha, beta):
     """
     L = 1.0
     for i in range(len(D)):
-        term1 = D[i] ** s[i]
+        s=numpy.sum(s[i])
+        term1 = D[i] ** s
         term2 = I[i] ** r[i]
-        term3 = (1 - D[i] - I[i]) ** (n[i] - s[i] - r[i])
-        product_term = np.prod([(alpha[j] / beta) ** s[i][j] for j in range(3)])
+        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))])
         L_i = term1 * term2 * term3 * product_term
         L *= L_i
     return L
 
-def probability_of_clinical_incidence(beta, w, mu, t, Q):
-    """
-    Calculate the probability of clinical incidence.
-
-    Parameters
-    ----------
-    w : float
-        Weighting factor.
-    mu : float
-        Mean of the distribution.
-    t : list
-        List of time intervals.
-    Q : function
-        Function Q(t) representing the probability distribution.
-    beta : float
-        The probability of detection.
+def likelihood_function(pars,n,s,r):
+   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)
+   
 
-    Returns
-    -------
-    list
-        A list of probabilities I_i for each time interval.
-    """
+def probability_of_clinical_incidence(beta, w, mu, t, Q):
+#Calculate the probability of clinical incidence.
+
+#    Parameters
+#    ----------
+#    w : float
+#        Weighting factor.
+#    mu : float
+#        Mean of the distribution.
+#    t : list
+#        List of screening times [0,..,tN] of length N+1 (where N is number intervals)
+#    Q : function
+#        Function Q(t) representing the probability distribution.
+#    beta : float
+#       The probability of detection.
+#
+#    Returns
+#    -------
+#    list
+#        A list of probabilities I_i for each time interval (of length N, 0 corresponds to interval [t0,t1))
+#    """
     I = []
-    for i in range(1, len(t) + 1):
-        sum_term = sum((1 - beta) ** (i - j - 1) * (Q(t[i - 1] - t[j]) - Q(t[i] - t[j])) for j in range(i))
-        I_i = w * mu * ((t[i] - t[i - 1]) / mu - beta * sum_term)
-        I.append(I_i)
-    return I
+    for i in range(1,len(t)):
+       sum_terms=[(1 - beta) ** (i - j - 1) * (Q(t[i - 1] - t[j]) - Q(t[i] - t[j])) for j in range(i)]
+       sum_term = np.sum(sum_terms)
+       dt=(t[i]-t[i-1])
+       I_i = w * (dt - beta * mu * sum_term)
+       I.append(I_i)
+    return np.array(I)
 
 def probability_of_preclinical_diagnosis(beta, w, mu, t, Q):
     """
@@ -235,14 +247,16 @@ def probability_of_preclinical_diagnosis(beta, w, mu, t, Q):
         A list of probabilities D_i for each time interval.
     """
     D = []
-    for i in range(1, len(t) + 1):
+    for i in range(1, len(t) ):
         if i == 1:
-            D_i = beta * w * mu
+           D_i = beta * w * mu
         else:
-            sum_term = sum((1 - beta) ** (i - j - 1) * Q(t[i - 1] - t[j - 1]) for j in range(1, i))
-            D_i = beta * w * mu * (1 - beta * sum_term)
+           sum_terms=[(1 - beta) ** (i - j - 1) * Q(t[i - 1] - t[j - 1]) for j in range(1, i)]
+           sum_term = np.sum(sum_terms) 
+           D_i = beta * w * mu * (1 - beta * sum_term)
         D.append(D_i)
-    return D
+    return np.array(D)
+
 
 # endregion
 # region Optimization
@@ -304,4 +318,4 @@ def analytic_simulation(sensitivity, specificity, risk, screening_interval, onse
   
     return np.array([cancer_quantity, screening_detected_q, interval_cancer_q, recall_q, recall_total])
 
-# endregion
+# endregion