El programa
El cálculo de O'Neill tiene dos mitades. Primero, un ajuste por mínimos cuadrados: se busca la curva suave (un polinomio) que mejor pasa por los 13 puntos del patrón. Después, una reorganización de ese mismo polinomio en piezas independientes —los polinomios ortogonales—, cuyos coeficientes son los λ. Cada lección añade una idea matemática nueva y la aplica a una muestra real de la hoja de cálculo de O'Neill.
El color te permitirá seguir con la vista en qué parte de cada fórmula entra el radio, que es una de tus preguntas. La otra pestaña, Mapa del radio, lo resume en un solo dibujo.
-
1
Disponible
Preparar los datos: de las ppm a 13 puntos
Normalizar al condrito CI, tomar logaritmos, elegir el radio como eje X y apartar el Eu.
-
2
Disponible
¿Qué es «ajustar»? Mínimos cuadrados con una recta
Residuos, por qué se suman sus cuadrados (con una analogía de muelles) y la receta de la recta. Medir el radio desde su media convierte la ordenada en λ0 y la pendiente en λ1. Incluye la primera mitad de la Fig. 1a de O'Neill.
-
3
Disponible
De la recta a la curva: el polinomio y su sistema de ecuaciones
Qué forma aporta cada término (r², r³, r⁴), de dónde salen 5 ecuaciones con 5 incógnitas y cómo las resuelve Excel con MINVERSA y MMULT, exactamente como la macro de O'Neill.
-
4
Disponible
Por qué los coeficientes «normales» a0…a4 no sirven para comparar
Con la muestra guía: a0 salta de 3.19 a −2367 al añadir términos, y los coeficientes están correlacionados al 99.99 %. Una simulación de análisis repetidos y la segunda mitad de la Fig. 1a de O'Neill (los patrones D y D†), explicada.
-
5
Disponible
Los polinomios ortogonales: fabricarlos con los radios
Qué significa «ortogonal» (piezas que no se pisan) y cómo salen, solo de los 13 radios, las constantes 1.05477, 1.00533, 1.12824… de la Tabla 2 (Bevington & Robinson, §7.3).
-
6
Disponible
Calcular los λ: tres caminos, un resultado
La fórmula directa λn = Σ y·fn / Σ fn², el sistema que resuelve la macro y la conversión desde un ajuste normal (ESTIMACION.LINEAL). Qué significa cada λ y cómo leer el diagrama λ2–λ1.
-
7
Disponible
¿Es bueno el ajuste? Residuos, Eu/Eu*, χ² e incertidumbres
Cocientes observado/calculado, la anomalía de Eu calculada en r = 1.066 Å, el χ² reducido, el s(ln[REE])* inferido y los factores ζ que dan el error de cada λ.
-
8
Disponible
Opcional · Hazlo tú con tus datos
El cálculo completo en Excel o en Python, y qué hacer si te falta un elemento o quieres dejar fuera otro (por ejemplo el Ce, algo habitual en circón): mantener los polinomios de la Tabla 2, como recomienda O'Neill, o recalcularlos, y qué cambia en cada caso.
-
9
Disponible
Opcional · Vectores petrogenéticos ψ
Por qué λ = λ⁰ + ψ: la fusión y la cristalización se describen con los mismos polinomios y se suman como flechas en el diagrama λ2–λ1.
Qué vas a necesitar
Nada de cálculo diferencial. Sumar, multiplicar, dividir y, en la Lección 1, el logaritmo neperiano, que explico desde cero. Donde el libro de Bevington usa derivadas o matrices, te lo contaré con una analogía y te daré el resultado comprobado con números.
Materiales
- O'Neill (2016), pp. 1464–1467: el método, las Tablas 1 y 2 y la Fig. 1.
- Apéndice 1 del material suplementario: por qué logaritmos y cómo se evalúan las incertidumbres (ecs. A1–A4).
- La hoja
LstSq_REE_OrthogPolynomial_April2015.xlsy su macroREE_OrthogRegress_5term. La he extraído para ver exactamente qué calcula; iré mostrando sus líneas en cada lección. REEshapes_reference_datasets_for_plots.xls: los 616 vidrios de Jenner & O'Neill (2012b) y los basaltos HSDP2 de Hawái.- Bevington & Robinson (2003), 3.ª ed.: cap. 6 (la recta), cap. 7 §7.1–7.3 (polinomios, matrices y «Independent parameters», pp. 116–131), §7.4 (logaritmos) y cap. 11 (χ²).
Curso completo. Las siete lecciones y las dos opcionales están en esta página. Si algo no queda claro, o quieres aplicarlo a tus propios datos, pregúntame.
Dónde entra el radio iónico
El radio aparece en seis lugares del cálculo, y en un séptimo si seguimos hasta los modelos petrogenéticos. La clave está en este dibujo: todo lo que depende solo de los radios se calcula una vez y vale para todas las muestras. Lo único que cambia de un análisis a otro son los valores y.
- 1Eje X. Cada yi se dibuja y se ajusta frente a su radio ri. Es el único papel del radio que se ve a simple vista. Lección 1
- 2Constantes de los polinomios ortogonales. 1.05477 es la media de los 13 radios (sin el Eu); 1.00533 y 1.12824, 0.99141… salen de exigir que cada pieza del polinomio sea independiente de las demás. Solo dependen de los radios: por eso la Tabla 2 sirve para cualquier muestra. Lección 5
- 3Tabla de valores fn(ri). Cada radio se sustituye en f1…f4: son las columnas f de la Tabla 1 de O'Neill. Por ejemplo, f1(La) = 1.160 − 1.05477 = 0.1052. Lección 5
- 4El ajuste. La macro suma rik·fn(ri) (una tabla de 5 × 5 que solo depende de los radios) y rik·yi (radios por datos), y resuelve el sistema. Con la fórmula directa: λn = Σ yi fn(ri) / Σ fn(ri)². Lecciones 3 y 6
- 5La curva calculada. Se evalúa en los 14 radios. En el del Eu (1.066 Å) da el Eu* con el que se calcula Eu/Eu*, aunque el Eu no haya entrado en el ajuste. Lección 7
- 6Incertidumbres. ζn = 1/√Σ fn(ri)² depende solo de los radios y convierte el error analítico en el error de cada λ: para un 1 %, s(λ) = 0.0028, 0.046, 0.99, 19.8 y 424 (Tabla 2). Lección 7
- 7Extensión. Los coeficientes de reparto mineral/fundido también son funciones suaves de r y se ajustan con los mismos polinomios: así nacen los vectores de proceso ψ. Lección 9
Dónde no entra el radio: en la normalización al condrito, en el logaritmo ni en la definición del χ² (solo indirectamente, a través de la curva calculada).
Lección 1 · de 7 (+ 2 opcionales)
Preparar los datos: de las ppm a 13 puntos
Antes de ajustar nada, O'Neill convierte cada análisis en 13 parejas de números (r, y). Aquí vemos de dónde sale cada una: la y, de la concentración; la r, del tamaño del ion.
- Idea matemática nueva
- Dividir y el logaritmo neperiano (ln)
- Dónde entra el radio
- ① Como eje X de cada punto
- Muestra guía
- Vidrio VG 189, dorsal Mesoatlántica
La muestra guía
Seguiremos un vidrio basáltico real: VG 189, de la dorsal Mesoatlántica (21.9° S, 7.98 % MgO). Es la primera fila de la hoja de cálculo de O'Neill, que reúne los 616 vidrios de Jenner & O'Neill (2012b). En la figura puedes cambiar a otras tres muestras que O'Neill usa en su Fig. 17d: una empobrecida (VG 304), una enriquecida (VG 4080) y un basalto de Hawái (8-2.70).
Paso 1 · Dividir por el condrito CI
Si dibujas las ppm tal cual (botón 1 · ppm), sale un zigzag: los elementos de número atómico par (Ce, Nd, Sm, Gd…) son siempre más abundantes que sus vecinos impares. No es geología: es cómo se fabricaron los núcleos en las estrellas (efecto Oddo–Harkins), y los condritos lo tienen igual.
Por eso se divide cada REE por su concentración en el condrito CI (columna «CI» de la Tabla 1 de O'Neill). El zigzag se cancela (botón 2 · ÷ CI) y queda lo que interesa: cuántas veces más —o menos— concentrado está cada elemento que en el material del que se formó la Tierra.
Paso 2 · Tomar el logaritmo neperiano (ln)
El logaritmo convierte multiplicar en sumar. Con el neperiano (ln), multiplicar por 2 equivale siempre a sumar 0.69; multiplicar por 10, a sumar 2.30; y multiplicar por e ≈ 2.718, a sumar exactamente 1.
| Veces el condrito | 1 | 2 | 2.718 | 10 | 20 | 100 |
|---|---|---|---|---|---|---|
| ln | 0 | 0.69 | 1 | 2.30 | 3.00 | 4.61 |
Es la misma idea que el pH o la magnitud de un terremoto: una escala de «veces». Los diagramas habituales usan log10; ln y log10 solo se diferencian en un factor fijo (ln = 2.303 × log10), como centímetros y pulgadas, así que la forma del patrón no cambia.
¿Por qué usa O'Neill el ln y no las concentraciones normalizadas tal cual? Lo razona en el Apéndice 1 del material suplementario:
- Los errores analíticos son porcentajes. Un ±5 % es ±1.4 en un valor de 28, pero ±0.7 en uno de 14. En ln, ese mismo ±5 % mide ≈ ±0.05 en todos los elementos. Activa las barras de error y pasa de 2 a 3: en ln todas miden igual. Así el ajuste puede dar a los 13 puntos el mismo peso, y el cálculo se simplifica muchísimo. (Bevington & Robinson, §7.4: al tomar ln, la incertidumbre de cada punto se convierte en su error relativo, σ/y.)
- El ajuste se vuelve un problema «lineal»: se resuelve con una receta directa de sumas y divisiones (lecciones 2 y 3), sin tanteos. Con las concentraciones sin logaritmo y errores porcentuales sería un problema no lineal (ecs. A3a y A3b del Apéndice).
- Una concentración nunca puede ser negativa; su ln, sí. Por eso es razonable suponer que los ln de medidas repetidas se reparten como una campana de Gauss, que es justo lo que supone el método de mínimos cuadrados. Las concentraciones, en rigor, no pueden hacerlo.
Y un regalo para más adelante: la fusión y la cristalización multiplican las concentraciones por factores; en ln esos factores se suman. Es lo que permitirá sumar «vectores de proceso» en la Lección 9.
Paso 3 · El eje X: radio iónico en vez de número atómico
El diagrama clásico coloca las REE por orden de número atómico Z, a distancias iguales. O'Neill usa otra regla: el radio iónico del catión 3+ en coordinación 8 (Shannon, 1976), la columna r de su Tabla 1. Son los números verdes de las fichas de arriba.
¿Por qué? Porque la química de las REE depende sobre todo de su tamaño. Los electrones 4f que se van añadiendo del La al Lu quedan escondidos bajo las capas 5s y 5p y no participan en los enlaces; lo único que cambia de un elemento al siguiente es que el ion se encoge (contracción lantánida). Además, los coeficientes de reparto mineral/fundido se describen en función del radio (modelo de deformación de red, lattice strain).
La contracción no es uniforme: se va frenando. Del La al Ce el radio baja 0.017 Å; del Yb al Lu, solo 0.008 Å. O'Neill la resume con una parábola:
Ojo con la orientación. En las figuras de O'Neill el radio crece hacia la derecha, así que el La queda a la derecha y el Lu a la izquierda: es un espejo del diagrama clásico. Un patrón enriquecido en REE ligeras «sube» hacia la derecha y tiene pendiente positiva; en la Lección 6 veremos que eso es λ1 > 0.
Para curiosos: ¿de verdad se ajusta mejor con el radio?
En los 616 vidrios de Jenner & O'Neill, el error típico que queda tras el ajuste baja así al usar 2, 3, 4 y 5 términos (O'Neill, p. 1471): con el radio, 11.7 → 3.9 → 2.4 → 2.14 %; con Z, 14.0 → 5.6 → 3.09 → 2.56 %. Con el radio hace falta un término menos para la misma calidad.
He repetido el cálculo con su hoja de datos y reproduzco esas cifras. Un matiz: las de Z salen colocando las 14 REE equiespaciadas, como en el diagrama clásico, sin hueco para el Pm. Con el Z real, dejando el hueco, la ventaja del radio es algo menor (13.2, 5.3, 2.7 y 2.16 %), pero se mantiene con 2 y 3 términos, que son los que importan en basaltos.
Paso 4 · Apartar el Eu (y no echar de menos al Pm)
El europio puede estar como Eu²⁺, que la plagioclasa acepta mucho mejor que el Eu³⁺. Eso produce anomalías que no forman parte de la forma suave del patrón, así que O'Neill excluye el Eu del ajuste y lo recupera al final para calcular Eu/Eu* (Lección 7). En la figura es el círculo hueco.
El Ce también puede tener anomalías (Ce⁴⁺), pero en basaltos O'Neill lo mantiene. El prometio no se analiza nunca porque no tiene isótopos estables, y como el eje es el radio y no una posición «de lista», su hueco no molesta.
Resultado: 13 puntos, cada uno con su r y su y.
El resultado: 13 parejas (r, y)
Esta tabla es todo lo que entra en el ajuste. Cambia con la muestra elegida en la figura (VG 189).
| REE | Z | ppm | CI (ppm) | [REE]/CI | y = ln(…) | r (Å) |
|---|---|---|---|---|---|---|
| Media de las 13 y (sin Eu) | ||||||
Tu primer λ. λ0 es simplemente la media de las 13 y (Tabla 2 de O'Neill: «average ln([REE]/[REE]CI)»). Como e3.19 ≈ 24, las REE de esta muestra están, de media, unas 24 veces por encima del condrito.
La hoja de O'Neill da λ0 = 3.19 para esta muestra. Coincide.
Lo mismo, dentro de la macro de O'Neill
La hoja LstSq_REE_OrthogPolynomial_April2015.xls hace todo esto con una macro de Visual Basic, REE_OrthogRegress_5term. Estas son las líneas que corresponden a esta lección (extracto; los comentarios en español son míos):
' Radios iónicos en Å, Shannon (1976) ← aquí entran los radios en la macro
RR(1) = 1.16 ' La
RR(2) = 1.143 ' Ce … y así hasta …
RR(14) = 0.977 ' Lu
' Condrito CI en ppm (Tabla 1)
CI(1) = 0.2472 ' La … y así hasta …
CI(14) = 0.0251 ' Lu
' Para cada REE (M = 1 a 14) de cada muestra:
If DataIn.Cells(Sample%, M%) <= 0 Then GoTo Skip ' vacío o ≤ 0: se salta
X(M%) = DataIn.Cells(Sample%, M%)
If Normalize = 1 Then X(M%) = Log(X(M%) / CI(M%)) Else X(M%) = Log(X(M%))
' ÷ CI y ln (con B3 = 0, solo ln)
If M% = 6 Then GoTo Skip ' el 6.º (Eu) no entra en el ajusteTrampa frecuente: en VBA, Log() es el logaritmo neperiano; en una celda de Excel, LOG() es en base 10 y el neperiano es LN(). Y la celda B3 de la hoja decide si se normaliza: 1 para concentraciones, 0 para coeficientes de reparto (solo ln, sin dividir).
Practica
1. Calcula la y del Yb de VG 189: 3.996 ppm, con CI = 0.1684 ppm.
2. VG 304 tiene λ0 = 2.65. ¿Cuántas veces el condrito es su concentración media de REE?
Resumen
- ÷ CI quita el zigzag de la nucleosíntesis y deja el enriquecimiento de cada REE respecto al condrito.
- ln convierte errores porcentuales en errores iguales para todos los elementos: los 13 puntos pesan lo mismo y el ajuste es lineal.
- Eje X = radio iónico: la variable que controla la química de las REE. En las figuras de O'Neill crece hacia la derecha, con el La a la derecha.
- Sin Eu: quedan 13 puntos (ri, yi). Su media ya es el primer coeficiente, λ0.
Dónde ha aparecido el radio en esta lección: solo en un sitio (① del mapa), como coordenada x de cada punto. Los y no dependen del radio en absoluto.
Próxima lección · ¿Qué es «ajustar»? Buscaremos la recta que mejor pasa por estos 13 puntos, veremos por qué se suman los cuadrados de las distancias (con una analogía de muelles) y descubriremos que, si medimos el radio desde su valor medio (1.0548 Å), la ordenada en el origen es exactamente λ0.
Fuentes de esta lección
- O'Neill, H. St. C. (2016). The smoothness and shapes of chondrite-normalized rare earth element patterns in basalts. Journal of Petrology 57, 1463–1508: pp. 1464–1465 (contracción lantánida, r frente a Z, Eu y Ce), Tabla 1, p. 1471 (r frente a Z).
- O'Neill (2016), Apéndice 1 del material suplementario, ecs. A1–A3 (por qué logaritmos).
- Macro
REE_OrthogRegress_5termde la hojaLstSq_REE_OrthogPolynomial_April2015.xls; datos de Jenner & O'Neill (2012b) y Feigenson et al. (2003) deREEshapes_reference_datasets_for_plots.xls. - Bevington, P. R. & Robinson, D. K. (2003). Data Reduction and Error Analysis for the Physical Sciences, 3.ª ed., §7.4, ecs. 7.49–7.50.
- Shannon, R. D. (1976). Acta Crystallographica A32, 751–767 (radios iónicos).
Lección 2 · de 7 (+ 2 opcionales)
¿Qué es «ajustar»? Mínimos cuadrados con una recta
Tenemos 13 puntos que no están alineados. ¿Cuál es la recta que mejor los representa? Para contestar hace falta una regla que puntúe cada recta posible. Esa regla es la suma de cuadrados, y al aplicarla aparecerán, casi sin querer, los dos primeros λ de O'Neill.
- Idea matemática nueva
- Residuo, sumatorio Σ y suma de cuadrados (χ²)
- Dónde entra el radio
- En las sumas Σr, Σr² y Σr·y; al centrar, en r̄ (②) y en r − r̄ (③)
- Muestra guía
- VG 189: los 13 puntos de la Lección 1
El residuo: cuánto falla la recta en cada punto
Dibuja cualquier recta sobre el patrón. En cada REE, la distancia vertical entre el punto y la recta es el residuo:
Es positivo si el punto queda por encima de la recta y negativo si queda por debajo. En la figura, cada residuo es un muelle; debajo tienes los mismos residuos en forma de barras.
Para puntuar una recta hay que resumir sus 13 residuos en un solo número, y para eso hay que sumar. La letra griega Σ (sigma) significa «suma para todos los puntos»: Σy es yLa + yCe + … + yLu, sin el Eu. Para VG 189, Σy = 41.437.
¿Por qué al cuadrado?
La suma de los residuos tal cual no sirve: los positivos y los negativos se compensan. De hecho, cualquier recta que pase por el centro de los datos (r̄, ȳ) tiene suma de residuos cero, esté inclinada como esté. Pruébalo: pulsa «Soltar los muelles», deja el pivote en r̄ y mueve solo la pendiente: la fuerza neta sigue en cero aunque la recta sea malísima.
Por eso se suman los cuadrados de los residuos (Bevington & Robinson, §6.2):
- Un cuadrado siempre es positivo, así que ya no hay compensaciones.
- Castiga más los fallos grandes: fallar por 0.2 cuesta cuatro veces lo que fallar por 0.1.
- Si los errores se reparten como una campana de Gauss, la recta con la menor suma de cuadrados es la más probable (máxima verosimilitud, ecs. 6.5–6.9).
- Lleva a una receta directa, sin tanteos. Con valores absolutos no la habría (p. 103).
Si antes de elevar al cuadrado divides cada residuo por el error analítico s, la suma se llama χ² (ji cuadrado): χ² = Σ (residuo / s)². Gracias al ln de la Lección 1, s es el mismo para todos los puntos, así que minimizar χ² y minimizar S es lo mismo: s solo cambia la escala (ec. A4 del Apéndice de O'Neill). El χ² nos servirá en la Lección 7 para juzgar si un ajuste es bueno.
La analogía de los muelles
En reposo, la varilla está en equilibrio, y eso impone dos condiciones. Míralas en los contadores de la figura, que se ponen a cero al soltar:
- No sube ni baja: los muelles que tiran hacia arriba compensan a los que tiran hacia abajo → Σ residuos = 0 (fuerza neta).
- No gira: los tirones a un lado y a otro se compensan, como en un balancín → Σ ri·residuoi = 0 (giro neto).
Escritas con sumas, son dos ecuaciones con dos incógnitas: la ordenada a y la pendiente b de la recta y = a + b·r. Se llaman ecuaciones normales (Bevington, ec. 6.11):
Σry = a·Σr + b·Σr2
La receta
Resolviendo ese par de ecuaciones se obtiene la receta de la recta cuando todos los puntos tienen el mismo error (Bevington, ec. 6.13):
b = (N·Σry − Σr·Σy) / Δ
a = (Σr2·Σy − Σr·Σry) / Δ
Con VG 189 (N = 13), las cuatro sumas son:
| Σr | Σr² | Σy | Σr·y |
|---|---|---|---|
| 13.712 | 14.510216 | 41.437392 | 43.628696 |
Dos avisos. Δ es una resta de dos números casi iguales (188.63 − 188.02): si redondeas Σr y Σr² a dos decimales, Δ sale 0.666 en vez de 0.614, un 8 % de error por un redondeo inocente. Hay que llevar todos los decimales (Bevington, p. 111). Y a = 4.934 es el valor de la recta en r = 0: el ln de la concentración de un ion de tamaño cero. No significa nada geológicamente.
Fíjate también en que N, Σr, Σr² y Δ solo dependen de los radios: valen lo mismo para las cuatro muestras de la figura y para las 616 de la hoja de O'Neill.
El truco: medir el radio desde su media
Cambiemos de regla: en vez de medir el radio desde cero, midámoslo desde su valor medio, r̄ = Σr / 13 = 1.054769 Å. La nueva variable es x′ = r − r̄. Por definición de media, los radios mayores que r̄ compensan exactamente a los menores, así que Σx′ = 0. Con eso, las dos ecuaciones normales se separan y cada incógnita se calcula sola:
¡Son los dos primeros coeficientes de O'Neill! La recta centrada es y = λ0 + λ1·(r − 1.05477), y r − 1.05477 es exactamente el primer polinomio ortogonal de su Tabla 2, f1. La pendiente es la misma en las dos versiones; lo que cambia es el significado de la ordenada: a′ es el nivel medio del patrón.
La trampa del origen
O'Neill lo ilustra en su Fig. 1a con cuatro patrones inventados: A = 5 − 3r, B = 3, C = 4 y D = 2 + 3r. Si comparas sus ordenadas en el origen (5, 3, 4 y 2), dirías que A es el más rico en REE… y es el más pobre. Con el radio centrado, las ordenadas son 1.84, 3, 4 y 5.16: el orden real.
Pruébalo: con el pivote en r = 0, mueve un poco la pendiente y verás que la recta sale disparada del gráfico. Encontrar la mejor recta así exige ir y volver entre los dos deslizadores. Con el pivote en r̄, basta con ajustar cada uno por separado. En la Lección 4 veremos que esto mismo, multiplicado por cinco términos, es el gran problema de los coeficientes «normales».
¿Es buena la recta?
Suelta los muelles con VG 189 y mira las barras de residuos: son negativas en los extremos (La −0.33, Lu −0.19) y positivas en el centro (Sm +0.22). No es ruido al azar: es una joroba, la curvatura del patrón que una recta no puede seguir.
En números: el residuo típico es 0.16 en ln (O'Neill lo expresaría como un 16 %), cuando el error analítico ronda el 2 %. Con s = 2 %, χ² = 705, y un buen ajuste debería dar un valor parecido al número de grados de libertad: 13 puntos − 2 parámetros = 11. La recta se queda muy corta: hace falta curvatura, y eso es la Lección 3.
Compruébalo en Excel
Con las 13 y (sin Eu) en B2:B14 y sus radios en C2:C14. En Excel en español, los argumentos se separan con punto y coma y los decimales llevan coma:
=PENDIENTE(B2:B14;C2:C14) → −1,656 pendiente b = λ₁
=INTERSECCION.EJE(B2:B14;C2:C14) → 4,934 la recta en r = 0
=PROMEDIO(B2:B14) → 3,187 ȳ = λ₀
=PROMEDIO(C2:C14) → 1,054769 r̄Y si en D2:D14 pones el radio centrado (=C2-PROMEDIO($C$2:$C$14)), =INTERSECCION.EJE(B2:B14;D2:D14) da 3,187: la ordenada se ha convertido en λ0.
¿Y la macro de O'Neill?
La macro no ajusta rectas: ajusta directamente cinco términos. Pero sus dos primeros coeficientes para VG 189, λ0 = 3.19 y λ1 = −1.66, son exactamente los que acabas de calcular con la recta centrada. Que no cambien al añadir más términos es la gracia de los polinomios ortogonales (Lección 5).
Practica
1. Para VG 4080, Σx′y = 0.52063 y Σx′² = 0.04722. ¿Cuánto vale λ1?
2. Una recta pasa exactamente por el centro de los datos (r̄, ȳ), con una pendiente cualquiera. ¿Cuánto vale la suma de sus residuos, sin elevarlos al cuadrado?
3. Dos basaltos tienen exactamente la misma forma de patrón, pero uno tiene todas las REE al doble de concentración que el otro. ¿Qué cambia?
Resumen
- Residuo: distancia vertical entre el punto y la recta (observado − calculado).
- Mínimos cuadrados: la mejor recta es la que hace mínima la suma de los residuos al cuadrado (el χ², si antes se dividen por el error). Los cuadrados evitan que se compensen positivos y negativos.
- Muelles: en reposo, Σ residuos = 0 y Σ r·residuo = 0. Son las dos ecuaciones normales, que dan una receta directa.
- Centrar el radio en r̄ = 1.05477 separa las incógnitas: la ordenada es la media (λ0) y la pendiente es λ1. La ordenada en r = 0 no significa nada.
- Para VG 189, la recta deja una joroba en los residuos: hace falta curvatura.
Dónde ha aparecido el radio en esta lección: en tres sitios. Primero, en las sumas Σr, Σr² y Σr·y de la receta. Segundo, en r̄ = 1.05477, la media de los radios, que es la primera constante de los polinomios ortogonales (② del mapa). Tercero, en r − r̄, que es el primer polinomio ortogonal, f1 (③ del mapa).
Próxima lección · De la recta a la curva. Añadiremos términos en r², r³ y r⁴ para seguir la joroba. Las dos ecuaciones normales se convertirán en cinco, que Excel resuelve con MINVERSA y MMULT, igual que la macro de O'Neill. Y veremos que la tabla de ese sistema depende solo de los radios.
Fuentes de esta lección
- Bevington & Robinson (2003), §6.2 (ec. 6.4 y p. 103: por qué los cuadrados; máxima verosimilitud, ecs. 6.5–6.9), §6.3 (ecs. 6.10–6.13), Tabla 6.1 (un ejemplo resuelto con las cuatro sumas), p. 111 (errores de redondeo) y ejercicio 6.7 (al centrar la variable, la ordenada es la media).
- O'Neill (2016), Fig. 1a y su pie (pp. 1465–1466) y Tabla 2 (f1 = rREE − 1.05477).
- O'Neill (2016), Apéndice 1, ecs. A3b y A4 (el χ² que se minimiza).
- Datos: vidrios VG 189, VG 304 y VG 4080 (Jenner & O'Neill, 2012b) y basalto 8-2.70 (Feigenson et al., 2003), de las hojas de O'Neill.
Lección 3 · de 7 (+ 2 opcionales)
De la recta a la curva: el polinomio y su sistema de ecuaciones
La recta dejaba una joroba en los residuos. Para seguirla hay que dejar que la línea se doble, y eso es un polinomio. El criterio sigue siendo el mismo (mínimos cuadrados), pero ahora hay cinco incógnitas y un sistema de cinco ecuaciones. Resolverlo es exactamente lo que hace la macro de O'Neill con dos funciones de Excel.
- Idea matemática nueva
- El polinomio y un sistema de ecuaciones escrito como tabla (matriz)
- Dónde entra el radio
- ④ La tabla del sistema solo contiene sumas de potencias de r
- Muestra guía
- VG 189, de la recta a la cuártica
Una línea que puede doblarse
Un polinomio es una suma de potencias de la variable, cada una multiplicada por su coeficiente:
Cada término que añades le da a la línea una posibilidad más de doblarse:
El mismo criterio, más incógnitas
Buscamos los coeficientes que hacen mínima la suma de los residuos al cuadrado, igual que con la recta. En la analogía de los muelles, la varilla rígida se convierte en una regla flexible. Al soltarla, se queda en reposo cuando cada forma de moverse está equilibrada, y hay una condición por término:
Σ r·residuo = 0 (inclinar)
Σ r2·residuo = 0 (doblar una vez)
Σ r3·residuo = 0
Σ r4·residuo = 0
Dicho de otra forma: en los residuos no puede quedar nada que alguna de las piezas 1, r, r², r³ o r⁴ sea capaz de explicar. Si quedara, el ajuste aún podría mejorar.
Sustituyendo residuo = y − (a0 + a1r + …) y agrupando, salen las cinco ecuaciones normales (Bevington, ec. 7.7). La primera es:
Las otras cuatro son iguales, pero multiplicando todo por r, r2, r3 y r4. En total aparecen las sumas de potencias del radio desde Σr0 = N hasta Σr8.
El sistema como una tabla (matriz)
Escribir cinco ecuaciones tan largas es incómodo, así que se ordenan en una tabla (Bevington, §7.2, ecs. 7.13–7.15): a la izquierda, los números que multiplican a cada incógnita (la matriz α); en medio, las incógnitas (a); a la derecha, los resultados (el vector β).
La tabla sigue al laboratorio: cambia el número de términos o la muestra y mira qué parte cambia. Pasa el ratón por una celda de α para ver sus gemelas.
Fíjate en tres cosas:
- α solo contiene sumas de potencias de los radios: en la fila j y la columna k está Σrj+k. No depende de la muestra: es la misma tabla para VG 189, para VG 304 y para los 616 vidrios. Es la marca ④ del mapa.
- β mezcla radios y datos: Σy, Σr·y, Σr2·y… Es lo único que cambia de una muestra a otra.
- Los números se repiten en cada diagonal que baja hacia la izquierda, porque en ella j + k no cambia.
Resolver: la tabla inversa
Con una sola ecuación, 5·a = 10, despejarías dividiendo: a = 10 / 5, que es lo mismo que multiplicar por el inverso, a = (1/5)·10. Con una tabla se hace igual: se calcula la tabla inversa α−1 y se multiplica por β (Bevington, ec. 7.19):
Nadie calcula a mano la inversa de una tabla de 5 × 5. Excel lo hace con MINVERSA (inversa) y MMULT (multiplicar tablas), y eso es exactamente lo que hace la macro de O'Neill, con los nombres en inglés MInverse y MMult.
Para VG 189 con cinco términos sale:
¡Números enormes, con signos alternos, para describir valores de y en torno a 3! En r̄ los cinco términos valen −2367, +9555, −14479, +9771 y −2478, y su suma es 3.35: gigantes que casi se anulan. Además, la tabla α es muy delicada, porque sus columnas se parecen muchísimo. Resolviéndola directamente con los 16 dígitos de un ordenador, a0 sale −2366.71; con aritmética exacta, −2366.66. Estos coeficientes son correctos como conjunto, pero uno a uno no significan nada. Es el problema de la Lección 4.
¿Cuántos términos?
Pulsa 1, 2, 3, 4 y 5 en la figura. Con VG 189:
| Términos | Residuo típico s | χ² reducido (s = 2 %) |
|---|---|---|
| 1 | 0.185 | 86 |
| 2 | 0.160 | 64 |
| 3 | 0.043 | 4.7 |
| 4 | 0.025 | 1.5 |
| 5 | 0.020 | 1.0 |
Con tres términos desaparece la joroba. Con cinco, el residuo típico (0.020, un 2 %) es del tamaño del error analítico y el χ² reducido vale ≈ 1: ya no queda nada que explicar salvo el ruido de la medida. En los 616 vidrios de Jenner & O'Neill, el término lineal explica el 83.4 % de la variación del patrón (aparte del nivel medio); con la curvatura se llega al 98.2 %, y los dos términos siguientes solo suben al 99.3 % y al 99.4 % (O'Neill, pp. 1470–1471). Por eso tres términos bastan para casi cualquier basalto.
Cada coeficiente «gasta» un punto: con 13 REE y 5 términos quedan 13 − 5 = 8 grados de libertad (ν). Con 13 términos, el polinomio pasaría exactamente por los 13 puntos… y dibujaría también el ruido analítico.
Dentro de la macro de O'Neill
La macro construye su tabla de 5 × 5 (la llama A) y su vector (Z) sumando REE a REE, y después resuelve. Extracto; los comentarios en español son míos y los radios (RR) van en verde azulado:
For M% = 1 To NoData% ' recorre las 14 REE (el Eu se salta)
A(1, 1) = A(1, 1) + 1 ' fila 1: Σ 1
A(1, 2) = A(1, 2) + (RR(M%) - 1.054769) ' fila 1: Σ f1(r)
' …
A(2, 1) = A(2, 1) + RR(M%) ' fila 2: Σ r
A(2, 2) = A(2, 2) + RR(M%) * (RR(M%) - 1.054769) ' fila 2: Σ r·f1(r)
' … así hasta A(5, 5)
Z(1, 1) = Z(1, 1) + X(M%) ' Σ y
Z(2, 1) = Z(2, 1) + X(M%) * RR(M%) ' Σ r·y
Z(3, 1) = Z(3, 1) + X(M%) * RR(M%) * RR(M%) ' Σ r²·y
' …
Next M%
MatInv = Application.WorksheetFunction.MInverse(A) ' tabla inversa A⁻¹
SolutionMatrix = Application.WorksheetFunction.MMult(MatInv, Z) ' A⁻¹ · ZHay una diferencia con nuestra tabla: en las columnas, O'Neill no pone r, r², r³ y r⁴, sino sus polinomios ortogonales: (r − 1.054769), (r − 1.005327)·(r − 1.128236), etc. Por eso lo que sale de MMult no son los a, sino directamente los λ. Las filas sí son potencias del radio, como en nuestra tabla. Es la receta de Bevington (ecs. 7.30–7.33), que desmontaremos en la Lección 6. Y como en nuestra α, todos los elementos de A contienen solo radios: la tabla es la misma para todas las muestras.
Compruébalo en Excel
Opción rápida, con ESTIMACION.LINEAL. Pon las 13 y en B2:B14, los radios en C2:C14 y, en D, E y F, sus potencias (=C2^2, =C2^3, =C2^4). Entonces:
=ESTIMACION.LINEAL(B2:B14;C2:F14) → a₄ a₃ a₂ a₁ a₀ (¡en orden inverso!)En Excel 365 el resultado ocupa solo cinco celdas en fila; en versiones antiguas, selecciona cinco celdas, escribe la fórmula y pulsa Ctrl+Mayús+Intro.
Opción «como la macro». Monta la tabla α con celdas =SUMAPRODUCTO(C2:C14^n) (con n = j + k, de 0 a 8), el vector β con =SUMAPRODUCTO(B2:B14*C2:C14^j) y resuelve con =MMULT(MINVERSA(α);β). Tendrás los mismos coeficientes, salvo en los últimos decimales: justo la delicadeza de la que hablábamos.
Practica
1. ¿Cuántas curvas (cambios de dirección) puede tener, como máximo, un polinomio de cuatro términos (hasta r³)?
2. La tabla α de VG 304, ¿es distinta de la de VG 189?
3. Con 13 REE y un polinomio de 5 términos, ¿cuántos grados de libertad quedan?
Resumen
- Un polinomio es una línea que puede doblarse: cada término añade una curva posible.
- El criterio sigue siendo mínimos cuadrados, con una condición de equilibrio por término: Σ rk·residuo = 0.
- Esas condiciones forman un sistema de ecuaciones que se escribe como una tabla: α·a = β.
- α solo contiene sumas de potencias de los radios (es igual para todas las muestras); β mezcla radios y datos.
- Se resuelve con la tabla inversa, a = α−1·β: MINVERSA y MMULT en Excel, igual que la macro.
- Con tres términos desaparece la joroba de VG 189; con cinco, los residuos quedan al nivel del error analítico.
- Los coeficientes a0…a4 salen gigantes y se cancelan entre sí: uno a uno no significan nada.
Dónde ha aparecido el radio en esta lección: en la tabla α (las sumas Σr0 … Σr8, solo radios: ④), en el vector β (Σrk·y) y al evaluar la curva en el radio de cada REE. En la macro, además, dentro de los polinomios de O'Neill (las constantes ②).
Próxima lección · Por qué los coeficientes «normales» no sirven para comparar. Veremos por qué a0 salta de 3.19 a −2367 al añadir términos, qué significa que dos coeficientes estén correlacionados al 99.99 % y la segunda mitad de la Fig. 1a de O'Neill (los patrones D y D†).
Fuentes de esta lección
- Bevington & Robinson (2003), §7.1 (ecs. 7.1–7.9: el polinomio y sus ecuaciones normales), §7.2 (ecs. 7.12–7.20: la forma de tabla y la solución con la inversa; Tablas 7.2 y 7.4, el cálculo en una hoja de cálculo) y p. 121 (grados de libertad, ν = N − m).
- O'Neill (2016), p. 1467 (el ajuste puede hacerse con LINEST), pp. 1470–1471 (varianza explicada por cada término) y Tabla 2 (significado de cada término).
- Macro
REE_OrthogRegress_5term: construcción de A y Z y uso de MInverse y MMult.
Lección 4 · de 7 (+ 2 opcionales)
Por qué los coeficientes «normales» no sirven para comparar
En la Lección 3 el ajuste funcionó: la cuártica sigue a VG 189 con residuos del 2 %. Pero sus coeficientes (a0 = −2367, a1 = 9059…) no sirven para describir ni para comparar patrones. Aquí vemos por qué: un síntoma, su causa y su consecuencia. Así entenderemos qué tienen de especial los λ.
- Idea matemática nueva
- Correlación entre coeficientes
- Dónde entra el radio
- El origen r = 0 está a 1.05 Å de los datos; entre 0.98 y 1.16 Å, r, r², r³ y r⁴ casi coinciden
- Muestra guía
- VG 189, y los patrones D y D† de la Fig. 1a
El síntoma: los a cambian al añadir términos
Ajusta VG 189 con 1, 2, 3, 4 y 5 términos y apunta los coeficientes. Los a no se parecen de una fila a otra: a0 pasa de 3.19 a 4.93, −52.5, 221 y −2367. Los λ de O'Neill, en cambio, no se mueven: cada término nuevo añade su λ y deja los anteriores como estaban.
| Términos | a0 | a1 | a2 | a3 | a4 |
|---|---|---|---|---|---|
| 1 | 3.187 | — | — | — | — |
| 2 | 4.934 | −1.656 | — | — | — |
| 3 | −52.46 | 106.30 | −50.60 | — | — |
| 4 | 221.30 | −666.93 | 676.15 | −227.29 | — |
| 5 | −2366.7 | 9059.2 | −13014.3 | 8326.8 | −2001.8 |
| Términos | λ0 | λ1 | λ2 | λ3 | λ4 |
|---|---|---|---|---|---|
| 1 | 3.187 | — | — | — | — |
| 2 | 3.187 | −1.656 | — | — | — |
| 3 | 3.187 | −1.656 | −50.60 | — | — |
| 4 | 3.187 | −1.656 | −50.60 | −227.29 | — |
| 5 | 3.187 | −1.656 | −50.60 | −227.29 | −2001.8 |
Con los a no se pueden comparar dos trabajos que usen distinto número de términos: el a0 de una recta y el a0 de una cuártica son cosas distintas. O'Neill lo resume en su Apéndice 1: si se añade un término al polinomio, todos los coeficientes cambian, y eso impide comparar útilmente los coeficientes de un patrón con los de otro.
Una curiosidad que aprovecharemos en la Lección 6: el último coeficiente de cada fila (recuadrado) coincide con el λ del mismo orden.
La causa: los a describen la curva en r = 0
Un polinomio y = a0 + a1r + a2r2 + … está escrito «desde el origen». a0 es el valor de la curva en r = 0; a1, su pendiente en r = 0; a2 tiene que ver con su curvatura en r = 0, y así sucesivamente (Bevington & Robinson, §7.3). Pero en r = 0 no hay datos. El ion más pequeño, el Lu, mide 0.977 Å, y todo el patrón ocupa solo 0.18 Å. El origen está a 1.05 Å del centro de los datos: casi seis veces la anchura del patrón.
Para saber cuánto vale la curva allí hay que prolongarla muchísimo, y cada polinomio la prolonga a su manera. Aleja el zoom:
Donde hay datos, las curvas de 3, 4 y 5 términos casi coinciden. Fuera, se separan enseguida: a 0.9 Å, solo 0.08 Å por debajo del Lu, ya discrepan en más de un factor 3 en concentración. En r = 0 valen −52.5, 221 y −2367, y esos números son los a0. Describen una zona sin datos, así que dependen de cómo se prolongue la curva y no de la muestra. Bevington lo explica igual: los coeficientes dependen del grado del polinomio porque el sistema de coordenadas se ha elegido sin tener en cuenta dónde están los datos.
Además, las piezas se parecen demasiado
Entre 0.977 y 1.160 Å, las piezas r, r2, r3 y r4 son casi la misma curva: una línea que sube. Estiradas a la misma escala, apenas se distinguen:
Para distinguir una pieza de otra, el ajuste tiene que usar coeficientes enormes y de signo contrario que casi se anulan. En r̄, los cinco términos de la cuártica valen −2367, +9555, −14479, +9771 y −2478, y suman 3.35. Es la «delicadeza» que vimos en la hoja Excel.
La consecuencia: coeficientes correlacionados
¿Qué pasa si repetimos el análisis? Cada medida lleva un error de un 2 %, así que otro análisis de VG 189 daría valores algo distintos, y otros coeficientes. La figura lo simula: añade a cada y un error al azar del 2 % y reajusta, 300 veces.
Coeficientes normales: a0 frente a a1
cada punto es un análisis simulado
Coeficientes de O'Neill: λ0 frente a λ1
los mismos análisis simulados
Mira las nubes. Con los a, los puntos forman una línea finísima: cuando a0 sale alto, a1 sale bajo, siempre en la misma proporción. Eso es una correlación de −0.998 con 2 términos, y de −0.99994 con 5: a0 y a1 no son dos datos independientes, sino casi el mismo dato contado dos veces. Con los λ, la nube es redonda y la correlación es ≈ 0: cada λ lleva su propia información.
Mira también las escalas. Con 5 términos y un 2 % de ruido, a0 varía con una desviación típica de unas 1100 unidades (puede salir cualquier cosa entre unos −6000 y +1000), aunque la curva apenas se mueve: entre un 0.9 y un 1.9 % en cada REE. λ0, en cambio, solo varía ±0.0055. Es lógico: λ0 es la media de 13 valores, y el error de una media de 13 valores es el error de cada uno dividido por √13: 0.02 / √13 = 0.0055.
Esto tiene una consecuencia práctica seria. Si alguien publica los a con sus errores por separado (a0 = −2367 ± 1100, a1 = 9059 ± 4100…), parece que no sabemos casi nada, cuando la curva está determinada al 2 %. Los errores de coeficientes correlacionados solo sirven acompañados de su tabla completa de covarianzas (Bevington, §7.2, ecs. 7.23–7.27). Bevington lo muestra con un ejemplo en la p. 126: con los errores individuales obtiene V = 2.45 ± 0.14; teniendo en cuenta las correlaciones, 2.45 ± 0.02.
Otro detalle en el que entra el radio: esa tabla de covarianzas es la inversa de α (ec. 7.25), y α solo contiene radios. Las correlaciones entre los a son, por tanto, las mismas para todas las muestras.
La segunda mitad de la Fig. 1a: D y D†
O'Neill lo resume con dos patrones inventados. D es una recta, y = 2 + 3r. D† es casi igual, pero con una ligera curvatura; ajustado con tres términos da y = 24.7 − 39.7r + 20r2. Mirando los coeficientes (2, 3 y 0 frente a 24.7, −39.7 y 20), nadie diría que se parecen. En la forma de O'Neill: D = 5.16 + 3f1 + 0f2 y D† = 5.16 + 3f1 + 20f2. Mismo nivel medio, misma pendiente: solo cambia la curvatura, y el número que la mide es exactamente ese 20.
| Patrón | a0 | a1 | a2 | λ0 | λ1 | λ2 |
|---|
Cambiar una sola propiedad del patrón (su curvatura) no debería cambiar los números que describen las otras (su nivel y su pendiente). Con los a las cambia; con los λ, no.
Lo que querríamos
Unos coeficientes que:
- no cambien al añadir o quitar términos;
- no estén correlacionados entre sí;
- tengan cada uno un significado claro: nivel medio, pendiente, curvatura…
En la Lección 2 ya lo conseguimos para la recta: bastó medir el radio desde su media (r − 1.05477), y la correlación entre ordenada y pendiente pasó de −0.998 a 0. Los polinomios ortogonales hacen lo mismo con todos los términos: cada pieza se «centra» respecto a las anteriores. Es la Lección 5.
Practica
1. Un colega publica a0 = −2100 para su basalto (ajuste de 5 términos) y tú obtienes a0 = −2400 para el tuyo. ¿Qué puedes concluir sobre sus REE?
2. Si el error analítico fuese del 5 % en vez del 2 %, ¿cuál sería la desviación típica de λ0 (con 13 REE)?
3. D y D† tienen el mismo λ0 y el mismo λ1. ¿En qué se diferencian?
Resumen
- Los coeficientes a cambian al añadir términos: a0 pasa de 3.19 a −2367 en VG 189. Los λ no cambian.
- La causa: los a describen la curva en r = 0, a 1.05 Å de los datos, adonde hay que prolongarla muchísimo. Además, las piezas r, r², r³ y r⁴ son casi idénticas entre 0.98 y 1.16 Å.
- La consecuencia: los a están correlacionados (−0.998 con 2 términos, −0.99994 con 5) y sus errores individuales son enormes aunque la curva esté determinada al 2 %.
- D y D†: con los a parecen patrones distintos; con los λ se ve que solo difieren en la curvatura (λ2 = 20).
- Queremos coeficientes estables, independientes y con significado: eso son los λ.
Dónde ha aparecido el radio en esta lección: en la distancia del origen a los datos (r̄ = 1.05 Å, casi seis veces la anchura del patrón), que es el «brazo de palanca» de los a; en que todos los radios se parecen (0.98–1.16 Å), lo que hace casi idénticas las piezas r…r⁴; y en la tabla de covarianzas, la inversa de α, que solo depende de los radios.
Próxima lección · Los polinomios ortogonales: fabricarlos con los radios. Veremos qué significa «ortogonal» (piezas que no se pisan) y cómo salen, solo de los 13 radios, las constantes 1.05477, 1.00533, 1.12824… de la Tabla 2, paso a paso con la receta de Bevington (§7.3).
Fuentes de esta lección
- O'Neill (2016), p. 1465 (los coeficientes del polinomio están muy correlacionados), Fig. 1a y su pie (patrones D y D†, p. 1466), y Apéndice 1 (al añadir un término cambian todos los coeficientes).
- Bevington & Robinson (2003), §7.2 (ecs. 7.23–7.27: la tabla de errores y covarianzas es la inversa de α; ejemplo de la p. 126 y Tabla 7.3) y §7.3, p. 127 (por qué los coeficientes dependen del grado del polinomio).
- Simulación de análisis repetidos: 300 ajustes de VG 189 con un error aleatorio normal del 2 % en cada ln([REE]/[REE]CI). Desviaciones típicas y correlaciones calculadas también de forma exacta a partir de la inversa de α.
Lección 5 · de 7 (+ 2 opcionales)
Los polinomios ortogonales: fabricarlos con los radios
En la Lección 4 pedimos coeficientes que no cambien al añadir términos y que no estén correlacionados. La solución de O'Neill, tomada de Bevington, es cambiar las piezas del polinomio: en vez de 1, r, r², r³ y r⁴, usa cinco piezas «ortogonales» fabricadas a medida para los 13 radios. Aquí veremos qué significa «ortogonal» y cómo salen, paso a paso, las constantes de su Tabla 2.
- Idea matemática nueva
- Ortogonalidad: dos listas de números cuya suma de productos es 0
- Dónde entra el radio
- ② las constantes β, γ, δ, ε y ③ la tabla de valores fn(ri)
- Muestra guía
- VG 189, solo para comprobar que los λ anteriores no se mueven
La idea: caminar en direcciones perpendiculares
Imagina los 13 valores y de una muestra como un punto, y cada pieza del polinomio como una dirección en la que se puede caminar. Ajustar es decidir cuánto caminar en cada dirección para llegar a ese punto, o lo más cerca posible. Esas distancias son los coeficientes.
Con 13 valores harían falta 13 dimensiones, que no se pueden dibujar. Pero con dos piezas todo cabe en un plano. Prueba a cerrar el ángulo entre las dos direcciones:
Esto no es solo una analogía. Las piezas 1 y r, evaluadas en los 13 radios, forman exactamente un ángulo de 3.27° (enseguida veremos cómo se calcula). Es la recta sin centrar de las Lecciones 2 y 4, y los números encajan: la correlación de −0.998 entre a0 y a1 es −cos(3.27°), y el error de a0 (0.097) era 17.5 veces el de λ0 (0.0055), justo 1/sen(3.27°). Con r, r2, r3 y r4 es aún peor: son cinco direcciones casi paralelas.
Con direcciones perpendiculares, en cambio, cada coeficiente dice cuánto caminar en su dirección, sin rodeos, y corregir uno no obliga a corregir los demás. Eso significa ortogonal: perpendicular. La recta centrada de la Lección 2 ya lo era (1 y r − r̄ forman 90°), y por eso allí la correlación entre ordenada y pendiente era 0.
La prueba de la perpendicularidad: multiplicar y sumar
En el papel, dos direcciones son perpendiculares cuando, al multiplicar sus componentes una a una y sumar, sale 0. El este es (1, 0) y el norte, (0, 1): 1 × 0 + 0 × 1 = 0. Con el este y el nordeste, (1, 0) y (1, 1), sale 1 × 1 + 0 × 1 = 1: no son perpendiculares.
Con las piezas del polinomio se hace igual, pero con 13 números: el valor de cada pieza en cada REE. Dos piezas fj y fk son ortogonales si (Bevington & Robinson, nota 1 de la p. 128):
Σ fj(ri) · fk(ri) = 0 (sumando las 13 REE, para j ≠ k)
Es decir: se multiplica el valor de las dos piezas en el La, luego en el Ce, y así hasta el Lu, y se suman los 13 productos. Si los positivos compensan exactamente a los negativos, las piezas no se pisan. Compruébalo:
Con cualquier pareja de piezas de O'Neill, las barras positivas y negativas se compensan y el ángulo es 90.0°. Con las piezas simples no hay compensación posible: todos los radios son positivos, así que todos los productos también lo son. Los ángulos entre ellas van de 3.3° (piezas vecinas, como r y r²) a 13.1° (1 y r⁴): las cinco caben en un cono de 13°.
Para curiosos: el ángulo se calcula como en el papel, coseno = suma de productos / (longitud × longitud), donde la longitud de una pieza es la raíz de la suma de sus cuadrados. Para 1 y r: 13.712 / (√13 × √14.5102) = 13.712 / (3.6056 × 3.8092) = 0.9984, que corresponde a 3.27°.
La receta: cada pieza nueva, perpendicular a todas las anteriores
La receta de Bevington & Robinson (§7.3) fabrica las piezas una a una. Cada pieza sube un grado y trae tantas constantes libres como piezas hay ya fabricadas, y cada una de esas piezas le impone una prueba:
| Pieza | Forma | Constantes | Debe ser ortogonal a | Nombre (Tabla 2) |
|---|---|---|---|---|
| f0 | 1 | — | — | Nivel medio |
| f1 | (r − β) | β | f0 | Pendiente |
| f2 | (r − γ1)(r − γ2) | γ1, γ2 | f0, f1 | Curvatura |
| f3 | (r − δ1)(r − δ2)(r − δ3) | δ1, δ2, δ3 | f0, f1, f2 | Inflexiones en los extremos (S) |
| f4 | (r − ε1)···(r − ε4) | ε1…ε4 | f0…f3 | Forma de W |
Siempre hay tantas incógnitas como pruebas, así que la receta nunca se atasca. Cada constante es un radio en el que la pieza vale cero.
¿Y por qué estas pruebas? Bevington lo razona así: al añadir una pieza, los coeficientes anteriores no deben cambiar. Si f1 no sumara 0 en las 13 REE, al añadirla subiría o bajaría el nivel medio, y λ0 tendría que corregirse. Superar la prueba garantiza que la pieza nueva no «pisa» a las anteriores. La figura lo comprueba con VG 189:
Paso 1 (ec. 7.35). La única β que supera la prueba es la media de los 13 radios: Σr / 13 = 13.712 / 13 = 1.054769 Å. Con cualquier otra, f1 no sumaría 0 y, al añadirla, el nivel medio λ0 se movería. Es la recta centrada de la Lección 2.
Paso 2 (ec. 7.36). Las dos pruebas solo se cumplen con γ1 = 1.00533 y γ2 = 1.12824 Å. Bevington escribe la segunda prueba con r en lugar de f1; como f2 ya suma 0, es lo mismo. La pieza f2 es una U: positiva en los extremos (La, Ce y Er–Lu) y negativa en el centro (Pr–Ho). Fíjate en que λ2 = −50.60 sale igual con cualquier γ: es el coeficiente del término más alto, que coincide con a2 (Lección 4). Lo que deciden las γ es si λ0 y λ1 se quedan quietos.
Pasos 3 y 4 (ec. 7.37). Igual, con tres y cuatro constantes: δ = 0.99141, 1.06055 y 1.14552; ε = 0.98482, 1.03052, 1.10441 y 1.15343 Å. Bevington lo despacha con una frase: la extensión a órdenes mayores es directa. f3 tiene forma de S y f4, de W. Son las constantes de la Tabla 2 de O'Neill.
Para curiosos: cómo se calculan las γ sin tanteo
Desarrolla f2 = r2 + b·r + c, donde b = −(γ1 + γ2) y c = γ1·γ2. Las dos pruebas se convierten en dos ecuaciones corrientes, con las sumas de potencias de los radios de la Lección 3:
Σr2 + b·Σr + 13·c = 0 → 14.510216 + 13.712·b + 13·c = 0
Σr3 + b·Σr2 + c·Σr = 0 → 15.405677 + 14.510216·b + 13.712·c = 0
Resolviéndolas, b = −2.133563 y c = 1.134247. Las γ son los dos valores de r en los que f2 vale 0, que da la fórmula de la ecuación de segundo grado: γ = (2.133563 ∓ √(2.133563² − 4 × 1.134247)) / 2 = (2.133563 ∓ 0.122909) / 2 = 1.005327 y 1.128236.
Con f3 y f4 se hace lo mismo con 3 y 4 ecuaciones. Desarrolladas, f3 = r³ − 3.19748 r² + 3.402002 r − 1.204446 y f4 = r⁴ − 4.27318 r³ + 6.839052 r² − 4.85869 r + 1.292811. Una curiosidad: entre cada dos ceros seguidos de una pieza cae exactamente un cero de la pieza anterior (β entre las dos γ, cada γ entre dos δ, cada δ entre dos ε).
Solo hacen falta los radios
Repasa la receta: en ningún momento ha hecho falta una concentración. Las pruebas solo usan los valores de las piezas en los 13 radios, así que las constantes dependen solo de los radios (y de qué REE se usan). Bevington lo advierte expresamente (p. 129): estos parámetros adicionales son funciones solo de la variable independiente, que aquí es r.
Por eso O'Neill las calcula una sola vez, las publica en su Tabla 2 y sirven para cualquier basalto. Su macro ni siquiera las calcula: las lleva escritas, con hasta nueve decimales. Esta es la línea que reconstruye la curva ajustada, partida en cuatro para que quepa:
' RR(M%) = radio de la REE número M; Q(1)…Q(5) = λ0…λ4
XCalc(M%) = Q(1) + Q(2) * (RR(M%) - 1.054769) _
+ Q(3) * (RR(M%) - 1.005327429) * (RR(M%) - 1.128236038) _
+ Q(4) * (RR(M%) - 1.060548105) * (RR(M%) - 1.145519887) * (RR(M%) - 0.991412204) _
+ Q(5) * (RR(M%) - 1.104414973) * (RR(M%) - 1.153426708) * (RR(M%) - 0.984820219) * (RR(M%) - 1.030518142)Es la fórmula de O'Neill escrita en Visual Basic, λ0 + λ1·f1 + λ2·f2 + λ3·f3 + λ4·f4, con las constantes de la Tabla 2 (el orden de los factores da igual).
Un detalle: las fórmulas de Bevington llevan un peso 1/σ² en cada punto. O'Neill da el mismo peso a las 13 REE, y así los pesos no influyen en las constantes. Si cada REE llevara su propio peso, las constantes dependerían también de los pesos.
La tabla de valores de cada pieza
Con las constantes fijadas, basta sustituir cada radio en cada pieza. Por ejemplo, para el La (r = 1.160 Å):
f1(La) = 1.160 − 1.054769 = 0.1052
f2(La) = (1.160 − 1.005327) × (1.160 − 1.128236) = 0.154673 × 0.031764 = 0.004913
Así sale toda la tabla, que es la de las columnas f de la Tabla 1 de O'Neill. Los valores de esta página se han calculado a partir de los radios y coinciden con los suyos cifra a cifra:
| REE | r (Å) | f1 | f2 (×10⁻³) | f3 (×10⁻³) | f4 (×10⁻⁶) |
|---|
Dos cosas saltan a la vista. La fila de sumas vale 0: cada pieza es ortogonal a f0 = 1. Y los valores son cada vez más pequeños (por eso O'Neill usa ×10⁻³ y ×10⁻⁶), porque cada pieza es un producto de diferencias de radios de unas centésimas de Å. Por eso los λ de orden alto son números grandes (λ4 = −2002 en VG 189): multiplican a una pieza diminuta. Un λ grande no significa una forma exagerada.
La fila del Eu se calcula igual (r = 1.066 Å), aunque el Eu no entra en el ajuste: servirá para calcular el Eu* en la Lección 7. Y otro dato que usaremos en las Lecciones 6 y 7: la suma de los cuadrados de cada columna (13 REE) vale 13; 0.04722; 1.029 × 10⁻⁴; 2.550 × 10⁻⁷ y 5.567 × 10⁻¹⁰.
¿Y si falta algún elemento?
Las constantes dependen de qué radios entran en la receta. Si se incluyera el Eu, o se quitara el Ce, saldrían otras. Pruébalo:
| Constante | Tabla 2 (13 REE) | Con tu selección | Diferencia |
|---|
¿Hay que recalcular las constantes cuando a una muestra le falta un elemento? O'Neill dice que no (p. 1466). Las covarianzas entre los λ solo son exactamente cero con las 13 REE, pero si falta alguna, en la práctica importa poco. Lo que de verdad importa es usar los mismos polinomios para todos los patrones, para que sus λ sean comparables. Su macro hace eso: si falta un dato, lo salta y sigue usando las constantes de la Tabla 2. El ajuste sigue siendo correcto (veremos por qué en la Lección 6); solo que los λ dejan de ser del todo independientes. Recalcular las constantes para otro conjunto de REE es posible, pero entonces los λ ya no se pueden comparar con los de O'Neill. Y con sentido común: si faltan todas las REE ligeras, la forma del patrón no queda definida. En la Lección 8 lo probaremos con VG 189 sin el Ce.
Practica
1. ¿Son ortogonales estas dos listas de tres números, A = (1, 2, −1) y B = (2, 0, 2)?
2. Calcula f2 para el Nd (r = 1.109 Å), con γ1 = 1.005327 y γ2 = 1.128236.
3. Si alguien incluyera el Eu (r = 1.066 Å) en la receta, con 14 REE, ¿cuánto valdría β? Recuerda que los 13 radios de O'Neill suman 13.712 Å.
4. Pasas de ajustar VG 189 a ajustar VG 4080, un basalto enriquecido. ¿Qué cambia?
Resumen
- Ortogonal significa perpendicular: dos piezas lo son si la suma de sus productos en las 13 REE es 0 (los productos positivos y negativos se compensan).
- Con piezas perpendiculares, cada coeficiente se lee sin rodeos. Las piezas simples 1, r, r², r³ y r⁴ son casi paralelas (de 3° a 13° entre sí), y por eso sus coeficientes se disparan y se correlacionan.
- La receta (Bevington, §7.3): cada pieza nueva tiene tantas constantes como piezas anteriores y debe ser ortogonal a todas. Salen β = 1.05477 (la media de los radios), γ = 1.00533 y 1.12824, las tres δ y las cuatro ε de la Tabla 2.
- Superar las pruebas equivale a que añadir una pieza no cambie los λ anteriores.
- Las constantes solo dependen de los radios (y de qué REE se usan): se calculan una vez y valen para todas las muestras. O'Neill recomienda no cambiarlas aunque falte algún elemento.
- Sustituyendo cada radio se obtiene la tabla fn(ri) (Tabla 1): cada columna suma 0 y los valores son cada vez más pequeños.
Dónde ha aparecido el radio en esta lección: en todo. ② Las constantes β, γ, δ y ε salen solo de los 13 radios. ③ La tabla fn(ri) es cada radio sustituido en cada pieza. Las concentraciones solo han aparecido para comprobar, con VG 189, que los λ anteriores no se mueven.
Próxima lección · Calcular los λ: tres caminos, un resultado. Con las piezas ya fabricadas, cada λ se calcula por su cuenta con una simple división: λn = Σ y·fn / Σ fn² (Bevington, ec. 7.38). Lo haremos con VG 189 y veremos que la macro llega al mismo resultado resolviendo un sistema de 5 ecuaciones, y que ESTIMACION.LINEAL también sirve.
Fuentes de esta lección
- Bevington & Robinson (2003), §7.3, pp. 127–130: ec. 7.28 (el polinomio de piezas ortogonales), nota 1 de la p. 128 (definición: Σ fj·fk = 0 para j ≠ k), ecs. 7.34–7.37 (las pruebas que fijan β, γ y δ, y la advertencia de que solo dependen de la variable independiente) y ec. 7.38 (adelanto de la Lección 6).
- O'Neill (2016), Tabla 1 (valores de f1…f4 en cada REE, reproducidos aquí), Tabla 2 (constantes y significado de cada término) y p. 1466 (la ortogonalidad es exacta solo con las 13 REE; hay que usar los mismos polinomios para todos los patrones).
- Macro REE_OrthogRegress_5term de la hoja de cálculo de O'Neill: las constantes van escritas en el código (tabla del sistema y línea XCalc).
- Cálculos propios: constantes recalculadas a partir de los 13 radios (coinciden con la Tabla 2 y con la macro), ángulos entre piezas y constantes para otros conjuntos de REE.
Lección 6 · de 7 (+ 2 opcionales)
Calcular los λ: tres caminos, un resultado
Con las piezas fabricadas en la Lección 5 falta el último paso: averiguar cuánto hay de cada pieza en una muestra, es decir, sus λ. Hay tres caminos: la fórmula directa de Bevington, el sistema que resuelve la macro de O'Neill y un ajuste corriente con ESTIMACION.LINEAL. Los tres dan lo mismo para VG 189: 3.187, −1.656, −50.60, −227.29 y −2001.8. Al final veremos qué significa cada λ y cómo se lee el diagrama λ2–λ1.
- Idea matemática nueva
- Proyección: la «sombra» de los datos sobre cada pieza
- Dónde entra el radio
- ④ la tabla del sistema de la macro, solo con radios, y los fn(ri) de la Lección 5 en cada fórmula
- Muestras
- VG 189, las otras tres muestras de ejemplo y los 616 vidrios OFB de la hoja de O'Neill
Camino 1: la sombra de los datos sobre cada pieza
En la Lección 5 imaginamos los 13 valores y de una muestra como un punto, y cada pieza como una dirección. Cuando las direcciones son perpendiculares hay un atajo para saber cuánto caminar en cada una: medir la sombra del punto sobre esa dirección, con la luz cayendo en perpendicular. Compruébalo:
La sombra solo depende del punto y de su dirección, no de las demás. La cuenta es la prueba del producto de la Lección 5 (multiplicar y sumar), dividida por la misma prueba de la dirección consigo misma. Con las piezas de O'Neill es igual, pero con 13 números en vez de 2 (Bevington & Robinson, ec. 7.38):
λn = Σ yi·fn(ri) / Σ fn(ri)²
Arriba se multiplica cada dato por el valor de la pieza en su REE y se suma; abajo, la suma de cuadrados de la pieza, que ya vimos en la Lección 5. Para λ0 la pieza es f0 = 1 y la fórmula da la media de las 13 y, como ya sabíamos. Pruébalo con cada muestra:
| Coeficiente | Σ y·fn | Σ fn² | λn |
|---|
De esta fórmula salen dos propiedades muy útiles:
- Cada λ se calcula por su cuenta. λ2 no necesita a λ1 ni a λ3, así que añadir o quitar términos no cambia los demás. O'Neill lo destaca (p. 1467): así se pueden comparar los λ de estudios poco precisos, que usan menos términos, con los de los análisis más precisos.
- Solo λ0 depende del nivel. Como f1…f4 suman 0 (Lección 5), sumar la misma cantidad a las 13 y no cambia Σ y·fn. Y multiplicar todas las concentraciones por un factor suma su logaritmo a todas las y. Si un magma se concentra o se diluye por igual en todas las REE, solo cambia λ0. Por eso, dice O'Neill (p. 1466), un error analítico sistemático que afecte igual a todas las REE solo pasa a λ0, y no a los coeficientes de forma.
Camino 2: el sistema de la macro
La macro de O'Neill no usa la fórmula directa. Hace lo mismo que vimos en la Lección 3: recorre las REE sumando una tabla A de 5 × 5 y un vector Z, y resuelve A·λ = Z con MINVERSA y MMULT (líneas resumidas):
For M% = 1 To NoData% ' recorre las REE (salta vacíos, ceros y el Eu)
A(1, 1) = A(1, 1) + 1
A(1, 2) = A(1, 2) + (RR(M%) - 1.054769) ' Σ f1
A(2, 2) = A(2, 2) + RR(M%) * (RR(M%) - 1.054769) ' Σ r·f1
' … así hasta 25 sumas: A(i, j) = Σ r^(i−1)·f(j−1)(r) ← solo radios
Z(1, 1) = Z(1, 1) + X(M%) ' Σ y
Z(2, 1) = Z(2, 1) + X(M%) * RR(M%) ' Σ r·y … hasta Σ r⁴·y
Next M%
MatInv = Application.WorksheetFunction.MInverse(A)
SolutionMatrix = Application.WorksheetFunction.MMult(MatInv, Z) ' = λ0 … λ4La tabla A solo contiene radios: es la estación ④ del mapa. El vector Z es exactamente el β de la Lección 3 (Σ rk·y), y las cinco ecuaciones son las 7.30–7.33 de Bevington. Así queda para VG 189:
| Camino | λ0 | λ1 | λ2 | λ3 | λ4 |
|---|
Con las 13 REE, la tabla es una escalera: por encima de la diagonal todo vale 0, porque cada pieza es ortogonal a las potencias más bajas de r. La primera ecuación, 13·λ0 = 41.437, da λ0 = 3.187 sin mirar las demás. La segunda, 13.712·λ0 + 0.04722·λ1 = 43.629, da λ1 = −1.656 en cuanto se conoce λ0. Y así hasta λ4: cada ecuación aporta el coeficiente siguiente, como explica Bevington (p. 130). La diagonal son las Σ fn² de la Lección 5.
¿Para qué un sistema, si la fórmula directa es tan sencilla? Porque el sistema sigue valiendo cuando falta un dato. Quita el Ce en la figura: la escalera se pierde, pero el sistema da el ajuste por mínimos cuadrados correcto con las mismas piezas de la Tabla 2, con λ muy parecidos a los de las 13 REE, como anunciaba O'Neill. La fórmula directa, en cambio, da disparates (λ4 = +40 062). Con 12 REE las piezas ya no son perpendiculares, y la sombra deja de ser el coeficiente, como en la primera figura con un ángulo distinto de 90°. Esta es la explicación que prometimos en la Lección 5: la macro puede saltarse un dato sin cambiar las constantes.
Para curiosos: dos detalles numéricos de la macro
El sistema tampoco necesita constantes perfectas. La macro redondea β a 1.054769 (el valor exacto es 1.0547692…), y con eso Σ f1 vale 0.000003 en vez de 0. La fórmula directa lo notaría en la cuarta cifra de λ1 (−1.6556 en vez de −1.6558); el sistema, no.
La macro guarda los radios en precisión simple (Dim RR(14) As Single), con unas siete cifras. Por eso su λ4 de VG 189 es −2001.816 y no −2001.807. La diferencia es irrelevante: la incertidumbre de λ4 con un 2 % de error analítico es ±848 (Lección 4).
Camino 3: con un ajuste corriente
a) ESTIMACION.LINEAL con las piezas como columnas. Es la vía más cómoda en Excel, y la propone el propio O'Neill (p. 1467): hacer una regresión múltiple usando los polinomios ortogonales como variables «independientes». Pon las 13 y (sin Eu) en B2:B14 y, en C2:F14, los valores de f1…f4 de cada REE, calculados con fórmulas a partir del radio (en la columna A):
C2: =A2-1.054769 f1
D2: =(A2-1.005327429)*(A2-1.128236038) f2
E2: =(A2-0.991412204)*(A2-1.060548105)*(A2-1.145519887) f3
F2: =(A2-0.984820219)*(A2-1.030518142)*(A2-1.104414973)*(A2-1.153426708) f4
=ESTIMACION.LINEAL(B2:B14;C2:F14) → λ₄ λ₃ λ₂ λ₁ λ₀ (en orden inverso, como en la Lección 3)Para VG 189 devuelve −2001.8, −227.29, −50.60, −1.656 y 3.187. Si usas solo C2:D14 obtienes λ2, λ1 y λ0, con los mismos valores: con piezas ortogonales, quitar términos no cambia los demás. Con las columnas r, r², r³ y r⁴, en cambio, los coeficientes cambiaban (Lección 4).
b) Reordenando los coeficientes de un ajuste normal. Si ya tienes los a de la Lección 3, O'Neill indica que basta reordenarlos (p. 1467). Se hace de arriba abajo, usando los polinomios desarrollados de la Lección 5 con más cifras:
| Fórmula | VG 189 |
|---|---|
| λ4 = a4 | −2001.8067 |
| λ3 = a3 + 4.27318004·λ4 | 8326.7924 − 8554.0803 = −227.288 |
| λ2 = a2 + 3.19748020·λ3 − 6.83905178·λ4 | −13014.3108 − 726.7487 + 13690.4595 = −50.600 |
| λ1 = a1 + 2.13356347·λ2 − 3.40200168·λ3 + 4.85868971·λ4 | 9059.2260 − 107.9583 + 773.2340 − 9726.1575 = −1.656 |
| λ0 = a0 + 1.054769·λ1 − 1.13424664·λ2 + 1.20444581·λ3 − 1.29281055·λ4 | −2366.65963 − 1.74653 + 57.39290 − 273.75603 + 2587.95678 = 3.1875 |
Funciona, pero fíjate en las cancelaciones: para obtener λ0 = 3.19 se suman números de miles. Con las constantes redondeadas a 4 decimales saldría λ0 = 3.206 y λ1 = −1.619, errores mayores que la incertidumbre analítica (0.0055 y 0.09). Si eliges este camino, usa todas las cifras. El camino a) no tiene ese problema.
Qué significa cada λ
Cada λ multiplica a una pieza, y el patrón es la suma de las cinco contribuciones. El tamaño del número no dice cuánto pesa en el patrón, porque las piezas son cada vez más pequeñas (Lección 5). Para VG 189:
| λ | Pieza | Qué mide (Tabla 2) | VG 189 | Qué hace en VG 189 |
|---|---|---|---|---|
| λ0 | 1 | Nivel medio: eλ₀ es la media geométrica de [REE]/[REE]CI | 3.187 | 24 veces el condrito |
| λ1 | f1, recta | Pendiente. Positiva: ligeras enriquecidas (el La, a la derecha, más alto) | −1.656 | La −16 %, Lu +14 % |
| λ2 | f2, U | Curvatura. Negativa: joroba (REE medias por encima de la recta); positiva: valle | −50.60 | Gd +20 %; La −22 %, Lu −20 % |
| λ3 | f3, S | Inflexiones en los extremos del patrón | −227.3 | La −5 %, Lu +5 % |
| λ4 | f4, W | Forma de W | −2002 | unos −2 % en La, Gd y Lu |
λ4 = −2002 cambia el patrón en un 2 %; λ1 = −1.66, en un 16 %. La figura descompone las cuatro muestras de ejemplo en sus contribuciones, como la Fig. 1b de O'Neill, y las sitúa en el diagrama λ2–λ1. Elige Punto libre y arrastra el punto por el diagrama para ver qué forma corresponde a cada zona.
El patrón
Lo que aporta cada término
Diagrama λ2–λ1
Cómo leer el diagrama: hacia la derecha, patrones más enriquecidos en REE ligeras (λ1 mayor); hacia arriba, menos joroba y más valle (λ2 mayor). Los vidrios OFB forman una banda inclinada, λ2 = −22.6 + 5.36·λ1: en ellos, la curvatura de los patrones va ligada a su pendiente (O'Neill, p. 1469). Los basaltos de Hawái quedan aparte, por debajo de la banda: para la misma pendiente, mucha más joroba. Lo que esas posiciones significan en términos de procesos (fusión, cristalización) es el tema de los vectores ψ de la Lección 9.
Separar el nivel de la forma tiene otra ventaja. En los OFB, λ0 sube a medida que el magma evoluciona y pierde MgO (λ0 = 5.100 − 0.287·[MgO]), mientras que λ1 y λ2 apenas cambian. O'Neill concluye que la forma de los patrones se hereda de los magmas parentales (p. 1468). ¿Y λ3 y λ4? Con la precisión analítica habitual suelen ser poco distinguibles de cero; lo veremos en la Lección 7.
Practica
1. Para VG 189, Σ y·f2 = −0.0052060 y Σ f2² = 0.00010288. ¿Cuánto vale λ2?
2. Un magma con el patrón de VG 189 se diluye a la mitad: todas las concentraciones de REE se reducen al 50 %. ¿Qué λ cambian?
3. ¿Por qué la macro resuelve un sistema de 5 ecuaciones en vez de usar la fórmula directa?
4. VG 4080 está arriba a la derecha en el diagrama (λ1 = 11.0, λ2 = 41.9). ¿Qué forma tiene su patrón?
Resumen
- Tres caminos, un resultado: la fórmula directa λn = Σ y·fn / Σ fn², el sistema A·λ = Z de la macro y ESTIMACION.LINEAL (con las piezas como columnas, o reordenando los a).
- La fórmula directa es una sombra: solo da el coeficiente porque las piezas son perpendiculares. Así, cada λ se calcula por su cuenta, y añadir o quitar términos no cambia los demás.
- Solo λ0 depende del nivel: diluir o concentrar un magma, o un error de calibración común a todas las REE, solo cambia λ0.
- El sistema de la macro es una escalera con las 13 REE y sigue dando el ajuste correcto aunque falte un dato; la fórmula directa, no.
- Reordenar los a funciona, pero exige muchas cifras: hay cancelaciones enormes.
- Cada λ mide una forma, y su tamaño no es su efecto (λ4 = −2002 ≈ 2 %). En el diagrama λ2–λ1, los OFB siguen λ2 = −22.6 + 5.36·λ1.
Dónde ha aparecido el radio en esta lección: ④ la tabla A de la macro solo contiene radios (sumas de potencias de r por las piezas), y su diagonal son las Σ fn² de la Lección 5. En la fórmula directa, el radio entra a través de los fn(ri) (③) y de Σ fn². Los datos solo entran en Σ y·fn y en Z.
Próxima lección · ¿Es bueno el ajuste? Residuos, Eu/Eu*, χ² e incertidumbres. Veremos cuánto se aparta cada REE de la curva (el cociente observado/calculado de la macro), cómo se obtiene la anomalía de Eu con la curva evaluada en r = 1.066 Å, qué es el χ² reducido, cómo infiere O'Neill la precisión de los análisis (2.14 %) y de dónde salen los errores de los λ (los factores ζ).
Fuentes de esta lección
- Bevington & Robinson (2003), §7.3, pp. 128–130: ecs. 7.30–7.33 (el sistema), p. 130 (cada ecuación da el coeficiente siguiente) y ec. 7.38 (la fórmula directa).
- O'Neill (2016): p. 1466 (un error sistemático común a todas las REE solo afecta a λ0), p. 1467 (quitar términos no cambia los demás; regresión múltiple con los polinomios ortogonales como variables; reordenar los coeficientes de un ajuste polinómico), p. 1468 (λ0 frente a MgO; la forma se hereda de los magmas parentales), Fig. 2f y p. 1469 (λ2 frente a λ1 en los OFB) y Fig. 1b (descomposición en términos).
- Macro REE_OrthogRegress_5term: construcción de A y Z y solución con MInverse y MMult; radios declarados en precisión simple.
- Datos del diagrama: hoja REEshapes_reference_datasets_for_plots.xls de O'Neill, con los λ calculados por su macro para 616 vidrios OFB (Jenner & O'Neill, 2012) y 119 basaltos del sondeo HSDP2 de Hawái (Feigenson et al., 2003).
Lección 7 · de 7 (+ 2 opcionales)
¿Es bueno el ajuste? Residuos, Eu/Eu*, χ² e incertidumbres
Ya tenemos los λ. Falta comprobar lo bien que describen el patrón y cuánto nos podemos fiar de ellos. La macro de O'Neill calcula tres cosas más: el cociente observado/calculado de cada REE (incluido el Eu, que da la anomalía Eu/Eu*) y un χ² reducido para cada muestra. Con muchas muestras, además, se puede deducir la precisión típica de los análisis. Con esa precisión y los factores ζ salen los errores de los λ.
- Idea matemática nueva
- χ² reducido y propagación de errores
- Dónde entra el radio
- ⑤ la curva evaluada en r = 1.066 Å (el Eu*) y ⑥ los factores ζ = 1/√Σfn², que solo dependen de los radios
- Muestras
- VG 189, las otras tres muestras de ejemplo y los 616 vidrios OFB de la hoja de O'Neill
Cuánto se aparta cada REE: el cociente observado/calculado
Con los λ ya calculados, la macro evalúa la curva en los 14 radios, incluido el del Eu (estación ⑤ del mapa). Después compara cada dato con la curva:
XCalc(M%) = Q(1) + Q(2) * (RR(M%) - 1.054769) + … ' la curva en el radio de la REE M (Lección 5)
DataOut.Cells(Sample%, 8 + M%) = Exp(X(M%) - XCalc(M%)) ' [REE]obs / [REE]calc
' … y pinta de rojo las celdas por debajo de 0.95 o por encima de 1.05Como y es un logaritmo, la diferencia y − ycalc es el logaritmo del cociente, y exp() lo deshace: 1.00 significa que el dato cae sobre la curva; 0.98, que está un 2 % por debajo. O'Neill llama suavidad a lo cerca que están estos cocientes de 1 (p. 1470). Mira cómo cambian al añadir términos:
La anomalía de Eu
El Eu queda fuera del ajuste «por la probabilidad de anomalías debidas al Eu2+» (Apéndice 1). Parte del Eu de un magma está como Eu2+, un ion mucho mayor (1.25 Å en coordinación 8, frente a 1.066 Å del Eu3+), parecido al Sr2+, que la plagioclasa acepta con facilidad. Si la plagioclasa ha cristalizado o se ha acumulado, el Eu se sale del patrón suave de las demás REE.
Por eso la curva, calculada con las otras 13 REE, dice cuánto Eu «debería» haber si no hubiera anomalía: es el Eu*, el valor de la curva en r = 1.066 Å. Y el cociente observado/calculado del Eu es precisamente Eu/Eu*. La hoja de O'Neill lo indica en la cabecera de esa columna: «this is Eu/Eu*». Para VG 189, y(Eu) = 3.2547 y la curva vale 3.3479 en 1.066 Å, así que Eu/Eu* = e3.2547 − 3.3479 = 0.911: le falta un 9 % de Eu.
En basaltos de patrón suave los dos Eu* casi coinciden: para VG 189, 0.911 con la curva y 0.912 con Sm y Gd. La diferencia está en de qué depende cada uno. El clásico depende a partes iguales de dos análisis. El de O'Neill es una media ponderada de todas las REE, con pesos que solo dependen de los radios: Sm 38 %, Gd 34 %, Tb 24 %, Nd 16 %, Dy 11 %, y pesos menores, algunos negativos, para las demás. Un error del 5 % en el Gd desplaza el Eu/Eu* clásico un 2.4 %, y el de O'Neill, un 1.6 %.
El χ²: la calidad del ajuste en un número
Para resumir los 13 residuos en un solo número, cada residuo se mide en unidades del error analítico esperado σ, se eleva al cuadrado y se suma (Apéndice 1, ec. A3b; es el χ² de la Lección 2):
χ² = Σ ( (yi − ycalc,i) / σ )²
Si la curva es buena y los residuos son solo ruido analítico, cada término vale en promedio 1… salvo que el ajuste ya ha «gastado» un punto por cada λ. De los 13 términos, solo cuentan 13 − 5 = 8, los grados de libertad ν. Por eso se esperan χ² ≈ ν y χ²ν = χ²/ν ≈ 1: el χ² reducido.
La macro usa σ = 1 % (0.01 en ln). Para VG 189 con 5 términos, la suma de residuos al cuadrado es 0.003301, así que χ² = 0.003301 / 0.01² = 33.0 y χ²ν = 33.0 / 8 = 4.1. Un χ²ν de 4 significa que los residuos son unas √4 = 2 veces mayores que el 1 % supuesto: el patrón se aparta de la curva en torno a un 2 %. No es un mal ajuste; es que el 1 % era optimista.
Muchas muestras: la precisión que infiere O'Neill
En una sola muestra, χ²ν fluctúa mucho por azar. Pero en una población de muestras analizadas del mismo modo, la media de χ²ν debe valer 1 si σ es la correcta. Entonces se puede despejar σ: se busca el valor que hace que la media sea 1 (Apéndice 1, ec. A4). O'Neill lo llama s(ln[REE])* inferido. Pruébalo con los 616 vidrios OFB de su hoja:
La σ que da media 1 es 2.14 % (con el ν = 9 de la macro; 2.27 % con ν = 8). Los análisis repetidos de dos patrones analizados en el mismo laboratorio, el vidrio vg2 y el basalto fundido bcr2g, tienen una dispersión de 1.9 ± 0.6 % y 2.1 ± 0.5 % (pp. 1468–1469). Es decir: lo que la cuártica no explica es solo ruido analítico. Los patrones de REE de estos basaltos son suaves hasta donde los análisis permiten verlo. Es la conclusión que da título al artículo.
La distribución sirve además de control. Si σ es la correcta, los χ²ν deben repartirse como la curva, y así ocurre casi siempre. Las pocas muestras de la cola derecha pueden ser fallos analíticos… o desviaciones reales de un patrón suave, petrológicamente interesantes. O'Neill recomienda volver a analizarlas (Apéndice 1).
Esto explica los porcentajes de la Lección 1: al ajustar los 616 vidrios con 2, 3, 4 y 5 términos, la σ inferida baja de 11.7 a 3.9, 2.4 y 2.14 % con el radio, y de 14.0 a 5.6, 3.09 y 2.56 % con Z (p. 1471). Esos porcentajes se reproducen dividiendo por 9 en todos los casos. Con el ν correcto de cada ajuste (13 − número de términos), los valores con el radio serían 10.6, 3.7, 2.4 y 2.27 %: la conclusión no cambia.
Los errores de los λ: los factores ζ
Cada λ es una suma ponderada de los 13 datos (Lección 6): λn = Σ y·fn / Σ fn². Si cada y lleva un error independiente de tamaño s*, el error de esa suma se calcula con la regla habitual de propagación de errores, y queda (Tabla 2 y p. 1466 de O'Neill):
s(λn) = ζn · s*, con ζn = 1 / √Σ fn(ri)²
Es la tabla de errores de la Lección 4 (Bevington, ecs. 7.23–7.27), que con piezas ortogonales solo tiene diagonal: 1/Σ fn² veces σ². Los ζ solo dependen de los radios: son la estación ⑥ del mapa.
| λ | Σ fn² | ζn | s(λ) con 1 % (Tabla 2) | s(λ) con 2.14 % |
|---|---|---|---|---|
| λ0 | 13 | 0.2774 | 0.0028 | 0.0059 |
| λ1 | 0.04722 | 4.602 | 0.046 | 0.098 |
| λ2 | 1.029 × 10⁻⁴ | 98.59 | 0.99 | 2.11 |
| λ3 | 2.550 × 10⁻⁷ | 1980 | 19.8 | 42.4 |
| λ4 | 5.567 × 10⁻¹⁰ | 42 382 | 424 | 907 |
Con estos errores se decide qué λ son significativos: con el criterio habitual (95 %), un λ lo es si se aparta de cero más de dos veces su error. O'Neill da un ejemplo con s* = 5 %: λ3 sería significativo si |λ3| > 2 × 5 × 19.8 = 198, y λ4, si |λ4| > 2 × 5 × 424 = 4240 (p. 1466). Compruébalo con cada muestra y con la precisión que quieras:
| λ | Valor | ± 2·s(λ) | |λ| frente a 2·s(λ) | ¿Significativo? |
|---|
Es lo que concluye O'Neill (Apéndice 1 y pp. 1469–1471). La mayor parte de la información de un patrón está en λ0, λ1 y λ2, más la anomalía de Eu si la hay. λ3 solo se distingue de cero gracias a la precisión del 2 % de los análisis modernos, y λ4 casi nunca.
Una última advertencia: los ζ solo propagan los errores independientes de cada REE. Un error que afecte igual a todas las REE de una muestra solo altera λ0 (Lección 6). O'Neill lo ilustra con 42 lavas de la erupción del Laki de 1783 (p. 1470). La dispersión de λ1…λ4 implica errores de menos del 1 %; la de λ0, un 10 %, debido sobre todo a la cantidad variable de fenocristales, que diluyen todas las REE por igual.
El recorrido completo
Con esta lección se cierra el cálculo de O'Neill:
- Preparar: ppm ÷ CI, ln y eje X = radio iónico; fuera el Eu → 13 puntos (Lección 1).
- Ajustar: mínimos cuadrados con un polinomio (Lecciones 2 y 3), pero con piezas ortogonales, fabricadas solo con los radios (Lecciones 4 y 5).
- Calcular los λ: fórmula directa, sistema de la macro o ESTIMACION.LINEAL (Lección 6).
- Comprobar: cocientes observado/calculado, Eu/Eu* en r = 1.066 Å, χ²ν y, con muchas muestras, s* (esta lección).
- Dar los errores: s(λn) = ζn·s* y decidir qué λ son significativos (esta lección).
Practica
1. En VG 189, y(Eu) = 3.2547 y la curva vale 3.3479 en r = 1.066 Å. ¿Cuánto vale Eu/Eu*?
2. Con σ = 1 %, una muestra da χ²ν = 9. ¿Qué precisión analítica sugiere, en %?
3. Tus análisis tienen una precisión del 3 %. ¿Es significativo un λ3 = −100? (ζ3 = 1980)
4. ¿Por qué ζ4 (42 382) es tanto mayor que ζ1 (4.6)?
Resumen
- La macro evalúa la curva en los 14 radios y da [REE]obs/[REE]calc = exp(y − ycalc). Pinta de rojo lo que se aparta más de un 5 %.
- Eu/Eu* es el cociente del Eu: la curva en r = 1.066 Å da el Eu*. Para VG 189, 0.911, casi igual que el clásico √(Sm·Gd), pero apoyado en todas las REE.
- χ² = Σ(residuo/σ)²; con la σ correcta, χ²ν = χ²/ν ≈ 1, con ν = 13 − 5 = 8 (la macro divide por 9).
- Con muchas muestras se infiere la precisión: s* = 2.14 % para los 616 vidrios OFB, igual que los análisis repetidos. Los patrones son suaves dentro del error analítico.
- s(λn) = ζn·s*, con ζn = 1/√Σ fn², que solo depende de los radios. Un λ es significativo si |λ| > 2·s(λ). λ3 lo es a menudo con un 2 % de precisión; λ4, casi nunca.
Dónde ha aparecido el radio en esta lección: ⑤ la curva se evalúa en los 14 radios; en el del Eu (1.066 Å) da el Eu*, aunque el Eu no entre en el ajuste. Y el Eu2+, con un radio mucho mayor, es la causa de la anomalía. ⑥ Los factores ζn = 1/√Σ fn(ri)² convierten el error analítico en el error de cada λ y solo dependen de los radios. Con esto están recorridas las seis estaciones del mapa.
Aquí termina el núcleo del curso. Quedan dos lecciones opcionales. Lección 8 · Hazlo tú con tus datos: el cálculo completo en Excel o en Python, y qué hacer si falta un elemento o quieres dejar fuera otro (por ejemplo el Ce en circón). Lección 9 · Vectores petrogenéticos ψ: cómo describe O'Neill la fusión y la cristalización con los mismos polinomios. Si quieres alguna, escríbeme «siguiente» o dime cuál.
Fuentes de esta lección
- O'Neill (2016), Apéndice 1 del material suplementario: tres tipos de incertidumbre (ec. A2), el χ² con logaritmos (A3b), cómo inferir s(ln[REE])* (A4), ν = 8 y Fig. A1; exclusión del Eu por el Eu2+; las muestras poco suaves merecen volver a analizarse.
- O'Neill (2016): p. 1466 (ζn, s(λn) = ζn·s* y el criterio de 2 errores, con los ejemplos 198 y 4240), p. 1469 (s* = 2.14 % frente a 2.0 % de vg2 y bcr2g; λ3 y λ4), p. 1470 (suavidad y complejidad; el ejemplo del Laki), p. 1471 (s* con 2–5 términos, r frente a Z; Fig. 3) y Tabla 2 (s(λ) con un 1 %).
- Macro REE_OrthogRegress_5term: curva en los 14 radios (XCalc), Exp(X − XCalc), celdas rojas fuera de 0.95–1.05, SumSq con σ = 0.01 y χ²ν = SumSq/(N − 4).
- Bevington & Robinson (2003): χ² (ec. 7.29), grados de libertad (p. 121) y tabla de errores (ecs. 7.23–7.27).
- Radio del Eu2+ y del Sr2+ en coordinación 8: Shannon (1976). Datos del histograma: χ²ν de los 616 vidrios OFB (Jenner & O'Neill, 2012) en la hoja REEshapes de O'Neill.
Lección 8 · opcional
Hazlo tú con tus datos
Ya conoces todo el cálculo; esta lección es práctica. Hay tres maneras de aplicarlo a tus análisis: la hoja de O'Neill con su macro, una plantilla de Excel sin macros y un script de Python (o pyrolite). Además, una prueba para comprobar que lo haces bien, y qué hacer si te falta un elemento o quieres dejar fuera otro, como el Ce en circón.
- Idea nueva
- Ninguna: aplicar lo aprendido, y saber qué cambia (y qué no) al quitar elementos
- Dónde entra el radio
- Si cambias qué REE entran, cambia la tabla del sistema (④); si además recalculas las constantes, cambian ② y ③
- Archivos
- Lambdas_ONeill_plantilla.xlsx y lambdas_oneill.py, que acompañan a esta lección en la conversación
La receta y la prueba de control
- Concentraciones en ppm de las 14 REE. Los valores por debajo del límite de detección se dejan en blanco, nunca a 0.
- Divide por el condrito CI de O'Neill (Tabla 1) y toma el logaritmo neperiano (LN en Excel, no LOG).
- Eje X: radio iónico del catión 3+ en coordinación 8 (Shannon, 1976), el de la Tabla 1.
- Deja fuera el Eu (y el Ce, u otras, si tienes motivo).
- Ajusta con las piezas de la Tabla 2: con la macro, con la plantilla o con el script.
- Revisa los cocientes observado/calculado (y las anomalías de las REE excluidas) y el χ² reducido.
- Da los errores de los λ con s(λn) = ζn·s*.
Antes de usar cualquier herramienta con tus datos, compruébala con VG 189. Estas son sus concentraciones en ppm:
| La | Ce | Pr | Nd | Sm | Eu | Gd | Tb | Dy | Ho | Er | Tm | Yb | Lu |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3.622 | 12.477 | 2.1646 | 12.125 | 4.4536 | 1.5340 | 5.7578 | 1.0734 | 6.8684 | 1.4817 | 4.2397 | 0.6578 | 3.9963 | 0.5711 |
Debes obtener λ = 3.187, −1.656, −50.60, −227.29 y −2001.8, Eu/Eu* = 0.911 y χ²ν = 4.13 con σ = 1 % (3.667 si divides por N − 4, como la macro). Si no te sale, revisa el logaritmo, el condrito y los radios antes que nada.
Camino A: la hoja de O'Neill con su macro
Es la opción más directa si trabajas en Excel y puedes ejecutar macros. En la hoja «Data In & Out» del archivo LstSq_REE_OrthogPolynomial_April2015.xls (instrucciones en la hoja «Instructions»):
- Borra las celdas de entrada y de salida desde la fila 6.
- Pega tus datos en ppm desde B6 (una muestra por fila, de La en B a Lu en O; el nombre, en A) y pon un 1 en B3 para normalizar al condrito.
- Pon en B4 el número de la última fila con datos, o uno algo mayor.
- Ejecuta la macro REE_OrthogRegress_5term.
| Columnas | Contenido |
|---|---|
| A · B–O | Nombre · La…Lu en ppm (entrada) |
| Q–U | λ0…λ4 |
| V | Número de REE usadas en el ajuste |
| W | χ² reducido con σ = 1 %, dividido por N − 4 (Lección 7) |
| Y–AL | [REE]obs/[REE]calc de La a Lu; la columna AD (Eu) es Eu/Eu*. En rojo, las que se apartan más de un 5 % |
| Filas 2–3 | W2: s* inferido de todas las filas (2.14 % con sus datos); Q3–U3: s(λ) para un 1 %; Y2–AL3: media y desviación típica de cada cociente |
Tres cosas que conviene saber. Si falta una REE, basta dejarla en blanco: la macro la salta y mantiene las constantes de la Tabla 2. Las propias instrucciones avisan de que entonces los λ «no serán del todo independientes». Una fila se ignora si le faltan a la vez La y Ce. Y para excluir el Ce hay que dejarlo en blanco, así que la macro no dará su cociente: Ce/Ce* hay que calcularlo aparte, con la curva en r = 1.143 Å (lo hacen la plantilla y el script).
Camino B: la plantilla de Excel sin macros
Lambdas_ONeill_plantilla.xlsx hace lo mismo que la macro solo con fórmulas, sin macros ni fórmulas matriciales: funciona en Excel y en LibreOffice. Tiene cuatro hojas:
- Muestra: una muestra con todos los pasos a la vista (y, piezas, tabla del sistema, su inversa, λ, curva, cocientes, χ²). Cada REE tiene su casilla «¿Entra?» con 1 o 0.
- Lote: hasta 200 muestras, una por fila, todas con las mismas REE en el ajuste (fila 6). Da N, λ0…λ4, χ²ν, s y los 14 cocientes, y arriba los errores de los λ.
- Constantes y Lote_calculos: radios, CI, constantes de la Tabla 2 y cálculos intermedios. No hace falta tocarlas.
El truco de la hoja Lote es el de la Lección 6: la tabla del sistema solo depende de los radios de las REE que entran (④), así que se calcula e invierte una sola vez para todo el lote. Cada muestra solo aporta su vector de sumas con los datos. Por eso todas las filas deben usar las mismas REE. Si a una muestra le falta una de las que entran, la fila lo avisa y no da λ; esa muestra se calcula en la hoja Muestra.
La plantilla resuelve G·λ = h con G = Σ w·fj·fk y h = Σ w·fj·y (w = 1 si la REE entra, 0 si no), usando INDEX(MINVERSA(G); fila; columna), que no necesita fórmulas matriciales. Con las 13 REE, G es diagonal y el resultado es la fórmula directa de la Lección 6; si falta alguna, es el mismo ajuste que hace la macro.
Camino C: Python
El script lambdas_oneill.py solo necesita numpy y pandas. Desde la línea de órdenes:
python lambdas_oneill.py --prueba # comprueba VG 189 contra O'Neill
python lambdas_oneill.py mis_datos.csv # CSV con columnas La…Lu en ppm
python lambdas_oneill.py mis_datos.csv --excluir Eu Ce # circón: fuera Eu y Ce
python lambdas_oneill.py mis_datos.csv --constantes recalcular --terminos 4O desde tu propio código (el cálculo central son cinco líneas):
F = piezas(radios_que_entran, TABLA2) # f0…f4 en cada radio (solo radios)
G = F @ F.T # tabla del sistema
h = F @ y # sumas con los datos
lam = np.linalg.solve(G, h) # λ0…λ4
obs_calc = np.exp(y_14 - lam @ piezas(radios_14, TABLA2)) # incluye Eu/Eu*, Ce/Ce*Ver el script completo (para copiarlo)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
lambdas_oneill.py — Coeficientes λ de O'Neill (2016) para patrones de REE.
Reproduce el cálculo de la macro REE_OrthogRegress_5term de la hoja de O'Neill:
1. y = ln([REE] / [REE]CI), con los valores CI de O'Neill (2016, Tabla 1).
2. Eje X: radio iónico del catión 3+ en coordinación 8 (Shannon, 1976).
3. Ajuste por mínimos cuadrados a λ0 + λ1·f1 + λ2·f2 + λ3·f3 + λ4·f4, con los
polinomios ortogonales de la Tabla 2 de O'Neill (u otros, ver «constantes»).
4. Curva calculada en los 14 radios: cocientes observado/calculado
(en los elementos excluidos del ajuste, son las anomalías: Eu/Eu*, Ce/Ce*…).
5. χ² reducido (σ = 1 % por defecto) y errores de los λ.
Solo necesita numpy y pandas.
USO DESDE LA LÍNEA DE ÓRDENES
python lambdas_oneill.py --prueba
Calcula VG 189 y lo compara con los valores de O'Neill.
python lambdas_oneill.py datos.csv --salida resultados.csv
datos.csv: una fila por muestra; columnas La, Ce, …, Lu en ppm (se admiten
otras columnas, por ejemplo el nombre de la muestra). Separador: coma o punto y coma.
Opciones:
--excluir Eu Ce REE que no entran en el ajuste (por defecto: Eu)
--terminos 5 número de términos (1 a 5; por defecto 5, hasta λ4)
--constantes tabla2 'tabla2' (O'Neill; recomendado) o 'recalcular'
--sigma 0.01 error supuesto de cada ln([REE]) para el χ² (O'Neill: 0.01)
--s 0.02 precisión analítica para los errores de los λ
USO COMO MÓDULO
from lambdas_oneill import ajustar, ajustar_tabla
r = ajustar({'La': 3.622, 'Ce': 12.477, …}, excluir=('Eu', 'Ce'))
r['lambdas'], r['obs_calc']['Ce'], r['chi2_red']
Equivalente con pyrolite (versión 0.3.7), para obtener los mismos λ:
import pyrolite.geochem
df.pyrochem.lambda_lnREE(params="ONeill2016", degree=5,
exclude=["Pm", "Eu"], sigmas=0.01, add_X2=True)
· params="ONeill2016": por defecto pyrolite usa una base calculada con las 14 REE,
incluido el Eu, y sus λ no coinciden exactamente con los de O'Neill.
· degree=5: en pyrolite, degree es el número de términos (degree=4 da λ0–λ3).
· Pasa siempre una lista nueva en exclude: la lista por defecto se modifica al usar
anomalies=[…] y el cambio persiste en la sesión.
· Una fila con algún valor ≤ 0 se descarta entera: convierte antes esos valores en NaN.
"""
import argparse
import sys
import numpy as np
import pandas as pd
REE = ['La', 'Ce', 'Pr', 'Nd', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu']
# Radios iónicos (Å) del catión 3+ en coordinación 8: Shannon (1976), Tabla 1 de O'Neill.
RADIOS = dict(zip(REE, [1.160, 1.143, 1.126, 1.109, 1.079, 1.066, 1.053,
1.040, 1.027, 1.015, 1.004, 0.994, 0.985, 0.977]))
# Condrito CI (ppm): O'Neill (2016), Tabla 1 (los mismos que usa su macro).
CI = dict(zip(REE, [0.2472, 0.6308, 0.0950, 0.4793, 0.15419, 0.0592, 0.2059,
0.0375, 0.2540, 0.0554, 0.1645, 0.0258, 0.1684, 0.0251]))
# Constantes de los polinomios ortogonales (Tabla 2 de O'Neill), con los decimales de su macro.
# Cada pieza es f_n(r) = producto de (r − constante).
TABLA2 = [
(), # f0 = 1
(1.054769,), # f1: β
(1.005327429, 1.128236038), # f2: γ1, γ2
(0.991412204, 1.060548105, 1.145519887), # f3: δ1, δ2, δ3
(0.984820219, 1.030518142, 1.104414973, 1.153426708), # f4: ε1…ε4
]
def piezas(r, constantes=TABLA2, terminos=5):
"""Matriz (terminos × len(r)) con los valores de f0…f_{terminos-1} en los radios r."""
r = np.atleast_1d(np.asarray(r, dtype=float))
return np.array([np.prod([r - c for c in constantes[n]], axis=0) if n else np.ones_like(r)
for n in range(terminos)])
def constantes_ortogonales(radios, terminos=5):
"""Constantes (ceros) de los polinomios ortogonales para un conjunto de radios.
Con los 13 radios de O'Neill (sin Eu) reproduce la Tabla 2. Usa la recurrencia de
tres términos, equivalente a la receta de Bevington & Robinson (2003, §7.3).
"""
x = np.asarray(radios, dtype=float)
p_prev, p = np.zeros_like(x), np.ones_like(x)
alfa, beta = [], [0.0]
for n in range(terminos - 1):
alfa.append(np.sum(x * p * p) / np.sum(p * p))
if n > 0:
beta.append(np.sum(p * p) / np.sum(p_prev * p_prev))
p_prev, p = p, (x - alfa[n]) * p - beta[n] * p_prev
def f(n, t):
a, b = np.ones_like(t), t - alfa[0]
if n == 0:
return a
for k in range(1, n):
a, b = b, (t - alfa[k]) * b - beta[k] * a
return b
# Los ceros de f_n están dentro del intervalo de los radios: se buscan por cambio de signo.
t = np.linspace(x.min() - 0.05, x.max() + 0.05, 20001)
consts = [()]
for n in range(1, terminos):
v = f(n, t)
idx = np.where(np.sign(v[:-1]) * np.sign(v[1:]) < 0)[0]
ceros = []
for i in idx: # refinamiento por bisección
lo, hi = t[i], t[i + 1]
for _ in range(60):
mid = 0.5 * (lo + hi)
if np.sign(f(n, np.array([lo]))[0]) * np.sign(f(n, np.array([mid]))[0]) <= 0:
hi = mid
else:
lo = mid
ceros.append(0.5 * (lo + hi))
consts.append(tuple(ceros))
return consts
def ajustar(ppm, excluir=('Eu',), terminos=5, constantes='tabla2', sigma=0.01, s_est=0.02):
"""Ajusta una muestra. ppm: dict o Series con las REE en ppm (faltan = NaN o ausentes).
constantes: 'tabla2' (O'Neill, recomendado para comparar con otros trabajos),
'recalcular' (polinomios ortogonales para las REE que entran), o una lista
de tuplas como TABLA2.
Devuelve un dict con: lambdas, errores (con s_est), n_ree, chi2_red (ν = N − términos),
chi2_red_macro (÷ (N − 4), como la macro con 5 términos), s_muestra y obs_calc
(cocientes observado/calculado en las 14 REE; en las excluidas, son la anomalía).
"""
ppm = pd.Series(ppm, dtype=float).reindex(REE)
y = pd.Series(np.nan, index=REE)
validos = ppm.notna() & (ppm > 0)
y[validos] = np.log(ppm[validos] / pd.Series(CI)[validos])
usar = [e for e in REE if e not in excluir and validos[e]]
if len(usar) < terminos + 1:
raise ValueError('Hacen falta al menos %d REE válidas para %d términos.' % (terminos + 1, terminos))
x = np.array([RADIOS[e] for e in usar])
yy = y[usar].values
if constantes == 'tabla2':
consts = TABLA2
elif constantes == 'recalcular':
consts = constantes_ortogonales(x, terminos)
else:
consts = constantes
F = piezas(x, consts, terminos) # terminos × N
G = F @ F.T # tabla del sistema (solo radios)
h = F @ yy # Σ f_n · y
lam = np.linalg.solve(G, h) # = G⁻¹ h (con las 13 REE, λn = Σ y fn / Σ fn²)
Ginv = np.linalg.inv(G)
todos = np.array([RADIOS[e] for e in REE])
ycalc = lam @ piezas(todos, consts, terminos) # curva en los 14 radios
res = yy - lam @ F
N = len(usar)
ss = float(np.sum(res ** 2))
obs_calc = pd.Series(np.exp(y.values - ycalc), index=REE)
corr = Ginv / np.sqrt(np.outer(np.diag(Ginv), np.diag(Ginv)))
return {
'lambdas': pd.Series(lam, index=['λ%d' % n for n in range(terminos)]),
'errores': pd.Series(s_est * np.sqrt(np.diag(Ginv)), index=['s(λ%d)' % n for n in range(terminos)]),
'corr_max': float(np.max(np.abs(corr - np.eye(terminos)))) if terminos > 1 else 0.0,
'n_ree': N,
'ree_usadas': usar,
'chi2_red': ss / sigma ** 2 / (N - terminos) if N > terminos else np.nan,
'chi2_red_macro': ss / sigma ** 2 / (N - 4) if (terminos == 5 and N > 4) else np.nan,
's_muestra': np.sqrt(ss / (N - terminos)) if N > terminos else np.nan,
'obs_calc': obs_calc,
'constantes': consts,
}
def ajustar_tabla(df, **kw):
"""Aplica ajustar() a cada fila de un DataFrame con columnas La…Lu en ppm."""
filas = []
for idx, fila in df.iterrows():
try:
r = ajustar(fila[[c for c in REE if c in fila.index]], **kw)
d = dict(r['lambdas'])
d.update(r['errores'])
d.update({'N': r['n_ree'], 'chi2_red': r['chi2_red'], 'chi2_red_macro': r['chi2_red_macro'],
's_muestra': r['s_muestra'], 'corr_max': r['corr_max']})
d.update({'%s obs/calc' % e: v for e, v in r['obs_calc'].items()})
except ValueError as err:
d = {'aviso': str(err)}
filas.append(pd.Series(d, name=idx))
return pd.DataFrame(filas)
VG189 = dict(zip(REE, [3.622024, 12.477229, 2.164558, 12.124874, 4.45359, 1.533984, 5.757827,
1.073365, 6.868444, 1.481695, 4.239659, 0.657756, 3.996338, 0.571062]))
def prueba():
"""Comprueba el cálculo con VG 189 (primera fila de la hoja de O'Neill)."""
r = ajustar(VG189)
ref = [3.1875, -1.6558, -50.600, -227.29, -2001.8]
print('VG 189 (13 REE, Tabla 2, 5 términos)')
for (k, v), rv in zip(r['lambdas'].items(), ref):
print(' %-3s = %11.4f (O\'Neill: %s)' % (k, v, rv))
print(' χ²ν (σ = 1 %%, ν = %d) = %.3f; como la macro (÷ %d) = %.3f (O\'Neill: 3.667)'
% (r['n_ree'] - 5, r['chi2_red'], r['n_ree'] - 4, r['chi2_red_macro']))
print(' Eu/Eu* = %.3f (O\'Neill: 0.911)' % r['obs_calc']['Eu'])
r2 = ajustar(VG189, excluir=('Eu', 'Ce'))
r3 = ajustar(VG189, excluir=('Eu', 'Ce'), constantes='recalcular')
print('Sin Ce, Tabla 2: ', np.round(r2['lambdas'].values, 4), ' Ce/Ce* = %.3f' % r2['obs_calc']['Ce'])
print('Sin Ce, recalculadas:', np.round(r3['lambdas'].values, 4), ' Ce/Ce* = %.3f' % r3['obs_calc']['Ce'])
ok = np.allclose(r['lambdas'].values, ref, rtol=2e-4, atol=1e-3)
print('Resultado:', 'coincide con O\'Neill' if ok else 'NO coincide: revisa los datos')
return ok
def main(argv=None):
ap = argparse.ArgumentParser(description="Coeficientes λ de O'Neill (2016) para patrones de REE.")
ap.add_argument('datos', nargs='?', help='CSV con columnas La…Lu en ppm')
ap.add_argument('--salida', default='resultados_lambdas.csv')
ap.add_argument('--excluir', nargs='*', default=['Eu'])
ap.add_argument('--terminos', type=int, default=5)
ap.add_argument('--constantes', choices=['tabla2', 'recalcular'], default='tabla2')
ap.add_argument('--sigma', type=float, default=0.01)
ap.add_argument('--s', type=float, default=0.02, dest='s_est')
ap.add_argument('--prueba', action='store_true')
a = ap.parse_args(argv)
if a.prueba or not a.datos:
return 0 if prueba() else 1
df = pd.read_csv(a.datos, sep=None, engine='python')
faltan = [e for e in REE if e not in df.columns]
if faltan:
print('Aviso: faltan columnas', faltan, '(se tratan como datos ausentes)')
res = ajustar_tabla(df, excluir=tuple(a.excluir), terminos=a.terminos,
constantes=a.constantes, sigma=a.sigma, s_est=a.s_est)
pd.concat([df, res], axis=1).to_csv(a.salida, index=False)
print('Escrito:', a.salida, '(%d muestras)' % len(df))
return 0
if __name__ == '__main__':
sys.exit(main())
Con pyrolite. Si ya usas pyrolite (versión 0.3.7), esta llamada reproduce exactamente los λ de O'Neill:
import pyrolite.geochem
ls = df.pyrochem.lambda_lnREE(params="ONeill2016", degree=5,
exclude=["Pm", "Eu"], sigmas=0.01, add_X2=True)Al revisar su código he encontrado cuatro detalles que cambian los resultados:
- params="ONeill2016". Por defecto (params=None), pyrolite fabrica las piezas con las 14 REE, incluido el radio del Eu, aunque luego no ajuste el Eu. Los λ cambian: para VG 189, λ0 = 3.199 y λ1 = −1.616 en vez de 3.187 y −1.656 (λ0 se desvía el doble de su error analítico).
- degree=5. En pyrolite, degree es el número de términos: degree=4 da λ0…λ3.
- exclude, siempre explícito y nuevo. La lista por defecto se modifica dentro de la función: una llamada con anomalies=["Ce", "Eu"] deja el Ce excluido en todas las llamadas siguientes de la sesión.
- Ceros. Una fila con algún valor ≤ 0 se descarta entera; convierte antes esos valores en NaN.
Además, los errores que da add_uncertainties=True están escalados por el χ²ν de cada muestra: usan la dispersión de esa muestra, no la precisión común s* de O'Neill. Ambas opciones son legítimas, pero no son lo mismo.
Si falta un elemento o quieres quitar otro
Hay dos opciones: mantener las piezas de la Tabla 2, como recomienda O'Neill y hacen su macro, la plantilla y pyrolite con params="ONeill2016", o recalcular las constantes para las REE que entran. Compáralas:
| λ | λ0 | λ1 | λ2 | λ3 | λ4 |
|---|
Lo que no cambia: la curva. Con las mismas REE y el mismo número de términos, los dos juegos de piezas generan los mismos polinomios, y el ajuste por mínimos cuadrados es único. Los cocientes, el χ² y las anomalías (Ce/Ce*, Eu/Eu*) son idénticos. Recalcular las constantes no «mejora» el ajuste.
Lo que cambia: los λ. Son las mismas curvas repartidas en piezas distintas. Con las piezas de la Tabla 2, los λ conservan su significado y se pueden comparar con los de O'Neill y con los de otras muestras. A cambio, si faltan varias REE, se correlacionan y pierden precisión. Sin el Ce, apenas se nota: los λ de VG 189 cambian menos de una desviación típica y la correlación máxima es 0.16. Sin La, Ce y Pr, la correlación sube a 0.996 y el error de λ1 se multiplica por 19. Con constantes recalculadas, los λ vuelven a ser independientes, pero describen otra cosa. Sin La, Ce y Pr, λ1 ya no es la pendiente de todo el patrón, sino la de Nd a Lu, y para VG 189 cambia de signo (+1.06 en vez de −1.66). Por eso O'Neill insiste en usar siempre los mismos polinomios (p. 1466). Y solo tiene sentido comparar λ de muestras con las mismas REE en el ajuste.
Cuidado con los extremos. Una REE excluida en medio del patrón (el Eu, el Ce si están el La y el Pr) queda entre datos: su anomalía se interpola y es robusta. Una excluida en un extremo se extrapola, y eso es muy inestable. Prueba «Eu, Ce, La y Pr» con VG 189, que no tiene anomalía de Ce. Con 2, 3, 4 y 5 términos, Ce/Ce* sale 0.67, 0.92, 0.99 y 1.47. Con el La y el Pr dentro sale 1.07, 1.06 y 1.04 con 3, 4 y 5 términos.
Errores frecuentes
- Logaritmo decimal en vez de neperiano: todos los λ salen divididos por ln 10 = 2.303.
- Otro condrito (por ejemplo McDonough & Sun, 1995): cambian λ0 y también la forma, porque los cocientes entre REE de los dos condritos no son idénticos. Solo son comparables los λ calculados con el mismo condrito (en pyrolite, «ChondriteREE_ON» es el de O'Neill).
- Otros radios (coordinación 6 o 9) o el número atómico Z: cambian todas las constantes.
- Incluir el Eu en el ajuste, o incluir una REE con una anomalía evidente.
- Ceros o valores bajo el límite de detección tratados como números: ln(0) no existe. Déjalos en blanco.
- Constantes redondeadas al convertir los a de un ajuste normal, u olvidar que ESTIMACION.LINEAL devuelve los coeficientes al revés (Lecciones 3 y 6).
- Comparar χ²ν de programas distintos sin mirar el ν: la macro divide por N − 4 y pyrolite, el script y la plantilla, por N − 5.
Practica
1. Calculas los λ de VG 189 con pyrolite, sin más opciones, y obtienes λ0 = 3.199 en vez de 3.187. ¿Qué ha pasado?
2. Alguien calcula los λ con log10 en vez de ln. ¿Por qué número habría que multiplicar sus λ para obtener los de O'Neill?
3. En un circón has quitado del ajuste el La, el Ce, el Pr y el Eu, y con 5 términos obtienes Ce/Ce* = 1.5. ¿Qué haces antes de interpretarlo?
Resumen
- Tres caminos: la macro de O'Neill (pega en B6:O…, pon 1 en B3 y la última fila en B4), la plantilla de Excel sin macros y el script de Python; o pyrolite con params="ONeill2016", degree=5 y exclude explícito.
- Comprueba siempre con VG 189: 3.187, −1.656, −50.60, −227.29, −2001.8.
- Si falta o quitas una REE, la curva no depende de las constantes: cocientes, anomalías y χ² son iguales con la Tabla 2 o con constantes recalculadas. Cambian los λ.
- Con la Tabla 2 los λ son comparables, pero se correlacionan si faltan varias REE. Con constantes recalculadas son independientes, pero describen otra cosa. O'Neill recomienda la Tabla 2.
- Las anomalías de REE en un extremo del patrón se extrapolan: comprueba cómo cambian con el número de términos.
Última lección opcional · Vectores petrogenéticos ψ. Veremos por qué O'Neill escribe λ = λ⁰ + ψ: los coeficientes de reparto de un mineral también se describen con los mismos polinomios en el radio, y los procesos (fusión, cristalización) se suman como flechas en el diagrama λ2–λ1.
Fuentes de esta lección
- O'Neill (2016): Tablas 1 y 2, p. 1466 (usar los mismos polinomios aunque falte algún elemento) y la hoja LstSq_REE_OrthogPolynomial_April2015.xls (hoja «Instructions», macro REE_OrthogRegress_5term y distribución de las columnas).
- pyrolite 0.3.7 (Williams et al., 2020, Journal of Open Source Software 5, 2314): código de lambda_lnREE, calc_lambdas y lambdas_ONeill2016, comprobado aquí con VG 189 y las otras tres muestras. Documentación de los λ en pyrolite.
- Anenburg, M. & Williams, M. J. (2022). Quantifying the tetrad effect, shape components, and Ce–Eu–Gd anomalies in rare earth element patterns. Mathematical Geosciences 54, 47–70. Artículo · ALambdaR.
- Ernst, D. M., Vogt, J., Bau, M. & Mues, M. (2025). Polynomial modelling of high-quality yet incomplete rare earth element data sets and a holistic assessment of REE anomalies. Scientific Reports. Artículo.
- Cálculos propios: plantilla y script verificados con VG 189 (13 REE, sin Ce y sin Lu) frente a la macro de O'Neill y a pyrolite; comparación de constantes fijas y recalculadas para varios conjuntos de REE.
Lección 9 · opcional
Vectores petrogenéticos ψ
Hasta ahora, los λ describían la forma de un patrón. En esta última lección veremos que también describen los procesos que cambian esa forma: la cristalización de un mineral, la fusión de la fuente o la mezcla de dos magmas. O'Neill lo resume en una suma, λ = λ⁰ + ψ, que en el diagrama λ2–λ1 se convierte en una flecha.
- Idea nueva
- Un proceso también tiene «forma»: sus coeficientes ψ se suman a los λ del magma inicial
- Dónde entra el radio
- ⑦ Los coeficientes de reparto D dependen del radio y se describen con las mismas piezas fn(r) (Tabla 7)
- Del artículo
- pp. 1485–1494: ecuaciones 2 a 8, Tablas 5 a 7 y Fig. 17
Un proceso se suma en logaritmos
Cuando un magma cristaliza, o cuando una roca funde, la concentración de cada REE en el líquido queda multiplicada por un factor. Por ejemplo, si la mitad de un magma cristalizara como clinopiroxeno, en el líquido restante el La se multiplicaría por 1.88 y el Lu por 1.38: el clinopiroxeno acepta el Lu mejor que el La. El factor depende del mineral y de cuánto ha cristalizado, pero no de cuánto La o Lu había al principio.
En logaritmos, multiplicar se convierte en sumar (lo anunciamos en la Lección 1):
[REE]⁰ es la concentración inicial: la del magma de partida o la de la fuente que funde. O'Neill llama al último término función del proceso petrogenético. Piensa en una cuenta bancaria: el saldo final es el saldo inicial más los movimientos, y los movimientos no dependen de lo que hubiera en la cuenta. Aquí, el saldo de cada REE es su y = ln([REE]/CI), y los movimientos, la función del proceso.
La función del proceso también cambia de forma suave con el radio (enseguida veremos por qué), así que se puede describir con las mismas piezas fn de la Tabla 2. Sus coeficientes se llaman ψ (psi) y forman el vector del proceso, Ψ = {ψ0, ψ1, ψ2…}. Es la ecuación 3 de O'Neill:
Los ψ se calculan igual que los λ: son la «sombra» de la función del proceso sobre cada pieza, ψn = Σ g·fn / Σ fn² (Lección 6), donde g es el valor de la función del proceso en cada REE. Y como la sombra de una suma es la suma de las sombras, los λ del magma final son los del inicial más los ψ, orden por orden:
Cada orden va por su lado: λ0 = λ0⁰ + ψ0, λ1 = λ1⁰ + ψ1, λ2 = λ2⁰ + ψ2… Sumar funcionaría con cualquier ajuste por mínimos cuadrados. Lo que aportan las piezas ortogonales es que cada ψn no depende de cuántos términos uses ni de los demás. Por eso puedes trabajar solo con λ1 y λ2 sin preocuparte de λ3 y λ4.
Compruébalo con números. Tomemos VG 304 como magma inicial y hagamos cristalizar la mitad como clinopiroxeno (coeficientes de O'Neill & Jenner, 2012). Multiplicamos cada REE por su factor, recalculamos los λ del patrón nuevo con la fórmula de la Lección 6 y los comparamos con la suma:
| λ0 | λ1 | λ2 | |
|---|---|---|---|
| VG 304, magma inicial (λ⁰) | 2.653 | −5.866 | −47.74 |
| Proceso: 50 % de clinopiroxeno (ψ) | +0.388 | +2.024 | +14.35 |
| Suma λ⁰ + ψ | 3.041 | −3.842 | −33.39 |
| Patrón final, ajustado de nuevo | 3.041 | −3.842 | −33.39 |
Coinciden en todas las cifras, y también en λ3 y λ4. De aquí salen tres consecuencias que usa O'Neill:
- La parte (ψ1, ψ2) es una flecha en el diagrama λ2–λ1, con una dirección y una longitud definidas: el vector del proceso.
- La flecha no depende del magma inicial, así que se puede trasladar a cualquier punto del diagrama.
- Los procesos sucesivos (fundir y luego cristalizar, por ejemplo) se encadenan como flechas, una tras otra.
Un matiz: la ortogonalidad de las piezas de la Tabla 2 es exacta con las 13 REE. Para los coeficientes de reparto casi nunca hay datos de las 13, pero según O'Neill la diferencia es despreciable si cubren del La al Lu (o al menos hasta el Yb).
Los coeficientes de reparto también son curvas en el radio
¿Por qué es suave la función del proceso? Porque depende de los coeficientes de reparto. El coeficiente de reparto de un elemento entre un mineral y el fundido, D, es el cociente entre su concentración en el mineral y en el líquido. Con D = 0.1, la concentración en el cristal es la décima parte que en el líquido: el elemento es incompatible y se acumula en el fundido. Con D = 5, el cristal tiene cinco veces más: lo acapara.
La teoría de la deformación de la red cristalina (lattice strain; Wood & Blundy, 2014) explica cómo depende D del radio. Cada posición del cristal donde pueden entrar las REE es como un zapato de una talla r0: un ion de ese tamaño entra cómodo; uno mayor o menor deforma la red y le cuesta más, tanto más cuanto más se aparte de r0. Es la ecuación 7 de O'Neill:
donde B mide lo rígida que es esa posición del cristal. El resultado es una curva suave en forma de campana, y como los radios de las REE cambian poco a poco del La al Lu, D también. En el granate, por ejemplo, la talla es pequeña: las REE pesadas, de radio pequeño, entran mucho mejor que las ligeras (DLu = 7.1 frente a DLa = 0.003). Es la estación ⑦ del mapa del radio.
O'Neill ajusta directamente D (no ln D) con las mismas piezas de la Tabla 2, y comprueba que bastan unos pocos términos para reproducir los datos dentro de su incertidumbre. Es su ecuación 6, y los coeficientes están en su Tabla 7:
Igual que con los λ, δ0 es el valor medio de D en las 13 REE, δ1 su pendiente y δ2 su curvatura. Ojo con los nombres: O'Neill reutiliza la letra δ, pero estos δ0…δ4 no tienen nada que ver con las constantes δ1, δ2 y δ3 de f3 (Lección 5).
| Mineral | Datos de | δ0 | δ1 | δ2 | DLa | DLu |
|---|---|---|---|---|---|---|
| Olivino | O'Neill & Jenner (2012) | 0.01 | −0.10 | 1.4 | 10⁻⁵ | 0.023 |
| Plagioclasa | O'Neill & Jenner (2012) | 0.02 | 0.10 | 1.3 | 0.033 | 0.014 |
| Ortopiroxeno | Yao et al. (2012) | 0.04 | −0.55 | 3.8 | 0.001 | 0.102 |
| Clinopiroxeno 1 | O'Neill & Jenner (2012) | 0.44 | −2.92 | −20.7 | 0.081 | 0.528 |
| Clinopiroxeno 2 | Sun & Liang (2012)* | 0.39 | −3.00 | −4.2 | 0.077 | 0.575 |
| Clinopiroxeno | Yao et al. (2012) | 0.30 | −2.44 | −6.6 | 0.047 | 0.433 |
| Anfíbol | Tiepolo et al. (2000)† | 0.89 | −4.09 | −61.3 | 0.162 | 0.949 |
| Granate | Van Westrenen & Draper (2007)* | 2.35 | −35.4 | 365 | 0.003 | 7.081 |
| Granate | McKenzie & O'Nions (1991) | 2.49 | −39.7 | 496 | 0.001 | 8.552 |
Extracto de la Tabla 7 de O'Neill (también da δ3 y δ4). * Aplicados a las fases del experimento 40.06 de Walter (1998), a 40 kbar. † Normalizados a DYb = 1; la Fig. 17a usa DYb = 0.7. La ecuación 7 relaciona ln D, y no D, con el radio; O'Neill comprueba que, aun así, D se ajusta bien con pocos términos. Le conviene ajustar D porque, en la cristalización fraccionada, D aparece sin logaritmo, como vas a ver ahora.
Cristalización fraccionada: flechas rectas
En la cristalización fraccionada simple (SFX en el artículo), los cristales se separan del líquido a medida que se forman. La ecuación de Rayleigh da la función del proceso (ecuación 4 de O'Neill):
F es la fracción de líquido que queda (F = 0.8 si ha cristalizado un 20 %). Como ln F es negativo, las REE con D < 1 aumentan en el líquido y las que tienen D > 1 disminuyen. Si cristalizan varios minerales a la vez, D es la media de sus D, pesada por la proporción de cada mineral, mX (ecuación 5).
Aquí está el atajo. Como D es una suma de piezas, δ0 + δ1·f1 + δ2·f2…, la función del proceso también lo es, y los ψ salen sin ajustar nada:
Con varios minerales, cada δ es la media pesada Σ mX·δX. De aquí salen dos propiedades muy útiles. La dirección de la flecha en el diagrama λ2–λ1 es ψ2/ψ1 = δ2/δ1: la fija el mineral, no la cantidad que cristaliza. Y su longitud es proporcional a ln F. Explóralo:
Coeficientes de reparto D
Función del proceso: ln F·(D − 1)
Diagrama λ2–λ1
Lo que muestra la figura, y O'Neill subraya (p. 1487):
- Solo el clinopiroxeno, el anfíbol y el granate cambian la forma de manera apreciable: son los únicos con δ1 y δ2 grandes. Aun así, hacen falta cantidades enormes. Con un 20 % de clinopiroxeno, una cantidad realista, la flecha mide (0.65, 4.6) y no se notaría entre la dispersión de los OFB. O'Neill usa un 80 % solo para que se vea.
- El granate va en otra dirección, casi perpendicular a la del clinopiroxeno y el anfíbol: baja λ2 (más joroba) porque se lleva las REE pesadas. Es la «firma del granate».
- El olivino y la plagioclasa apenas cambian la forma, solo λ0. Sus D son casi 0, así que la función del proceso es casi igual en todas las REE: −ln F. Con un 20 % de olivino, λ0 sube 0.22, mientras que λ1 cambia 0.02 y λ2, −0.3. Por eso, cuando un basalto evoluciona por cristalización fraccionada en un sistema cerrado, λ0 sube a medida que baja el MgO mientras la forma apenas cambia (Fig. 18a-b de O'Neill; Lección 6).
Lo mismo ocurre al revés: si una roca acumula cristales sin REE (olivino y, casi, plagioclasa), el líquido se diluye por igual en todas las REE y solo cambia λ0, en el logaritmo de la fracción de líquido que hay en la mezcla (ecuación 8 de O'Neill). La acumulación de olivino o de plagioclasa no mueve una muestra en el diagrama λ2–λ1.
En los volcanes, el magma no suele cristalizar en un sistema cerrado: la cámara se rellena, entra en erupción y cristaliza muchas veces (replenish–tap–fractionate, RTX). O'Neill muestra que ese proceso amplifica el efecto del clinopiroxeno (Fig. 17b) y propone que los patrones «en plato» de algunas boninitas podrían deberse a un RTX con anfíbol (p. 1490).
Fusión parcial: trayectorias curvas
Para la fusión, O'Neill usa el modelo más sencillo, la fusión en equilibrio (batch melting; ecuación 2):
Aquí F es la fracción de fundido y [REE]⁰, la concentración en la fuente. D⁰ es el D medio de la fuente, pesado por su moda (la proporción de cada mineral). P es el D medio de lo que funde, pesado por la reacción de fusión (qué proporción de cada mineral pasa al líquido). Las modas y reacciones que usa O'Neill están en su Tabla 5:
| Fuente | Olivino | Opx | Cpx | Espinela o granate | Referencia | |
|---|---|---|---|---|---|---|
| Lherzolita con espinela, 10 kbar | Moda | 0.51 | 0.275 | 0.19 | 0.025 | O'Neill (2016) |
| Reacción | −0.167 | 0.652 | 0.466 | 0.049 | Niu (1997) | |
| Lherzolita con granate, 30 kbar | Moda | 0.52 | 0.145 | 0.30 | 0.035 | Walter (1998) |
| Reacción | 0.05 | −0.15 | 0.96 | 0.14 | Walter (2014) | |
| Lherzolita con granate, 40 kbar | Moda | 0.53 | 0.03 | 0.30 | 0.14 | Walter (1998) |
| Reacción | 0.25 | −0.89 | 1.23 | 0.41 | Walter (1998) |
Fracciones en masa; las reacciones, por unidad de masa de fundido. Un número negativo significa que ese mineral se forma mientras la roca funde. Coeficientes de reparto: olivino, O'Neill & Jenner (2012); opx y cpx, Yao et al. (2012); granate, Van Westrenen & Draper (2007); la espinela no admite REE (D = 0). Son las combinaciones de la Fig. 17c.
Ahora D aparece dentro de un logaritmo, sumado a F, y ya no hay atajo. Para cada F hay que calcular la función del proceso en las 13 REE y ajustarla con las piezas de la Tabla 2, como si fuera una muestra más: es la «fuerza bruta» que propone O'Neill. Al cambiar F, la punta del vector dibuja una curva:
D⁰ y P de la fuente
Función del proceso: −ln(D⁰ − F·P + F)
Trayectorias en el diagrama λ2–λ1
Tres cosas que se ven en la figura:
- Con poca fusión, el fundido se enriquece mucho en REE ligeras (ψ1 grande). Con F = 1 %, el La se multiplica por 41 a 51, según la fuente, y el Lu por 7.7 (espinela), por 2.5 (granate, 30 kbar) o por 0.9 (granate, 40 kbar).
- El granate hunde λ2: retiene las REE pesadas y deja una joroba. Cuanto más granate tiene la fuente (0.035 a 30 kbar, 0.14 a 40 kbar), más abajo va la trayectoria.
- Con mucha fusión, el fundido hereda la forma de la fuente. Cuando F es mucho mayor que los D, D⁰ − F·P + F se parece a F en todas las REE, y la función del proceso se acerca a −ln F, igual en todas. En la facies de espinela ocurre a partir de F ≈ 10 %: el vector se acerca a Ψ ≈ {−ln F, 0, 0} (pp. 1493–1494). Con granate, no: a 40 kbar y con F = 20 %, ψ2 todavía vale −43.
La prueba de fuerza bruta sale bien: con 5 términos, la función del proceso se ajusta en todos los casos con desviaciones de menos del 1.2 %, del orden de la precisión de un buen análisis. Con 3 términos llegan al 9 % en la fuente de 40 kbar.
Un resultado de O'Neill con ψ0. Si Ψ ≈ {−ln F, 0, 0}, entonces λ0 = λ0⁰ − ln F, y el grado de fusión sale de comparar el λ0 del magma parental con el de la fuente. Para el magma parental de los OFB, O'Neill obtiene λ0 = 2.36 ± 0.1, de la tendencia de λ0 frente al MgO a 9.5 % de MgO. Con las fuentes de la Tabla 6 (λ0⁰ = 0.47, 0.60 o 0.86), F = 0.15, 0.17 o 0.22. Su estimación es F ≈ 0.19, con un intervalo del 95 % de ±0.03 (p. 1494), de acuerdo con el modelo independiente de Niu & Batiza (1991) aplicado a 190 segmentos de dorsal (0.176 ± 0.020).
El método de la flecha. Así propone O'Neill usar el diagrama, con papel de calco «o su equivalente electrónico» (pp. 1490–1491). Supón una composición para la fuente y dibuja una flecha desde ella hasta la muestra. Traslada la flecha, sin girarla ni estirarla, hasta que su base quede en (0, 0). Si la punta cae cerca de una trayectoria, ya tienes un modelo de fusión y un F. Después, compruébalo con λ0: ψ0 debe ser la diferencia entre el λ0 de la muestra y el de la fuente. También se puede hacer al revés: si conoces F (por ejemplo, por λ0), la flecha desde (0, 0) hasta ese punto de la trayectoria, trasladada hasta la muestra, te dice dónde estaría su fuente. Pruébalo en la figura con el botón «Método de la flecha»: la fuente S y la muestra M son inventadas, y tienes que encontrar la fuente y el F que las unen.
Mezclas: curvas, no rectas
La mezcla de dos magmas suma concentraciones, no logaritmos. El La de una mezcla a partes iguales es la media del La de los dos magmas, pero su logaritmo no es la media de los logaritmos. Por eso las mezclas no dan rectas en el diagrama λ2–λ1. Además, en cada REE manda el magma más rico: con solo un 20 % del basalto de Hawái 8-2.70 (9.1 ppm de La) en VG 304 (1.42 ppm), el La de la mezcla se duplica (2.96 ppm), mientras que el Lu apenas cambia (0.432 frente a 0.472 ppm). La figura reproduce los dos ejemplos de la Fig. 17d de O'Neill:
Los patrones
Curvas de mezcla en el diagrama λ2–λ1
Aun así, las curvas de mezcla se pueden dibujar y son informativas; solo hay que calcularlas con las concentraciones, no interpolando entre los λ de los extremos. O'Neill también sitúa en la Fig. 17d la corteza continental superior y los sedimentos subducidos, porque son posibles contaminantes.
Practica
1. Con los coeficientes del granate de McKenzie & O'Nions (1991), δ1 = −39.7 y δ2 = 496, ¿cuánto vale ψ2 si cristaliza un 10 % de granate?
2. Un basalto ha acumulado un 20 % de cristales de olivino. ¿Qué coeficientes cambian respecto al líquido original?
3. Una fuente tiene λ0⁰ = 0.86 y el magma parental que sale de ella, λ0 = 2.36. Con la aproximación Ψ ≈ {−ln F, 0, 0}, ¿qué grado de fusión F le corresponde?
4. ¿Por qué la mezcla de dos magmas no da una recta en el diagrama λ2–λ1?
Resumen
- Un proceso multiplica cada REE por un factor. En logaritmos, suma la función del proceso, que también se describe con las piezas de la Tabla 2. Sus coeficientes son los ψ, y λ = λ⁰ + ψ, orden por orden.
- Los coeficientes de reparto dependen del radio de forma suave (deformación de la red) y se describen con las mismas piezas: D = δ0 + δ1f1 + δ2f2 + … (Tabla 7).
- Cristalización fraccionada: ψ0 = ln F·(δ0 − 1), ψ1 = ln F·δ1, ψ2 = ln F·δ2. Flechas rectas con dirección δ2/δ1. Solo el cpx, el anfíbol y el granate cambian la forma; el olivino y la plagioclasa solo cambian λ0, también cuando se acumulan.
- Fusión parcial: sin atajo; los ψ se obtienen ajustando la función del proceso para cada F, y dibujan trayectorias curvas. En la facies de espinela, con F > 10 %, Ψ ≈ {−ln F, 0, 0}: el fundido hereda la forma de la fuente y F ≈ eλ0⁰ − λ0. Para los OFB, F ≈ 0.19.
- Mezclas: curvas, no rectas, porque se suman concentraciones y no logaritmos.
El radio, de principio a fin
Esta era la pregunta con la que empezamos: dónde usa O'Neill los radios iónicos. Estas son las siete estaciones del mapa:
- 1Eje X. Cada y = ln([REE]/CI) se dibuja y se ajusta frente al radio del catión 3+ en coordinación 8 (Shannon, 1976). Lección 1
- 2Constantes de los polinomios ortogonales. β = 1.05477 es la media de los 13 radios; γ, δ y ε salen de exigir que las piezas no se pisen. Solo dependen de los radios. Lecciones 2 y 5
- 3Tabla de valores fn(ri). Cada radio sustituido en cada pieza. Lección 5
- 4El ajuste. La tabla del sistema solo contiene radios; con piezas ortogonales, λn = Σ y·fn / Σ fn². Lecciones 3 y 6
- 5La curva calculada. Evaluada en los 14 radios; en el del Eu (1.066 Å) da el Eu*. Lección 7
- 6Incertidumbres. ζn = 1/√Σ fn² convierte el error analítico en el error de cada λ. Lección 7
- 7Procesos. Los coeficientes de reparto dependen del radio y se describen con las mismas piezas fn(r): de ahí los vectores ψ de cristalización y de fusión. Lección 9
En las seis primeras, lo que aporta el radio es igual para todas las muestras: se calcula una sola vez. La séptima añade la física de los cristales, y también ella se reduce a las mismas piezas.
Fin del curso. Has recorrido todo el cálculo de O'Neill: de las ppm a los λ, sus errores, lo que significan y cómo se relacionan con los procesos. Para aplicarlo a tus análisis, vuelve a la Lección 8. Si quieres que revisemos juntos tus datos o cualquier punto del curso, pregúntame.
Fuentes de esta lección
- O'Neill (2016), pp. 1485–1494: ecuaciones 2 a 8, Tablas 5, 6 y 7 y Fig. 17 (vectores de cristalización, RTX, trayectorias de fusión y mezclas); estimación de FOFB, p. 1494.
- Datos citados por O'Neill en sus Tablas 5 y 7: McKenzie & O'Nions (1991), Niu (1997), Walter (1998, 2014), Tiepolo et al. (2000), Van Westrenen & Draper (2007), O'Neill & Jenner (2012), Sun & Liang (2012) y Yao et al. (2012). Fuentes de la Tabla 6: Salters & Stracke (2004), Workman & Hart (2005) y O'Neill & Palme (2008).
- Shaw, D. M. (2006). Trace Elements in Magmas: A Theoretical Treatment. Cambridge University Press (ecuaciones de fusión y cristalización).
- Wood, B. J. & Blundy, J. D. (2014). Trace element partitioning: the influences of ionic radius, cation charge, pressure, and temperature. Treatise on Geochemistry, 2.ª ed., vol. 3 (deformación de la red).
- Muestras de la mezcla: 8-2.70 (HSDP2, Feigenson et al., 2003) y VG 304 y VG 4080 (Jenner & O'Neill, 2012b), de las hojas de O'Neill.
- Cálculos propios: la comprobación de λ = λ⁰ + ψ con VG 304, los vectores de cristalización, las trayectorias de fusión (reproducen las curvas continuas de la Fig. 17c) y las curvas de mezcla (Fig. 17d), con las Tablas 1, 2, 5 y 7.