Simulation de QA/QC sur duplicatas géochimiques avec incertitude contrôlée¶
Cette section présente trois graphiques types couramment utilisés dans l’évaluation des duplicatas en contrôle de la qualité (QA/QC).
Nous explorerons ensemble leurs points forts et leurs limites, en mettant en évidence ce qu’ils permettent de détecter — ou non — dans les écarts entre duplicatas.
📊 1. Nuage de points des duplicatas¶
Ce graphique compare directement les deux séries simulées de duplicatas.
- La ligne noire en pointillé représente l’égalité parfaite (Duplicata 1 = Duplicata 2).
- Les bandes colorées indiquent les tolérances acceptables de ±10%, ±20% et ±30%.
- Les points rouges indiquent les cas où l’écart dépasse ±10%, signalant un problème potentiel de reproductibilité.
Cela permet une évaluation visuelle immédiate du respect des critères QA/QC selon les tolérances définies. Cependant, il n’est pas très informatif. On a recourd généralement aux graphiques des points 2 et 3.
📈 2. Différence relative (%) selon la moyenne des duplicatas¶
Ce graphique montre la différence relative entre les deux duplicatas en pourcentage, en fonction de leur moyenne :
- Il met en évidence les écarts systématiques ou aléatoires.
- Les lignes pointillées indiquent les niveaux de tolérance.
- Les points rouges signalent les duplicatas hors tolérance de ±10 %.
Ce graphique permet d’identifier si les écarts entre duplicatas sont constants ou proportionnels à l’intensité des valeurs — un effet souvent appelé effet multiplicatif.
On observe généralement que les faibles teneurs présentent des erreurs relatives plus élevées que les fortes teneurs.
Il est donc essentiel de porter une attention particulière à la zone autour de la teneur de coupure, où ces erreurs peuvent avoir un impact significatif sur les décisions d’exploitation.
📐 3. Courbe HARD (Half Absolute Relative Difference)¶
Le graphique HARD trace la courbe cumulative de l’erreur relative :
- Sur l’axe vertical, on mesure l’écart relatif (|D1 − D2| / (D1 + D2)).
- L’axe horizontal correspond au rang normalisé des points (i.e., leur position dans la distribution triée).
- Le point rouge représente un objectif typique (par ex. 90% des duplicatas dans ±10%).
Ce graphique est souvent utilisé pour évaluer la performance globale du protocole de QA/QC.
Source
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider
def generate_correlated_lognormal_series_mv_variable_noise(
n_points=100,
base_median=2,
sigma_base=0.4,
correlation=0.95,
p=0.0, # nouveau paramètre contrôle la pente du sigma selon la valeur
):
mu = np.log(base_median)
# Tirage base_vals avec sigma fixe non nul (exemple 0.2)
base_vals = np.exp(np.random.normal(mu, 0.2, size=n_points))
safe_vals = np.clip(base_vals, 1e-3, None)
sigma_vals = sigma_base # scalaire
corr_mat = np.array([[1.0, correlation],
[correlation, 1.0]])
dup1 = np.empty(n_points)
dup2 = np.empty(n_points)
for i in range(n_points):
cov = corr_mat * sigma_vals**2 # <-- ici on utilise sigma_vals comme scalaire
noise = np.random.multivariate_normal(mean=[0, 0], cov=cov)
dup1[i] = np.exp(mu + noise[0])
dup2[i] = np.exp(mu + noise[1])
dup1 += np.random.normal(0, p, size=n_points)
dup2 += np.random.normal(0, p, size=n_points)
return dup1, dup2
def plot_lognormal_variable_noise(median=2.0, sigma=0.4, corr=0.9, p=0.0):
np.random.seed(42)
dup1, dup2 = generate_correlated_lognormal_series_mv_variable_noise(
n_points=200,
base_median=median,
sigma_base=sigma,
correlation=corr,
p=p,
)
fig = plt.figure(figsize=(18, 6), constrained_layout=True)
max_val = min(max(np.quantile(dup1, 0.95), np.quantile(dup2, 0.95)) * 1.1, 100)
# --- Scatter duplicatas ---
ax1 = fig.add_subplot(1, 3, 1)
ax1.scatter(dup1, dup2, alpha=0.6, label="Points")
lims = [0, max_val]
ax1.plot(lims, lims, 'k--', label="y = x")
tolerances = [0.1, 0.2, 0.3]
colors = ['r', 'orange', 'purple']
counts_out = []
for tol, col in zip(tolerances, colors):
lower = 1 - tol
upper = 1 + tol
ax1.plot(lims, [lims[0]*upper, lims[1]*upper], color=col, linestyle='-', alpha=0.6, label=f"±{int(tol*100)}%")
ax1.plot(lims, [lims[0]*lower, lims[1]*lower], color=col, linestyle='-', alpha=0.6)
out_of_bounds = (dup2 < dup1 * lower) | (dup2 > dup1 * upper)
counts_out.append(np.sum(out_of_bounds))
tol_max = 0.1
lower_max = 1 - tol_max
upper_max = 1 + tol_max
out_max = (dup2 < dup1 * lower_max) | (dup2 > dup1 * upper_max)
ax1.scatter(dup1[out_max], dup2[out_max], color='red', s=80, label='Hors ±10%')
median_r = round(median, 2)
sigma_r = round(sigma, 2)
corr_r = round(corr, 2)
p_r = round(p, 2)
ax1.set_xlabel("Duplicata 1")
ax1.set_ylabel("Duplicata 2")
ax1.set_title(f"Séries lognormales corrélées\nMédiane={median_r}, Sigma={sigma_r}, Corr={corr_r}, p={p_r}\n"
f"Hors ±10%: {counts_out[0]} | ±20%: {counts_out[1]} | ±30%: {counts_out[2]} sur {len(dup1)} points")
ax1.grid(True)
ax1.legend(loc='best')
ax1.set_xlim(0, max_val)
ax1.set_ylim(0, max_val)
# --- Différence relative ---
ax2 = fig.add_subplot(1, 3, 2)
mean_vals = (dup1 + dup2) / 2
diff_rel = 100 * (dup1 - dup2) / mean_vals
ax2.scatter(mean_vals, diff_rel, alpha=0.6, color='blue', label="Différence relative (%)")
for tol, col in zip(tolerances, colors):
ax2.axhline(y=tol*100, color=col, linestyle='--', alpha=0.6, label=f'±{int(tol*100)}%')
ax2.axhline(y=-tol*100, color=col, linestyle='--', alpha=0.6)
out_diff_max = (diff_rel > 10) | (diff_rel < -10)
ax2.scatter(mean_vals[out_diff_max], diff_rel[out_diff_max], color='red', s=80, label='Hors ±10%')
ax2.set_xlabel("(Duplicata 1 + Duplicata 2) / 2")
ax2.set_ylabel("Différence relative (%)")
ax2.set_title("Différence relative entre duplicatas")
ax2.grid(True)
ax2.legend(loc='best')
ax2.set_xlim(0, max_val)
ax2.set_ylim(-50, 50)
# --- HARD plot ---
ax3 = fig.add_subplot(1, 3, 3)
N = len(dup1)
hard_vals = np.sort(np.abs(dup1 - dup2) / (dup1 + dup2))
ranks = np.arange(1, N + 1) / (N + 1)
ax3.plot(ranks, hard_vals, color='black', linewidth=2, label='Cible')
ax3.plot(0.9, 0.1, 'o', color='red', markersize=10, label='Point critique')
ax3.set_xlabel('Rang/(N+1)')
ax3.set_ylabel('Graphique HARD')
ax3.set_ylim(0, max(0.3, hard_vals.max() * 1.1)) # plus flexible
ax3.grid(True)
ax3.legend(loc='best')
plt.show()
interact(
plot_lognormal_variable_noise,
median=FloatSlider(min=0.1, max=10, step=0.001, value=0.9, description="Médiane"),
sigma=FloatSlider(min=0.05, max=1.0, step=0.001, value=1.4, description="Sigma base"),
corr=FloatSlider(min=0.95, max=0.999, step=0.001, value=0.996, description="Corrélation", readout_format=".3f"),
p=FloatSlider(min=0, max=2.0, step=0.05, value=0, description="Bruit décroissant p")
)Loading...