Skip to content

Z-value added to cross validation #521

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
SimonMolinsky merged 1 commit into from
Mar 19, 2026
Merged

Z-value added to cross validation #521

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ Changes - from version >= 1.x
* [experimental] LSA-method Ordinary Kriging tests and experiments
* [enhancement] Users can set negative predictions to zero in Area-to-Area, Area-to-Point, and Centroid-based Poisson Kriging
* [docs] Added block-to-block distance equation, and corrected semivariance equation in functions docstrings
* [enhancement] Z-value added to cross validation

2025-12-26
----------
Expand Down
Binary file modified docs/build/doctrees/api/evaluate/evaluate.doctree
View file
Open in desktop
Binary file not shown.
Binary file modified docs/build/doctrees/environment.pickle
View file
Open in desktop
Binary file not shown.
135 changes: 0 additions & 135 deletions docs/build/html/_modules/pyinterpolate/distance/block.html
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -486,87 +486,6 @@ <h1>Source code for pyinterpolate.distance.block</h1><div class="highlight"><pre
<span class="p">)</span>


<span class="c1"># def _calc_b2b_dist_from_dataframe(</span>
<span class="c1"># ps_blocks: Union[pd.DataFrame, gpd.GeoDataFrame],</span>
<span class="c1"># lon_col_name: Union[str, Hashable],</span>
<span class="c1"># lat_col_name: Union[str, Hashable],</span>
<span class="c1"># val_col_name: Union[str, Hashable],</span>
<span class="c1"># block_id_col_name: Union[str, Hashable],</span>
<span class="c1"># verbose=False</span>
<span class="c1"># ) -&gt; pd.DataFrame:</span>
<span class="c1"># r&quot;&quot;&quot;</span>
<span class="c1"># Function calculates distances between the blocks&#39; point supports.</span>
<span class="c1">#</span>
<span class="c1"># Parameters</span>
<span class="c1"># ----------</span>
<span class="c1"># ps_blocks : Union[pd.DataFrame, gpd.GeoDataFrame]</span>
<span class="c1"># DataFrame with point supports and block indexes.</span>
<span class="c1">#</span>
<span class="c1"># lon_col_name : Union[str, Hashable]</span>
<span class="c1"># Longitude or x coordinate.</span>
<span class="c1">#</span>
<span class="c1"># lat_col_name : Union[str, Hashable]</span>
<span class="c1"># Latitude or y coordinate.</span>
<span class="c1">#</span>
<span class="c1"># val_col_name : Union[str, Hashable]</span>
<span class="c1"># The point support values column.</span>
<span class="c1">#</span>
<span class="c1"># block_id_col_name : Union[str, Hashable]</span>
<span class="c1"># Column with block names / indexes.</span>
<span class="c1">#</span>
<span class="c1"># verbose : bool, default = False</span>
<span class="c1"># Show progress bar.</span>
<span class="c1">#</span>
<span class="c1"># Returns</span>
<span class="c1"># -------</span>
<span class="c1"># block_distances : DataFrame</span>
<span class="c1"># Indexes and columns are block indexes, cells are distances.</span>
<span class="c1">#</span>
<span class="c1"># &quot;&quot;&quot;</span>
<span class="c1"># calculated_pairs = set()</span>
<span class="c1"># unique_blocks = list(ps_blocks[block_id_col_name].unique())</span>
<span class="c1">#</span>
<span class="c1"># col_set = [lon_col_name, lat_col_name, val_col_name]</span>
<span class="c1">#</span>
<span class="c1"># results = []</span>
<span class="c1">#</span>
<span class="c1"># for block_i in tqdm(unique_blocks, disable=not verbose):</span>
<span class="c1"># for block_j in unique_blocks:</span>
<span class="c1"># # Check if it was estimated</span>
<span class="c1"># if not (block_i, block_j) in calculated_pairs:</span>
<span class="c1"># if block_i == block_j:</span>
<span class="c1"># results.append([block_i, block_j, 0])</span>
<span class="c1"># else:</span>
<span class="c1"># i_value = ps_blocks[</span>
<span class="c1"># ps_blocks[block_id_col_name] == block_i</span>
<span class="c1"># ]</span>
<span class="c1"># j_value = ps_blocks[</span>
<span class="c1"># ps_blocks[block_id_col_name] == block_j</span>
<span class="c1"># ]</span>
<span class="c1"># value = _calculate_block_to_block_distance(</span>
<span class="c1"># i_value[col_set].to_numpy(),</span>
<span class="c1"># j_value[col_set].to_numpy()</span>
<span class="c1"># )</span>
<span class="c1"># results.append([block_i, block_j, value])</span>
<span class="c1"># results.append([block_j, block_i, value])</span>
<span class="c1"># calculated_pairs.add((block_i, block_j))</span>
<span class="c1"># calculated_pairs.add((block_j, block_i))</span>
<span class="c1">#</span>
<span class="c1"># # Create output dataframe</span>
<span class="c1"># df = pd.DataFrame(data=results, columns=[&#39;block_i&#39;, &#39;block_j&#39;, &#39;z&#39;])</span>
<span class="c1"># df = df.pivot_table(</span>
<span class="c1"># values=&#39;z&#39;,</span>
<span class="c1"># index=&#39;block_i&#39;,</span>
<span class="c1"># columns=&#39;block_j&#39;</span>
<span class="c1"># )</span>
<span class="c1">#</span>
<span class="c1"># # sort</span>
<span class="c1"># df = df.reindex(columns=unique_blocks)</span>
<span class="c1"># df = df.reindex(index=unique_blocks)</span>
<span class="c1">#</span>
<span class="c1"># return df</span>


