2e:null 2f:["$","title",null,{"children":"kingman | Kevin R. Thornton | Observable"}] 34:[] 32:{"id":"1881e422cee5c984","slug":"kingman","title":"kingman","thumbnail":"35fb5c537545d63d59f7d0e1c77fa70216563f1afc45a872f5ecce8d39dbf092","access_level":"public","create_time":"$D2021-11-18T18:24:07.248Z","update_time":"$D2021-11-28T20:54:52.941Z","view_time":null,"first_public_time":"$D2021-11-18T18:37:49.943Z","has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"10f23817881a3073","login":"molpopgen","name":"Kevin R. Thornton","avatar_url":"https://avatars-next.observableusercontent.com/266416fe400ffcd4186b26579ed97451084368283db1df9f3d1c5f8c205ea8c8","tier":"starter_2024"},"version":382,"publish_version":382,"role":null,"authored":false,"authors":[{"type":"individual","id":"eb56f99e4b557bf7","login":"molpopgen","name":"Kevin R. Thornton","avatar_url":"https://avatars-next.observableusercontent.com/266416fe400ffcd4186b26579ed97451084368283db1df9f3d1c5f8c205ea8c8","tier":null}],"members":[{"type":"individual","id":"eb56f99e4b557bf7","login":"molpopgen","name":"Kevin R. Thornton","avatar_url":"https://avatars-next.observableusercontent.com/266416fe400ffcd4186b26579ed97451084368283db1df9f3d1c5f8c205ea8c8","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":{"id":"6cf580b27e276df3","slug":"utilities","title":"Utilities","owner":{"id":"10f23817881a3073","login":"molpopgen","name":"Kevin R. Thornton"}},"type":"notebook","body":{"title":"kingman","cells":[{"id":0,"value":"# kingman\n\nUtilities for doing stuff with the Kingman coalescent.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":267,"value":"## Expected time to most recent common ancestor (MRCA)","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":16,"value":"plot_etmrca_vs_sample_size = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n x: {\n grid: true,\n label: \"Sample size (haploid lineages).\",\n },\n y: {\n grid: true,\n label: \"E(TMRCA), units of 2Ne generations\",\n },\n marks: [\n Plot.line(etmrca_vs_sample_size, {\n x: \"n\", \n y: \"etmrca\",\n }),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":13,"value":"etmrca_vs_sample_size = {\n let rv = [];\n for(let n = 2 ; n < 50 ; ++n) {\n rv.push({'n': n, 'etmrca': etmrca(n)});\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":271,"value":"# Expeted total time on tree (aka sum of branch lengths, aka TTOT)","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":30,"value":"plot_ettot_vs_sample_size = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n x: {\n grid: true,\n label: \"Sample size (haploid lineages).\",\n },\n y: {\n grid: true,\n label: \"E(TTOT), units of 2Ne generations\",\n },\n marks: [\n Plot.line(ettot_vs_sample_size, {\n x: \"n\", \n y: \"ettot\",\n }),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":27,"value":"ettot_vs_sample_size = {\n let rv = [];\n for(let n = 2 ; n < 50 ; ++n) {\n rv.push({'n': n, 'ettot': ettot(n)});\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":275,"value":"## Expected site (or allele) frequency spectum, aka SFS or AFS","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":43,"value":"plot_example_fs = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n x: {\n grid: true,\n label: \"Number of occurrences of derived character state\",\n },\n y: {\n grid: true,\n label: \"Expected number of mutations\",\n },\n marks: [\n Plot.barY(example_fs, {\n x: \"i\",\n y: \"expected_number_mutations\",\n fill: \"pink\",\n }),\n Plot.text(example_fs, {x: \"i\", y: d => d.expected_number_mutations/2., text: d => `θ/${d.i}`}),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":41,"value":"example_fs = neutral_unfolded_fs(10, 50.);","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":185,"value":"ftmrca_data = {\n let data = [];\n let n = [2, 5, 10, 20, 50, 100];\n for (const ni of n) {\n let temp = ftmrca(ni, 0, 5, 1e-3);\n for(const t of temp) {\n data.push({\"n\": ni, \"t\": t.t, \"density\": t.density});\n }\n }\n return data;\n};","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":279,"value":"## Distribution of TMRCA","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":195,"value":"plot_ftmrca_data = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n x: {\n grid: true,\n label: \"TMRCA\",\n },\n y: {\n grid: true,\n label: \"density\",\n },\n marks: [\n Plot.line(ftmrca_data, {\n x: \"t\", \n y: \"density\",\n stroke: \"n\",\n title: (o) => `n = ${o.n}`,\n strokeWidth: 3,\n }),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":229,"value":"fttot_data = {\n let data = [];\n \n let n = [2, 5, 10, 20, 50, 100];\n for (const ni of n) {\n let temp = fttot(ni, 0, 20, 1e-3);\n for(const t of temp) {\n data.push({\"n\": ni, \"t\": t.t, \"density\": t.density});\n }\n }\n\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":283,"value":"## Distribution of total time","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":232,"value":"plot_fttot_data = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n\n x: {\n grid: true,\n label: \"Total time on tree.\",\n },\n y: {\n grid: true,\n label: \"density\",\n },\n marks: [\n Plot.line(fttot_data, {\n x: \"t\", \n y: \"density\",\n stroke: \"n\",\n title: (o) => `n = ${o.n}`,\n strokeWidth: 3,\n }),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":287,"value":"## Distribution of the number of segregating sites","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":309,"value":"plot_P_s_equals_k_2d = Plot.plot({\n style: \"overflow: visible;\",\n style: \"font-size: 11.5px\",\n x: {\n grid: true,\n label: \"Number of segregating sites.\",\n },\n y: {\n grid: true,\n label: \"Probability\",\n },\n marks: [\n Plot.dot(plot_P_s_equals_k_data, {\n x: \"k\", \n y: \"density\",\n stroke: \"n\",\n title: (o) => `n = ${o.n}`,\n strokeWidth: 1.5,\n }),\n \n Plot.ruleX(plot_P_s_equals_k_data, {\n x: \"k\", \n y: \"density\",\n stroke: \"n\",\n title: (o) => `n = ${o.n}`,\n strokeWidth: 1.5,\n strokeOpacity: 0.025,\n }),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":293,"value":"viewof max_k = Inputs.range([10, 100], {step: 1, label: \"Max number of segregating sites.\", value: 50})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":295,"value":"viewof n_min = Inputs.range([2, 10], {step: 1, label: \"Minimum sample size\", value: 2})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":301,"value":"viewof n_max = Inputs.range([n_min+1, 30], {step: 1, label: \"Maximum sample size\", value: 15})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":303,"value":"viewof theta = Inputs.range([1, 10], {step: 1e-2, label: \"θ\", value: 5})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":305,"value":"plot_P_s_equals_k_data = {\n let data = [];\n for (let n = n_min ; n <= n_max ; ++n) {\n let temp = p_S_equals_k(n, theta, 0, max_k);\n for (const t of temp) {\n data.push({'n': n, 'k': t.k, 'density': t.density});\n }\n }\n return data;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":222,"value":"## Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":10,"value":"etmrca = function(n) {\n let _etmrca = 0.0;\n for(let i = 2 ; i < n + 1 ; ++i) {\n _etmrca += 2.0/(i*(i - 1));\n }\n return _etmrca;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":24,"value":"// Also the denominator of Watterson's estimate\n// of theta = 4Ne\\mu\nettot = function(n) {\n let e = 0.0;\n for (let i = 1 ; i < n; ++i) {\n e += 1./i;\n }\n return e;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":36,"value":"neutral_unfolded_fs = function(n, theta) {\n let rv = [];\n for (let i = 1; i < n; ++i) {\n rv.push({'i': i, \"expected_number_mutations\": theta/i});\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":137,"value":"i_choose_2 = function(i) { return i*(i-1)/2.0; }","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":134,"value":"ftmrca = function(n, mint, maxt, dt) {\n // n - 2 convolutions of exponential deviates.\n // See Wakeley, p. 77 for detaills.\n let rv = [];\n let t = mint;\n while (t <= maxt) {\n let d = 0.0;\n for(let i = 2 ; i <= n ; ++i) {\n let ic2 = i_choose_2(i);\n \n let temp = ic2 * Math.exp(-ic2*t);\n let d2 = 1.0; \n for(let j = 2 ; j <= n ; ++j) {\n if (i != j)\n {\n let jc2 = i_choose_2(j);\n d2 *= jc2/(jc2 - ic2);\n }\n }\n d += temp * d2;\n }\n rv.push({\"t\": t, \"density\": d});\n t += dt;\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":227,"value":"fttot = function(n, mint, maxt, dt) {\n // Eqn 3.34 of Wakeley, p. 78\n let rv = [];\n let t = mint;\n while (t <= maxt) {\n let density = ((n-1)/2)*Math.exp(-t/2.)*(1.0 - Math.exp(-t/2.))**(n-2);\n rv.push({'t': t, 'density': density});\n t += dt;\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":247,"value":"function choose (n, k) {\n // BAD! This will blow up eventually\n if (k === 0) return 1;\n return (n * choose(n-1, k-1)) / k;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":246,"value":"p_S_equals_k = function(n, theta, k_min, k_max) {\n // Tavare, 1984.\n // Wakeley, p. 94\n let rv = [];\n let k = k_min;\n while (k <= k_max) {\n let density = 0.;\n for(let i=2; i <= n ; ++i) {\n let temp = ((-1)**i)*choose(n-1, i-1);\n temp *= (i-1)/(theta + i - 1);\n temp *= (theta/(theta + i - 1))**k;\n density += temp;\n }\n rv.push({\"k\":k, \"density\": density});\n k += 1;\n }\n return rv;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W34","stdlib":"1","titling":"auto"}} 33:"$32:body"