The expression is the Binet formula for the Fibonacci numbers. Compute the two roots and round the result to the nearest integer:
from math import sqrt
def binet(n):
phi = (1 + sqrt(5)) / 2
psi = (1 - sqrt(5)) / 2
return round((phi**n - psi**n) / sqrt(5))list50 = [binet(n) for n in range(1, 51)]Use the first 10, 20, 30 and 40 entries of the already computed list in a two-by-two array of axes. Each subplot needs an index sequence starting at one and having the same length as its data:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 2, figsize=(8, 6), layout="constrained")
for ax, stop in zip(axes.flat, [10, 20, 30, 40]):
ax.plot(range(1, stop + 1), list50[:stop])
ax.set(title=f"n = {stop}", xlabel="n", ylabel="f(n)")
plt.show()list50, even though the requested four panels stop at 40.
TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
MultiPathogenData.csv dataset. The file was not supplied; no numerical output is fabricated.The first instruction asks for maximization, whereas the later parts refer to minimization. For the literal maximization task, minimize the negative objective:
import numpy as np
from scipy import optimize as opt
def f(x):
return 1 + np.exp(-x*x/5) * np.sin(5*x)
maximum = opt.minimize(lambda z: -f(z[0]), x0=[-1.0])
print(maximum.x[0], f(maximum.x[0]))opt.minimize(lambda z: f(z[0]), x0=[-1.0]), which instead returns a local minimum near , . The optimizer's stopping criterion establishes local numerical convergence, not a certificate of a global optimum.Plot the entire stated interval:
import matplotlib.pyplot as plt
x_grid = np.linspace(-4, 4, 1500)
plt.plot(x_grid, f(x_grid))
plt.xlabel("x")
plt.ylabel("f(x)")
plt.show()
Local extrema, the global minimum and the constrained minimum of the oscillatory objective
. Details. 
Use reproducible random starts and retain the best successful result:
rng = np.random.default_rng(2024)
starts = rng.uniform(-4, 4, size=100)
results = [opt.minimize(lambda z: f(z[0]), x0=[x0]) for x0 in starts]
successful = [result for result in results if result.success]
best = min(successful, key=lambda result: result.fun)
print(best.x[0], best.fun)Here an analytic check is possible. For , the envelope gives , above the candidate. On , . In the remaining interval , the derivative condition is equivalent, with , to . The left side decreases strictly from to while the right side increases, giving one minimum. This certifies the selected well as global rather than relying only on the restart experiment.
def transformed(y):
return f(1.5 + np.exp(y[0]))
result = opt.minimize(transformed, x0=[0.0])
x_interior = 1.5 + np.exp(result.x[0])
x_best = min([1.5, x_interior], key=f)
print(x_best, f(x_best))Constraint reparameterization enforces feasibility but does not remove nonconvexity. In this example, up to the next negative-sine lobe, the derivative has one minimum in that lobe, and for the envelope gives . Thus later oscillations cannot improve on the displayed constrained minimum.
Write for the true DNA methylation status and for its observation. The observation likelihood function is a binary symmetric channel, soFor the loop suggested in the practical, set and add when , or otherwise. The following independent implementation combines these terms without multiplying many small numbers:
import numpy as np
def methylation_score(pattern, observed):
"""Conditional-product score used by the practical's natural loop filling."""
pattern=np.asarray(pattern,dtype=int);observed=np.asarray(observed,dtype=int)
if pattern[0]!=1 or pattern[-1]!=1:return -np.inf
neighbour_sum=pattern[:-2]+pattern[2:]
p=np.array([.1,.5,.9])[neighbour_sum]
prior_terms=np.where(pattern[1:-1]==1,p,1-p)
observation_terms=np.where(pattern==observed,.8,.2)
return float(np.log(prior_terms).sum()+np.log(observation_terms).sum())
log_likelihood = methylation_score
print(f"{log_likelihood(obs, obs):.2f}")The conditional-product log score at the observed pattern is . The endpoint observations are included once each, while the endpoint statuses are fixed rather than assigned an additional prior probability.
There is a modeling qualification: multiplying the specified neighbour conditional probabilities gives a pseudolikelihood, not generally the joint probability distribution having those conditional probabilities. Thus the function above defines the practical's conditional-product target; the pure observation log-likelihood at would instead be .
If the neighbour probabilities are intended literally as full conditionals of a Markov random field, a coherent alternative has edge weights and for , with both endpoints fixed. Given two equal neighbours, the matching centre has conditional odds ; mixed neighbours give equal weights. Consequently its unnormalized posterior probability has log scoreFor this alternative, before the constant associated with conditioning the endpoints. It can be implemented as follows, but its numerical results must not be mixed with the conditional-product results in the remaining parts:
def joint_methylation_score(pattern,observed):
"""Coherent nearest-neighbour model with the specified full conditionals."""
pattern=np.asarray(pattern,dtype=int);observed=np.asarray(observed,dtype=int)
if pattern[0]!=1 or pattern[-1]!=1:return -np.inf
edges=np.where(pattern[1:]==pattern[:-1],.75,.25)
emissions=np.where(pattern==observed,.8,.2)
return float(np.log(edges).sum()+np.log(emissions).sum())A single randomly chosen internal-site flip is a symmetric proposal: the reverse flip has the same probability . Therefore the Metropolis–Hastings algorithm accepts withComparing logarithms avoids overflow. On rejection the old state must be retained, including a repeated copy in the stored Markov chain; retaining only accepted states would sample a different distribution. This implementation also stores the actual initial state at row zero:
def sample_methylation(observed, nsteps=5000, seed=2024, score=methylation_score):
rng=np.random.default_rng(seed)
state=rng.integers(0,2,size=len(observed));state[[0,-1]]=1
states=np.empty((nsteps,len(state)),dtype=int);values=np.empty(nsteps)
states[0]=state;values[0]=score(state,observed)
accepted=0
for t in range(1,nsteps):
candidate=state.copy();site=rng.integers(1,len(state)-1)
candidate[site]=1-candidate[site]
trial=score(candidate,observed)
if np.log(rng.random())<min(0.,trial-values[t-1]):
state=candidate;values[t]=trial;accepted+=1
else:values[t]=values[t-1]
states[t]=state
return values,states,accepted/(nsteps-1)
lls, meths, acceptance = sample_methylation(obs, nsteps=5000, seed=2024)
import matplotlib.pyplot as plt
plt.plot(lls)
plt.xlabel("Stored state index")
plt.ylabel("Conditional-product log score")
plt.show()For this random-number convention, the acceptance fraction is approximately . The trace climbs substantially during its first several thousand updates, so a small fixed burn-in such as 100 iterations would discard too little of the initial transient. Five thousand single-site proposals are only about fifty updates per site on average.
Conditional-product sampling trace and retained-sample methylation profile, with a provisional burn-in of 3000 stored states
. Details. 
For a binary state, its expected value is its probability of being one. Hence the Monte Carlo method estimate of the posterior methylation profile under the chosen target isThe averaging axis must be zero: average over sampled patterns while keeping the site index. A provisional choice for this particular trace is , retaining 2000 states:
def mean_methylation(meths, burnin):
if not 0 <= burnin < len(meths):
raise ValueError("burn-in must leave at least one state")
return np.asarray(meths)[burnin:].mean(axis=0)
burnin = 3000
profile_estimate = mean_methylation(meths, burnin)
print(lls[burnin:].mean())
plt.plot(profile_estimate, label="retained-sample mean")
plt.scatter(np.arange(len(obs)), .2*obs, s=10, label="0.2 × observation")
plt.xlabel("Site index")
plt.ylabel("Methylation probability")
plt.legend()
plt.show()The retained mean log score is , and the profile is shown below the trace in part (b). The fixed endpoints have estimate exactly one; adjacent sites are often smoothed into methylated or unmethylated stretches despite individual conflicting observations. The rescaled observations are only a visual comparison and are not another estimate of the posterior probability.
This is a finite-chain estimate for the conditional-product target defined in part (a). A burn-in chosen by inspection does not establish convergence, and the long excursions in this short Markov chain leave substantial uncertainty in individual-site estimates; several longer runs would be needed for precise inference. To sample the coherent joint alternative instead, pass
score=joint_methylation_score and inspect that new trace before choosing its own burn-in.The retained mean conditional-product log score, , is higher than the observed-pattern score, , because sampled patterns can improve agreement between neighbouring CpG sites while tolerating some observational errors. The observation-matching pattern maximizes the observation likelihood function alone, whereas this score also rewards the assumed spatial structure, and its posterior mean is not its maximum.
Generate independent Bernoulli random variables above the diagonal and set everything else to zero. Taking the strict upper triangle of a random square array gives each allowed entry probability of being one:
import numpy as np
rng = np.random.default_rng(2024)
def rand_triu(n):
return np.triu(rng.integers(0, 2, size=(n, n)), k=1)
example = rand_triu(10)
print(example)One reproducible strictly upper triangular matrix is:
0 1 0 0 0 0 1 1 1 1
0 0 1 0 0 0 1 0 0 0
0 0 0 1 1 0 0 1 1 0
0 0 0 0 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0Every matrix returned in this way is a nilpotent matrix; randomness affects its nilpotency index, not the fact of the nilpotent property.
Starting with the identity matrix, multiply by once per iteration and test the resulting power. Using integer objects keeps the calculation exact and avoids fixed-width integer overflow accidentally turning a nonzero power into zero:
def kzero(matrix,kmax=100):
a=np.asarray(matrix,dtype=object)
if a.ndim!=2 or a.shape[0]!=a.shape[1]:raise ValueError('square matrix required')
power=np.eye(len(a),dtype=object)
for k in range(1,kmax+1):
power=power@a
if np.equal(power,0).all():return k
return None
print([kzero(rand_triu(10)) for _ in range(5)])Continuing the random-number stream from the displayed example gives . These are the least positive zero powers, rather than the number of nonzero powers; in particular, a zero matrix has nilpotency index one.
For a strictly upper triangular matrix, a nonzero entry of can only lie at least positions above the diagonal. Thus , and every one of these examples must return an integer at most ten. The code interprets “up to” as including
kmax; if the final strict inequality in the question is intended instead, replace range(1,kmax+1) by range(1,kmax).TODO: The numerical/data-specific answer awaits the original
clustering_dataset.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
clustering_dataset.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
clustering_dataset.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
clustering_dataset.csv dataset. The file was not supplied; no numerical output is fabricated.Use the supplied BLOSUM62 entry for every diagonal move in the Smith-Waterman algorithm, while retaining the stipulated linear gap penalty of . If and are the current amino acids, the dynamic programming recurrence isThe zero option permits a new local start. The following original implementation also records the traceback coordinates needed in part (c); its tie ordering is left, up, diagonal, and an equal overall maximum replaces the previous best cell, as in the supplied convention:
def local_alignment(query,reference,alphabet,scores,gap=-1):
lookup={aa:i for i,aa in enumerate(alphabet)}
table=np.asarray(scores,dtype=int)
height,width=len(query)+1,len(reference)+1
H=np.zeros((height,width),dtype=int)
back=np.zeros((height,width),dtype=np.int8)
best=(0,0,0)
for i,a in enumerate(query,1):
for j,b in enumerate(reference,1):
# Order implements the printed LEFT, UP, DIAGONAL tie priority.
candidates=(0,H[i,j-1]+gap,H[i-1,j]+gap,
H[i-1,j-1]+int(table[lookup[a],lookup[b]]))
value=max(candidates);H[i,j]=value
back[i,j]=0 if value==0 else candidates.index(value)
if value>=best[0]:best=(value,i,j)
score,i,j=best;query_end,ref_end=i,j
aq=[];ar=[]
while H[i,j]>0:
step=back[i,j]
if step==1:aq.append('-');ar.append(reference[j-1]);j-=1
elif step==2:aq.append(query[i-1]);ar.append('-');i-=1
elif step==3:aq.append(query[i-1]);ar.append(reference[j-1]);i-=1;j-=1
qa=''.join(reversed(aq));ra=''.join(reversed(ar));length=len(qa)
identities=sum(a==b and a!='-' for a,b in zip(qa,ra))
similarities=sum(a!='-' and b!='-' and table[lookup[a],lookup[b]]>0 for a,b in zip(qa,ra))
gaps=sum(a=='-' or b=='-' for a,b in zip(qa,ra))
return dict(query=qa,reference=ra,score=int(score),query_start=i+1,query_end=query_end,
reference_start=j+1,reference_end=ref_end,query_length=len(query),reference_length=len(reference),
alignment_length=length,identities=int(identities),similarities=int(similarities),gaps=int(gaps),
identity_percent=100*identities/length if length else 0.,
similarity_percent=100*similarities/length if length else 0.,gap_percent=100*gaps/length if length else 0.)
queries = {"baboon": baboon, "squirrel_monkey": squirrel_monkey,
"dingo": dingo, "goldfish": goldfish}
alignments = {
name: local_alignment(seq, reference_epo_receptor_human, amino_acids, bmat)
for name, seq in queries.items()
}
for name, result in alignments.items():
print(name, "score", result["score"])
print(result["query"])
print(result["reference"])The resulting local sequence alignments, with the query above the reference, are:
baboon; score 314
SLLTPSDLDPLILTLSLILVLIVVLLTVLALLSHRRALKQKIWPGIPSPESEFEGLFTTHKG-FQL
SLLTPSDLDPLILTLSLILVVILVLLTVLALLSHRRALKQKIWPGIPSPESEFEGLFTTHKGNFQLsquirrel monkey; score 278
VCFWEEAASAGMGPDNYSFSYQL-DNEPWKR-CRLHQTPTARGVVRFWCSLPRADTSS
VCFWEEAASAGVGPGNYSFSYQLED-EPWK-LCRLHQAPTARGAVRFWCSLPTADTSSdingo; score 301
MNHLW--THLWPG-VGSLCLLLAGAAWASLPKP-L-DPKFESKAALLAARAPEELLCFTERLEDLVCFWEE
MDHL-GAS-LWP-QVGSLCLLLAGAAWAP-P-PNLPDPKFESKAALLAARGPEELLCFTERLEDLVCFWEEgoldfish; score 190
L--LVLSIGLILC--L-LSLTVI-L-SHQKFLLK-KLWPDIPTPEHKFPGLFTVYKGDFKE-WMSH-NNGSM
LDPLILTLSLILVVILVL-LTVLALLSHRR-ALKQKIWPGIPSPESEFEGLFTTHKGNF-QLWL-YQNDGCLThe several short gaps are a consequence of the specified inexpensive linear gap penalty; substituting a different default or an affine gap penalty would solve a different scoring problem.
Use the number of alignment columns, including gap columns, as the denominator. A column contributes to sequence identity when its two nongap letters are equal, to sequence similarity when its two nongap letters have a strictly positive BLOSUM62 score, and to the gap fraction when either letter is a gap. A zero substitution score is not counted as similar.
The implementation in part (a) computes each numerator directly during its final pass through the aligned strings:
for name, a in alignments.items():
print(name, a["alignment_length"], a["identities"],
a["similarities"], a["gaps"])
print(a["identity_percent"], a["similarity_percent"], a["gap_percent"])Identical columns also count as similar for these data, so these percentages are overlapping summaries rather than three categories summing to 100%.
The maximum-scoring cell records the inclusive end position in each sequence. During traceback, moving left consumes a reference residue, moving up consumes a query residue, and moving diagonally consumes both; when traceback stops, the remaining prefix lengths are one less than the corresponding one-based starts. Consequently
i+1 and j+1 are the correct starts returned in part (a), without counting any inserted gaps as sequence residues.The following report joins these coordinates to the alignments and sequence identity/sequence similarity values already calculated:
for name, a in alignments.items():
print(f"{name}: query length={a['query_length']}, "
f"reference length={a['reference_length']}, "
f"alignment columns={a['alignment_length']}, score={a['score']}")
print(f"Query {a['query_start']} {a['query']} {a['query_end']}")
print(f"Ref {a['reference_start']} {a['reference']} {a['reference_end']}")
print(f"Identity {a['identity_percent']:.2f}%; "
f"similarity {a['similarity_percent']:.2f}%; "
f"gaps {a['gap_percent']:.2f}%")All reference sequences have length 508. The following are one-based inclusive coordinates; the aligned strings and all percentages are given in parts (a) and (b).
For example, the goldfish local sequence alignment excludes its first query residue and aligns residues 2–63. As a check, removing gaps from each aligned string gives exactly the sequence slice determined by its reported interval.
Convert the supplied multiple sequence alignment to an ordinary array so that reducing over rows returns one value per column. With nine sequences, fewer than 30% gaps means at most two gaps in a column:
from collections import Counter
A = np.asarray(msa_matrix)
gap_fraction = np.mean(A == ".", axis=0)
match_mask = gap_fraction < .30
match_columns = np.flatnonzero(match_mask)
print(match_columns + 1) # human-readable, one-based alignment columnsThere are 72 match states: alignment columns 1–21 and 29–79, using one-based numbering. The seven remaining columns, 22–28, form one contiguous insertion region. The strict inequality is applied to the actual fraction, rather than to rounded percentages.
For each match state, count the nongap amino acids and divide by the number of nongap observations in that column. A gap is not one of the twenty emitted amino acids; including it in the denominator would make the emission probabilities sum to less than one:
def emission_counts(column):
counts = Counter(str(x) for x in column if x != ".")
n = sum(counts.values())
if n == 0:
raise ValueError("an all-gap column has no empirical emission distribution")
return counts, n
match_emissions = {}
for state, j in enumerate(match_columns, 1):
counts, n = emission_counts(A[:, j])
match_emissions[state] = {aa: count/n for aa, count in counts.items()}
print(state, j+1, match_emissions[state])The table gives an exact compact representation of every emission probability: divide each displayed count by the denominator in its row; all unlisted amino acids have probability zero. The state number follows the order of match states, while the column refers to the original multiple sequence alignment.
| State | Alignment column | Nongap counts | |
|---|---|---|---|
| 1 | 1 | 9 | I:1, K:1, R:1, S:2, V:3, Y:1 |
| 2 | 2 | 9 | D:1, E:2, G:1, I:1, L:1, Q:2, S:1 |
| 3 | 3 | 9 | E:2, I:1, K:1, Q:1, R:3, T:1 |
| 4 | 4 | 9 | L:8, M:1 |
| 5 | 5 | 9 | E:3, G:1, L:2, N:1, Q:1, S:1 |
| 6 | 6 | 9 | A:1, D:2, G:1, K:1, N:1, S:2, T:1 |
| 7 | 7 | 9 | A:1, E:3, I:2, L:1, M:1, V:1 |
| 8 | 8 | 9 | I:2, L:7 |
| 9 | 9 | 9 | D:3, E:1, K:1, L:2, N:1, Q:1 |
| 10 | 10 | 9 | D:1, E:3, G:1, N:1, P:1, S:2 |
| 11 | 11 | 9 | A:5, D:1, G:1, P:1, T:1 |
| 12 | 12 | 9 | I:1, L:7, Y:1 |
| 13 | 13 | 9 | E:1, F:1, G:1, H:2, M:1, Q:2, T:1 |
| 14 | 14 | 9 | F:9 |
| 15 | 15 | 9 | F:1, L:8 |
| 16 | 16 | 9 | E:1, I:1, L:1, N:3, Q:1, R:1, Y:1 |
| 17 | 17 | 9 | D:6, E:2, S:1 |
| 18 | 18 | 9 | C:7, F:1, S:1 |
| 19 | 19 | 9 | K:1, N:6, R:2 |
| 20 | 20 | 9 | E:1, F:1, I:5, L:1, V:1 |
| 21 | 21 | 9 | K:3, P:1, Q:1, S:4 |
| 22 | 29 | 7 | K:4, P:1, R:1, S:1 |
| 23 | 30 | 7 | E:2, F:1, G:1, I:1, K:1, V:1 |
| 24 | 31 | 7 | D:3, E:1, G:1, K:1, S:1 |
| 25 | 32 | 8 | D:1, E:1, G:1, K:1, N:1, R:3 |
| 26 | 33 | 8 | D:7, N:1 |
| 27 | 34 | 8 | A:1, P:6, S:1 |
| 28 | 35 | 9 | C:1, H:1, I:1, K:1, M:1, R:1, T:3 |
| 29 | 36 | 9 | E:1, F:1, G:1, L:1, Q:3, W:1, Y:1 |
| 30 | 37 | 9 | I:5, L:1, T:2, Y:1 |
| 31 | 38 | 9 | A:1, G:1, P:3, S:4 |
| 32 | 39 | 9 | I:1, K:6, L:1, R:1 |
| 33 | 40 | 9 | F:1, N:1, Q:4, T:2, V:1 |
| 34 | 41 | 9 | A:1, F:1, I:1, V:6 |
| 35 | 42 | 9 | I:1, L:7, V:1 |
| 36 | 43 | 9 | A:1, C:1, E:1, N:1, R:1, S:3, V:1 |
| 37 | 44 | 9 | D:3, E:5, Q:1 |
| 38 | 45 | 9 | A:1, C:7, I:1 |
| 39 | 46 | 9 | K:2, L:3, R:4 |
| 40 | 47 | 9 | A:3, E:1, L:1, N:3, T:1 |
| 41 | 48 | 9 | I:4, L:1, T:1, V:3 |
| 42 | 49 | 9 | L:9 |
| 43 | 50 | 9 | D:1, E:1, H:1, Q:1, S:1, T:1, V:1, Y:2 |
| 44 | 51 | 9 | E:1, K:1, Q:2, R:1, V:4 |
| 45 | 52 | 9 | L:9 |
| 46 | 53 | 9 | G:9 |
| 47 | 54 | 9 | P:8, T:1 |
| 48 | 55 | 9 | W:9 |
| 49 | 56 | 9 | A:2, C:7 |
| 50 | 57 | 9 | A:6, C:1, L:1, T:1 |
| 51 | 58 | 9 | A:1, D:3, F:1, S:2, V:1, Y:1 |
| 52 | 59 | 9 | K:4, R:2, S:3 |
| 53 | 60 | 9 | A:2, F:1, I:2, L:2, V:2 |
| 54 | 61 | 9 | A:7, G:1, Y:1 |
| 55 | 62 | 9 | E:3, G:3, Q:3 |
| 56 | 63 | 9 | I:3, M:3, T:1, Y:2 |
| 57 | 64 | 9 | F:3, L:2, M:3, N:1 |
| 58 | 65 | 9 | H:1, I:3, L:1, V:4 |
| 59 | 66 | 9 | K:1, M:1, P:1, R:4, S:1, T:1 |
| 60 | 67 | 9 | D:1, E:3, Q:4, Y:1 |
| 61 | 68 | 9 | I:1, L:7, M:1 |
| 62 | 69 | 9 | D:1, E:4, Q:4 |
| 63 | 70 | 9 | K:7, R:2 |
| 64 | 71 | 9 | A:1, I:2, M:1, N:1, R:1, Y:3 |
| 65 | 72 | 9 | A:1, C:1, D:1, E:2, I:2, V:2 |
| 66 | 73 | 9 | K:9 |
| 67 | 74 | 9 | H:4, L:1, Q:2, S:1, V:1 |
| 68 | 75 | 9 | E:8, Q:1 |
| 69 | 76 | 9 | Q:3, R:1, S:2, T:3 |
| 70 | 77 | 9 | C:2, E:3, G:1, R:1, S:1, V:1 |
| 71 | 78 | 9 | A:1, D:2, E:3, I:1, K:1, Y:1 |
| 72 | 79 | 9 | I:2, L:3, M:1, N:1, Q:1, V:1 |
These are unregularized maximum-likelihood estimates of the emission probabilities. No additional artificial observations are added. Although only emissions are requested here, a full profile hidden Markov model would also require transitions and deletion mechanisms to represent the gaps that remain in some match columns.
For each alignment column, ignore gaps, form the empirical distribution of its amino acids, and calculate its Shannon entropy in bits. The requested alignment column information content isThus it is the relative entropy from the column distribution to the uniform distribution on twenty amino acids, not the ordinary Shannon entropy itself. A column containing one residue type has the maximum value bits.
information = []
for j in range(A.shape[1]):
counts, n = emission_counts(A[:, j])
p = np.array(list(counts.values()), dtype=float) / n
H = -np.sum(p * np.log2(p))
information.append(float(np.log2(20) - H))
print(np.round(information, 4))The values for all 79 alignment columns are as follows, in bits. Each line lists its one-based column interval followed by the corresponding values in order:
01–10: 1.9025, 1.5964, 1.9025, 3.8187, 1.9025, 1.5964, 1.9025, 3.5577, 1.9025, 1.9025
11–20: 2.4420, 3.3355, 1.5964, 4.3219, 3.8187, 1.6803, 3.0975, 3.3355, 3.0975, 2.4420
21–30: 2.5692, 2.8219, 3.5106, 3.3219, 4.3219, 4.3219, 2.8000, 2.0703, 2.6574, 1.8003
31–40: 2.1938, 1.9163, 3.7784, 3.2606, 1.6803, 1.6803, 2.6642, 2.5692, 2.8753, 2.2631
41–50: 2.8753, 3.3355, 1.6803, 2.9703, 3.3355, 2.7914, 2.2086, 2.5692, 4.3219, 1.3742
51–60: 2.2631, 4.3219, 4.3219, 3.8187, 4.3219, 3.5577, 2.8753, 1.9025, 2.7914, 2.0409
61–70: 3.3355, 2.7370, 2.4309, 2.4309, 2.5692, 2.0409, 2.5692, 3.3355, 2.9298, 3.5577
71–79: 1.9025, 1.8187, 4.3219, 2.2631, 3.8187, 2.4309, 1.9025, 1.9025, 1.9025A high score indicates conservation among the observed nongap letters, but does not itself imply low gap fraction. For example, an insertion column represented by only one observed nongap residue has maximal information by this calculation; its tiny sample size still makes that estimate uncertain.
Columns 22–28 form one run between alignment columns 21 and 29, corresponding to match states and . A self-repeating insertion state can emit a variable number of residues, so these seven alignment columns require one populated insertion state, not seven successive mandatory states.
Pool all nongap residues in this run to estimate its emission probabilities:
insert_columns = np.flatnonzero(~match_mask)
runs = np.split(insert_columns, np.flatnonzero(np.diff(insert_columns) > 1) + 1)
insert_emissions = []
for run in runs:
if len(run):
counts, n = emission_counts(A[:, run].ravel())
insert_emissions.append({aa: count/n for aa, count in counts.items()})
print(run+1, n, insert_emissions[-1])There are 24 emitted residues, givingAll other amino acids have probability zero in this empirical emission probabilities. A conventional full architecture may provide an optional insertion state at every boundary, including the two ends, giving possible insertion positions; only one has an observed insertion here and is required by this alignment's populated insertion region.
Take to be the total instantaneous mutation rate away from a base, so each specified alternative has rate . Under the Jukes-Cantor model, the uniform component is stationary and the three contrast components decay as , givingAn implementation is:
def jukes(a, b, t, mu):
decay = np.exp(-4*mu*t/3)
return .25 + .75*decay if a == b else .25*(1-decay)
print(f"{jukes(2,2,1,.5):.3f} {jukes(2,3,1,.5):.3f}")The two probabilities are and . The probabilities in any row sum to one, at the transition matrix is the identity, and at long times every entry tends to .
The Jukes-Cantor model is also parameterized by a rate for each specified replacement, in which case the exponential is and . The total-rate convention is used consistently below; it yields an optimum within the stated search interval, whereas the corresponding per-replacement rate is about .
Initialize all eight leaves, not seven: every observed base contributes a one-hot vector in the order A, C, G, T. The reference to leaf rows in the prose is inconsistent with its eight-leaf tree.
def initialise_prob(leafList):
alphabet = {base: j for j, base in enumerate("ACGT")}
result = np.zeros((2*len(leafList)-1, 4), dtype=float)
for i, base in enumerate(leafList):
result[i, alphabet[base]] = 1.
return result
leaves = ["G", "A", "T", "G", "G", "G", "C", "G"]
prob = initialise_prob(leaves)
print(prob[:len(leaves)])The initialized leaf rows are:
0 0 1 0
1 0 0 0
0 0 0 1
0 0 1 0
0 0 1 0
0 0 1 0
0 1 0 0
0 0 1 0The remaining seven internal-node rows initially contain zeros. Their eventual entries are conditional likelihoods of descendant observations given the ancestral base, rather than normalized posterior probabilities of that ancestral base.
For node with children , the phylogenetic pruning algorithm combines the two conditionally independent descendant subtrees:The likelihood function is under a uniform ancestral prior. In the supplied numbering, children have smaller indices than their parents, so increasing node order is a valid bottom-up traversal:
def felsenstein(prob, mu, L, left, right, left_time, right_time):
# All four entries are overwritten at every internal node.
for node in range(L, 2*L-1):
Tl = np.array([[jukes(a,b,left_time[node],mu)
for b in range(4)] for a in range(4)])
Tr = np.array([[jukes(a,b,right_time[node],mu)
for b in range(4)] for a in range(4)])
prob[node] = (Tl @ prob[left[node]]) * (Tr @ prob[right[node]])
return float(np.log(prob[-1].mean()))
print(felsenstein(prob, .5, len(leaves), left, right, left_time, right_time))At , the log-likelihood is . Internal rows must remain unnormalized in this implementation, because their magnitudes contribute to the tree likelihood function; normalization without retaining the removed factors would lose that information. The final root normalization required for a posterior probability is a separate step in part (e).
Evaluate the log-likelihood function on the specified interval, retaining the rate with the largest value. A grid spacing of finds ; refining locally avoids relying on a coarse-grid location when reporting three significant figures:
def evaluate(mu):
p = initialise_prob(leaves)
return felsenstein(p, mu, len(leaves), left, right, left_time, right_time)
coarse = np.linspace(.1, 1., 901)
centre = coarse[np.argmax([evaluate(mu) for mu in coarse])]
fine = np.linspace(max(.1, centre-.001), min(1., centre+.001), 2001)
mu_max = float(fine[np.argmax([evaluate(mu) for mu in fine])])
print(f"mu = {mu_max:.3g}; log-likelihood = {evaluate(mu_max):.8f}")The maximum-likelihood estimate is to three significant figures, with refined value about and log-likelihood . Under the alternative per-replacement convention this same fitted transition matrix has ; these are equivalent parameterizations, not competing fitted trees.
Recompute the root conditional likelihoods at the fitted rate, then multiply by the uniform root prior and normalize. The factors cancel between numerator and denominator:
prob = initialise_prob(leaves)
felsenstein(prob, mu_max, len(leaves), left, right, left_time, right_time)
root_posterior = prob[-1] / prob[-1].sum()
print(dict(zip("ACGT", root_posterior)))At the refined optimum, the posterior distribution isThus G is the most probable ancestral base, but it is not inferred with certainty. This conditions on the fitted rate and does not additionally average over uncertainty in that rate.
TODO: The numerical/data-specific answer awaits the original
SerologyData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
SerologyData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
SerologyData.csv dataset. The file was not supplied; no numerical output is fabricated.TODO: The numerical/data-specific answer awaits the original
SerologyData.csv dataset. The file was not supplied; no numerical output is fabricated.The Gillespie algorithm draws a waiting time with total event rate , then chooses infection with probability and recovery otherwise. Each event conserves , and is absorbing in this closed stochastic SIR model.
Choose a horizon of 100 days, which is longer than all extinction times in the runs reported here. This independent implementation also avoids executing a sampled event beyond the horizon:
import numpy as np
import matplotlib.pyplot as plt
EPIDEMIC_RNG = np.random.default_rng(2345)
def sir_sim(t_max,X0,beta,gamma,rng=None):
if rng is None:rng=EPIDEMIC_RNG
state=np.array(X0,dtype=int);population=state.sum();time=0.
times=[time];states=[state.copy()]
while time<t_max and state[1]>0:
infection=beta*state[0]*state[1]/population;recovery=gamma*state[1]
rate=infection+recovery
if rate==0:
times.append(float(t_max));states.append(state.copy());break
next_time=time+rng.exponential(1/rate)
if next_time>t_max:
times.append(float(t_max));states.append(state.copy());break
if rng.random()<infection/rate:state+=np.array([-1,1,0])
else:state+=np.array([0,-1,1])
time=next_time;times.append(time);states.append(state.copy())
return np.array(times),np.array(states)
rng = np.random.default_rng(2345)
trajectories = [sir_sim(100, [99,1,0], 1, .25, rng) for _ in range(10)]
for times, states in trajectories:
plt.step(times, states[:,1], where="post", alpha=.8)
plt.xlabel("Time (days)")
plt.ylabel("Infected individuals")
plt.show()A step plot represents the piecewise constant infected population between events. The trajectories show both substantial outbreaks and early epidemic extinction despite the large basic reproduction number .
Ten original stochastic SIR trajectories, with frequency-dependent transmission and one initial infection
. Details. 
One simulation contributes its largest observed infected population and its last recorded time. Repeating this operation produces one row per realization:
def sir_sim_stats(n_sim,t_max,X0,beta,gamma,rng=None):
result=np.empty((n_sim,2))
for run in range(n_sim):
t,s=sir_sim(t_max,X0,beta,gamma,rng)
result[run]=s[:,1].max(),t[-1]
return resultThe maximum must be taken over the infected column,
states[:,1], rather than across all compartments. With the implementation in part (a), the final time is the epidemic extinction time if extinction occurs first, and otherwise the specified horizon; an outbreak still active at that horizon has a censored extinction time.Continue the random-number stream used for the ten trajectories, then generate 1000 further independent realizations with the same model parameters. Plot each statistic in its own histogram:
stats = sir_sim_stats(1000, 100, [99,1,0], 1, .25, rng)
fig, axes = plt.subplots(1, 2, figsize=(9, 4), layout="constrained")
axes[0].hist(stats[:,0], bins=np.arange(.5,65.6,2))
axes[1].hist(stats[:,1], bins=np.arange(0,61,2))
axes[0].set(xlabel="Peak infected population", ylabel="Number of outbreaks")
axes[1].set(xlabel="Final recorded time (days)", ylabel="Number of outbreaks")
plt.show()Both histograms have a group of early failures separated from a group of large outbreaks. For this reproducible run the mean peak is 33.73 individuals, and the mean final time is 22.165 days; the largest final time is 56.916 days, so all 1000 runs reached epidemic extinction before the 100-day horizon.
Peak size and extinction time for 1000 original stochastic SIR realizations, with operational failure cutoffs
. Details. 
“Failed to take off” needs an operational definition; the two separated groups in the histograms suggest a low peak or an early end. For example, classify a run as a failure if its peak is below ten infected individuals, or alternatively if its extinction time is below ten days:
failed_by_peak = stats[:,0] < 10
failed_by_time = stats[:,1] < 10
print(failed_by_peak.mean(), failed_by_time.mean())Both definitions give in this run. Here they also classify the same 255 realizations, although agreement is not guaranteed for other simulations or cutoffs. These are empirical failure fractions for the stated rules, not exact model constants.
The early branching process approximation provides a useful check: while susceptible depletion is negligible, gives eventual extinction probability about from one initial infection. Allowing for the initial susceptible fraction gives a rate ratio ; both approximations are compatible with this simulation, but neither defines the finite-population cutoff exactly.
For an undirected graph with adjacency matrix , degrees of vertices and edges, network modularity isEach undirected edge contributes twice to the double sum, explaining the denominator . Include the terms with : even though the observed graph has no self-loops, the null-model contribution on the diagonal is part of this definition.
An independent implementation forms the modularity matrix and sums its within-community blocks. Initializing the best value to negative infinity also handles a candidate collection whose scores are all negative:
def modularity(adjacency,communities):
a=np.asarray(adjacency);degrees=a.sum(axis=1);twom=degrees.sum()
if twom==0:raise ValueError('positive edge count required')
B=a-np.outer(degrees,degrees)/twom
return sum(B[np.ix_(group,group)].sum() for group in communities)/twom
best_score = -np.inf
best_partition = None
for candidate in community_partitions:
value = modularity(adjacency_matrix, candidate)
if value > best_score:
best_score, best_partition = value, candidate
print(best_partition, best_score)The best partition is , with . There are seven edges, and each selected community has three internal edges and total degree seven, so independentlyThe supplied candidates include the optimum, but do not exhaust even the balanced two-community set partitions: there are of those, rather than six. To substantiate a best partition among all community counts, the following short recursion enumerates all 203 unlabeled set partitions of the six vertices:
def partitions(items):
if not items:
yield [];return
first,*rest=items
for part in partitions(rest):
yield [[first]]+[g.copy() for g in part]
for j in range(len(part)):
new=[g.copy() for g in part];new[j]=[first]+new[j];yield new
candidates = list(partitions(list(range(len(adjacency_matrix)))))
best_partition = max(candidates, key=lambda p: modularity(adjacency_matrix, p))
print(len(candidates), best_partition, modularity(adjacency_matrix, best_partition))This exhaustive calculation returns the same two triangles and the same network modularity. Enumerating every set partition becomes impractical rapidly as the number of vertices increases, but is small enough here to verify the global optimum.
Articles by others on the same topic
There are currently no matching articles.