<span class="c1"># noinspection PyUnresolvedReferences</span>
<span class="k">def</span> <span class="nf">_calc_b2b_dist_from_ps</span><span class="p">(</span><span class="n">ps_blocks</span><span class="p">:</span> <span class="s1">&#39;PointSupport&#39;</span><span class="p">)</span> <span class="o">-&gt;</span> <span class="n">Dict</span><span class="p">:</span>
<span class="w"> </span><span class="sa">r</span><span class="sd">&quot;&quot;&quot;</span>
Expand Down Expand Up @@ -739,60 +658,6 @@ <h1>Source code for pyinterpolate.distance.block</h1><div class="highlight"><pre



<span class="c1"># def _calculate_block_to_block_distance(ps_block_1: np.ndarray,</span>
<span class="c1"># ps_block_2: np.ndarray) -&gt; float:</span>
<span class="c1"># r&quot;&quot;&quot;</span>
<span class="c1"># Function calculates distance between two blocks&#39; point supports.</span>
<span class="c1">#</span>
<span class="c1"># Parameters</span>
<span class="c1"># ----------</span>
<span class="c1"># ps_block_1 : numpy array</span>
<span class="c1"># Point support of the first block.</span>
<span class="c1">#</span>
<span class="c1"># ps_block_2 : numpy array</span>
<span class="c1"># Point support of the second block.</span>
<span class="c1">#</span>
<span class="c1"># Returns</span>
<span class="c1"># -------</span>
<span class="c1"># weighted_distances : float</span>
<span class="c1"># Weighted point-support distance between blocks.</span>
<span class="c1">#</span>
<span class="c1"># Notes</span>
<span class="c1"># -----</span>
<span class="c1"># The weighted distance between blocks is derived from the equation given</span>
<span class="c1"># in publication [1] from References. This distance is weighted by</span>
<span class="c1">#</span>
<span class="c1"># References</span>
<span class="c1"># ----------</span>
<span class="c1"># .. [1] Goovaerts, P. Kriging and Semivariogram Deconvolution in the</span>
<span class="c1"># Presence of Irregular Geographical Units.</span>
<span class="c1"># Math Geosci 40, 101–128 (2008).</span>
<span class="c1"># https://doi.org/10.1007/s11004-007-9129-1</span>
<span class="c1">#</span>
<span class="c1"># TODO</span>
<span class="c1"># ----</span>
<span class="c1"># * Add reference equation to the special part of the documentation.</span>
<span class="c1"># &quot;&quot;&quot;</span>
<span class="c1">#</span>
<span class="c1"># a_shape = ps_block_1.shape[0]</span>
<span class="c1"># b_shape = ps_block_2.shape[0]</span>
<span class="c1"># ax = ps_block_1[:, 0].reshape(1, a_shape)</span>
<span class="c1"># bx = ps_block_2[:, 0].reshape(b_shape, 1)</span>
<span class="c1"># dx = ax - bx</span>
<span class="c1"># ay = ps_block_1[:, 1].reshape(1, a_shape)</span>
<span class="c1"># by = ps_block_2[:, 1].reshape(b_shape, 1)</span>
<span class="c1"># dy = ay - by</span>
<span class="c1"># aval = ps_block_1[:, -1].reshape(1, a_shape)</span>
<span class="c1"># bval = ps_block_2[:, -1].reshape(b_shape, 1)</span>
<span class="c1"># w = aval * bval</span>
<span class="c1">#</span>
<span class="c1"># dist = np.sqrt(dx ** 2 + dy ** 2, dtype=float, casting=&#39;unsafe&#39;)</span>
<span class="c1">#</span>
<span class="c1"># wdist = dist * w</span>
<span class="c1"># distances_sum = np.sum(wdist) / np.sum(w)</span>
<span class="c1"># return distances_sum</span>


<span class="k">def</span> <span class="nf">select_neighbors_in_range</span><span class="p">(</span><span class="n">data</span><span class="p">:</span> <span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">,</span>
<span class="n">current_lag</span><span class="p">:</span> <span class="nb">float</span><span class="p">,</span>
<span class="n">previous_lag</span><span class="p">:</span> <span class="nb">float</span><span class="p">):</span>
Expand Down
23 changes: 16 additions & 7 deletions docs/build/html/_modules/pyinterpolate/evaluate/cross_validation.html
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<title>pyinterpolate.evaluate.cross_validation &#8212; pyinterpolate 1.1.0 documentation</title>
<title>pyinterpolate.evaluate.cross_validation &#8212; pyinterpolate 1.2.0 documentation</title>



Expand Down Expand Up @@ -38,7 +38,7 @@
<link rel="preload" as="script" href="../../../_static/scripts/bootstrap.js?digest=8878045cc6db502f8baf" />
<link rel="preload" as="script" href="../../../_static/scripts/pydata-sphinx-theme.js?digest=8878045cc6db502f8baf" />

<script src="../../../_static/documentation_options.js?v=fc837d61"></script>
<script src="../../../_static/documentation_options.js?v=6efca38a"></script>
<script src="../../../_static/doctools.js?v=9bcbadda"></script>
<script src="../../../_static/sphinx_highlight.js?v=dc90522c"></script>
<script src="../../../_static/clipboard.min.js?v=a7894cd8"></script>
Expand Down Expand Up @@ -111,7 +111,7 @@



<p class="title logo__title">pyinterpolate 1.1.0 documentation</p>
<p class="title logo__title">pyinterpolate 1.2.0 documentation</p>

</a></div>

Expand Down Expand Up @@ -514,7 +514,7 @@ <h1>Source code for pyinterpolate.evaluate.cross_validation</h1><div class="high
<span class="sd"> * Mean Prediction Error,</span>
<span class="sd"> * Mean Kriging Error: ratio of variance of prediction errors to</span>
<span class="sd"> the average variance error of kriging,</span>
<span class="sd"> * array with: ``[coordinate x, coordinate y, prediction error, kriging estimate error]``</span>
<span class="sd"> * array with: ``[coordinate x, coordinate y, prediction error, kriging estimate error, z-value, z-ci-min, z-ci-max]``</span>

<span class="sd"> References</span>
<span class="sd"> ----------</span>
Expand Down Expand Up @@ -552,8 +552,6 @@ <h1>Source code for pyinterpolate.evaluate.cross_validation</h1><div class="high
<span class="sd"> &gt;&gt;&gt; print(validation_results[1]) # mean kriging error</span>
<span class="sd"> 1.6386630811210166</span>
<span class="sd"> &quot;&quot;&quot;</span>
<span class="c1"># TODO:</span>
<span class="c1"># Use (2) to calc Z-score</span>
<span class="c1"># TODO:</span>
<span class="c1"># Validation tutorials</span>
<span class="c1"># TODO:</span>
Expand Down Expand Up @@ -606,9 +604,20 @@ <h1>Source code for pyinterpolate.evaluate.cross_validation</h1><div class="high
<span class="n">preds</span> <span class="o">=</span> <span class="n">preds</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

<span class="n">prediction_error</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">preds</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
<span class="n">p1</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">preds</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>

<span class="n">z</span> <span class="o">=</span> <span class="n">prediction_error</span> <span class="o">/</span> <span class="p">(</span><span class="n">p1</span> <span class="o">*</span> <span class="n">preds</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
<span class="n">z_ci_min</span> <span class="o">=</span> <span class="n">z</span> <span class="o">-</span> <span class="mi">2</span><span class="o">*</span><span class="n">p1</span>
<span class="n">z_ci_max</span> <span class="o">=</span> <span class="n">z</span> <span class="o">+</span> <span class="mi">2</span><span class="o">*</span><span class="n">p1</span>

<span class="n">coordinates_and_errors</span><span class="o">.</span><span class="n">append</span><span class="p">(</span>
<span class="p">[</span><span class="n">preds</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">preds</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="n">prediction_error</span><span class="p">,</span> <span class="n">preds</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span>
<span class="p">[</span><span class="n">preds</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span>
<span class="n">preds</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span>
<span class="n">prediction_error</span><span class="p">,</span>
<span class="n">preds</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span>
<span class="n">z</span><span class="p">,</span>
<span class="n">z_ci_min</span><span class="p">,</span>
<span class="n">z_ci_max</span><span class="p">]</span>
<span class="p">)</span>

<span class="n">output_arr</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">coordinates_and_errors</span><span class="p">)</span>
Expand Down
2 changes: 1 addition & 1 deletion docs/build/html/api/evaluate/evaluate.html
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -516,7 +516,7 @@ <h2>Cross-validation<a class="headerlink" href="#cross-validation" title="Link t
<li><p>Mean Prediction Error,</p></li>
<li><p>Mean Kriging Error: ratio of variance of prediction errors to
the average variance error of kriging,</p></li>
<li><p>array with: <code class="docutils literal notranslate"><span class="pre">[coordinate</span> <span class="pre">x,</span> <span class="pre">coordinate</span> <span class="pre">y,</span> <span class="pre">prediction</span> <span class="pre">error,</span> <span class="pre">kriging</span> <span class="pre">estimate</span> <span class="pre">error]</span></code></p></li>
<li><p>array with: <code class="docutils literal notranslate"><span class="pre">[coordinate</span> <span class="pre">x,</span> <span class="pre">coordinate</span> <span class="pre">y,</span> <span class="pre">prediction</span> <span class="pre">error,</span> <span class="pre">kriging</span> <span class="pre">estimate</span> <span class="pre">error,</span> <span class="pre">z-value,</span> <span class="pre">z-ci-min,</span> <span class="pre">z-ci-max]</span></code></p></li>
</ul>
</div></blockquote>
</dd>
Expand Down
2 changes: 1 addition & 1 deletion docs/build/html/searchindex.js
View file
Open in desktop

Large diffs are not rendered by default.

17 changes: 13 additions & 4 deletions src/pyinterpolate/evaluate/cross_validation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ def validate_kriging(
* Mean Prediction Error,
* Mean Kriging Error: ratio of variance of prediction errors to
the average variance error of kriging,
* array with: ``[coordinate x, coordinate y, prediction error, kriging estimate error]``
* array with: ``[coordinate x, coordinate y, prediction error, kriging estimate error, z-value, z-ci-min, z-ci-max]``

References
----------
Expand Down Expand Up @@ -127,8 +127,6 @@ def validate_kriging(
1.6386630811210166
"""
# TODO:
# Use (2) to calc Z-score
# TODO:
# Validation tutorials
# TODO:
# Areal kriging validation
Expand Down Expand Up @@ -180,9 +178,20 @@ def validate_kriging(
preds = preds[0]

prediction_error = row[-1] - preds[0]
p1 = np.sqrt(preds[1])

z = prediction_error / (p1 * preds[0])
z_ci_min = z - 2*p1
z_ci_max = z + 2*p1

coordinates_and_errors.append(
[preds[2], preds[3], prediction_error, preds[1]]
[preds[2],
preds[3],
prediction_error,
preds[1],
z,
z_ci_min,
z_ci_max]
)

output_arr = np.array(coordinates_and_errors)
Expand Down
17 changes: 17 additions & 0 deletions tests/test_evaluate/test_cross_validation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
import math

import numpy as np

from pyinterpolate.semivariogram.experimental.classes.experimental_variogram import ExperimentalVariogram
Expand Down Expand Up @@ -76,3 +78,18 @@ def test_with_separate_geometry():

assert validation_results_sep[0] == validation_results[0]
assert validation_results_sep[1] == validation_results[1]


def test_z_scores():
validation_results = validate_kriging(
theoretical_model=THEORETICAL_MODEL,
points=ARMSTRONG_DATA,
no_neighbors=4,
progress_bar=False
)
validation_arr = validation_results[-1]

m = np.mean(validation_arr[:, 4])

assert math.isclose(m, 0, abs_tol=0.01)

Loading