proper co-kriging formulation needs each variable to have its own drift columns (and its own unbiasedness row).
Here is what the current code actually builds vs. what the correct system is:
Current code — shared drift (one μ for all variables)
matsize = npp + ndrift + unbias
obs₁ obs₂ F 1
┌──────┬──────┬─────┬───┐
obs₁ │ C₁₁ │ C₁₂ │ F₁ │ 1 │ = c₀₁
obs₂ │ C₂₁ │ C₂₂ │ F₂ │ 1 │ = c₀₂
F^T │ F₁ᵀ │ F₂ᵀ │ 0 │ 0 │ = f₀ ← one shared μ
1^T │ 1ᵀ │ 1ᵀ │ 0 │ 0 │ = 1 ← one shared μ_u
└──────┴──────┴─────┴───┘
Both variables share one set of Lagrange multipliers. This is only correct if you assume the same trend parameters apply across all variables simultaneously — essentially "co-kriging with a common external drift."
Correct co-kriging — independent drift per variable
matsize = npp + nvar*(ndrift + unbias)
obs₁ obs₂ F₁ 1 F₂ 1
┌──────┬──────┬──────┬───┬──────┬───┐
obs₁ │ C₁₁ │ C₁₂ │ F₁ │ 1 │ 0 │ 0 │ = c₀₁
obs₂ │ C₂₁ │ C₂₂ │ 0 │ 0 │ F₂ │ 1 │ = c₀₂
F₁ᵀ │ F₁ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ = f₀ ← μ₁ for var 1
1₁ᵀ │ 1ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ = 1
F₂ᵀ │ 0 │ F₂ᵀ │ 0 │ 0 │ 0 │ 0 │ = f₀ ← μ₂ for var 2
1₂ᵀ │ 0 │ 1ᵀ │ 0 │ 0 │ 0 │ 0 │ = 1
└──────┴──────┴──────┴───┴──────┴───┘
Each variable has its own block-diagonal drift section and its own Lagrange multipliers. Variable 1's drift constraint involves only variable 1's observations; variable 2's constraint involves only variable 2's.
Also see Geostatistics Modeling Spatial Uncertainty
The generic form
← obs₁ → ← obs₂ → ←F₁ᵛ→ ←1₁→ ←F₂ᵛ→ ←1₂→ ←── Fˢ ──→
n₁ cols n₂ cols p_v 1 p_v 1 p_s cols
┌──────────┬──────────┬─────┬───┬─────┬───┬──────────┐
obs₁ n₁ │ C₁₁ │ C₁₂ │ F₁ᵛ │ 1 │ 0 │ 0 │ F₁ˢ │ = c₀₁
obs₂ n₂ │ C₂₁ │ C₂₂ │ 0 │ 0 │ F₂ᵛ │ 1 │ F₂ˢ │ = c₀₂
(F₁ᵛ)ᵀ p_v │ (F₁ᵛ)ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ 0 │ = f₀ᵛ
1₁ᵀ 1 │ 1ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ 0 │ = 1
(F₂ᵛ)ᵀ p_v │ 0 │ (F₂ᵛ)ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ = f₀ᵛ
1₂ᵀ 1 │ 0 │ 1ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ = 1
(Fˢ)ᵀ p_s │ (F₁ˢ)ᵀ │ (F₂ˢ)ᵀ │ 0 │ 0 │ 0 │ 0 │ 0 │ = f₀ˢ
└──────────┴──────────┴─────┴───┴─────┴───┴──────────┘
This can be implemented through the array obs%drift(ndrift, nobs) and block%drift(nvar,ndrift,nblock). Therefore, block%drift(ndrift,nblock) needs to be extended.
Practical implications
| Scenario |
Which model is correct |
| Single variable |
Either (shared = per-variable trivially) |
| Co-kriging, same physical trend drives all variables |
Shared may be acceptable as approximation |
| Co-kriging, variables have independent trends |
Must use per-variable |
| Ordinary co-kriging (unbias=1, ndrift=0) |
Per-variable gives proper Σλ_k=1 per variable |
proper co-kriging formulation needs each variable to have its own drift columns (and its own unbiasedness row).
Here is what the current code actually builds vs. what the correct system is:
Current code — shared drift (one μ for all variables)
matsize = npp + ndrift + unbiasBoth variables share one set of Lagrange multipliers. This is only correct if you assume the same trend parameters apply across all variables simultaneously — essentially "co-kriging with a common external drift."
Correct co-kriging — independent drift per variable
matsize = npp + nvar*(ndrift + unbias)Each variable has its own block-diagonal drift section and its own Lagrange multipliers. Variable 1's drift constraint involves only variable 1's observations; variable 2's constraint involves only variable 2's.
Also see Geostatistics Modeling Spatial Uncertainty
The generic form
This can be implemented through the array obs%drift(ndrift, nobs) and block%drift(
nvar,ndrift,nblock). Therefore, block%drift(ndrift,nblock) needs to be extended.Practical implications