瀏覽代碼

Compute proportion of information used for theta plots

tforest 5 月之前
父節點
當前提交
87bef76e28
共有 1 個文件被更改,包括 11 次插入11 次删除
  1. 11 11
      swp2.py

+ 11 - 11
swp2.py 查看文件

262
             # divide by N0
262
             # divide by N0
263
             y[i] = y[i]/N0
263
             y[i] = y[i]/N0
264
             x[i] = x[i]/N0
264
             x[i] = x[i]/N0
265
-        sum_theta_i = 0
266
-        print(epoch, x, y)
267
-        for i in range(2, len(y)-1):
268
-            sum_theta_i=y[i] / (i-1)
269
-        prop = []
270
-        for k in range(2, len(y)-1):
271
-            prop.append(y[k+1] / (k - 1) / sum_theta_i)
272
-        #print(epoch, prop)
273
         plt.plot(x, y, 'o', linestyle = "-", alpha=0.75, lw=2, label = str(epoch)+' BrkPt | Lik='+greatest_likelihood)
265
         plt.plot(x, y, 'o', linestyle = "-", alpha=0.75, lw=2, label = str(epoch)+' BrkPt | Lik='+greatest_likelihood)
274
         if theta_scale:
266
         if theta_scale:
275
             plt.xlabel("Coal. time")
267
             plt.xlabel("Coal. time")
294
     # number of monomorphic sites
286
     # number of monomorphic sites
295
     L = L_stored
287
     L = L_stored
296
     S0 = L-S
288
     S0 = L-S
297
-    print("SFS", SFS_stored)
298
-    print("S", S, "L", L, "S0=", S0)
289
+    # print("SFS", SFS_stored)
290
+    # print("S", S, "L", L, "S0=", S0)
299
     # compute Ln
291
     # compute Ln
300
     Ln = log_facto(S+S0) - log_facto(S0) + np.log(float(S0)/(S+S0)) * S0
292
     Ln = log_facto(S+S0) - log_facto(S0) + np.log(float(S0)/(S+S0)) * S0
301
     for xi in range(0, len(SFS_stored)):
293
     for xi in range(0, len(SFS_stored)):
302
         p_i = SFS_stored[xi] / float(S+S0)
294
         p_i = SFS_stored[xi] / float(S+S0)
303
         Ln += np.log(p_i) * SFS_stored[xi] - log_facto(SFS_stored[xi])
295
         Ln += np.log(p_i) * SFS_stored[xi] - log_facto(SFS_stored[xi])
304
     res = Ln
296
     res = Ln
305
-    print(res)
297
+    # print(res)
306
     # basic plot likelihood
298
     # basic plot likelihood
307
     plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
299
     plt.figure(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
308
     plt.rcParams['font.size'] = '18'
300
     plt.rcParams['font.size'] = '18'
362
                 N0 = y[0]
354
                 N0 = y[0]
363
         for i in range(len(y)):
355
         for i in range(len(y)):
364
             y[i] = y[i]/N0
356
             y[i] = y[i]/N0
357
+        # compute the proportion of information used at each bin of the SFS
358
+        sum_theta_i = 0
359
+        for i in range(2, len(y)-1):
360
+            sum_theta_i+=y[i] / (i-1)
361
+        prop = []
362
+        for k in range(2, len(y)-1):
363
+            prop.append(y[k] / (k - 1) / sum_theta_i)
364
+        # plot 
365
         plt.plot(x, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks')
365
         plt.plot(x, y, 'o', linestyle="dotted", alpha=0.75, lw=2, label = str(epoch)+' brks')
366
         plt.xlabel("# breaks")
366
         plt.xlabel("# breaks")
367
         plt.ylabel("theta")
367
         plt.ylabel("theta")