如何在随机体积中使用 LatinHyper Cube 采样来采样点

Question

我有一个由XYZ坐标组成的数据集，可以将其渲染为3D条形图，如下图所示。

我的主要问题是：如何在所述体积内生成 X 点（X 是用户输入）以确保尽可能的最佳表示？

我的第一个方法是使用 LatinHyper Cube 采样方法，但顾名思义，它假设一个立方设计空间。这里的体积可以是任何形状。

然后我尝试应用此处描述的方法：https://math.stackexchange.com/questions/2174751/generate-random-points-within-n-Dimension-ellipsoid我希望该方法可以调整为任何形状而不是椭球体。事实并非如此（或者至少，我没有这样做）。

作为一个“额外”问题，我对 3D 直方图方法不满意，我宁愿有一个“适当”的体积。使用

griddata

意味着立方底，但此处情况并非如此。

编辑：

我设法将数据绘制为 3D 表面，这非常适合我使用以下内容：

data_2d = [
    [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 0],
    [0, 0, 6, 8, 10, 12, 14, 16, 18, 20, 22],
    [0, 0, 0, 0, 0, 12, 14, 16, 18, 20, 22],
    [0, 0, 0, 0, 0, 12, 14, 16, 18, 20, 0],
    [0, 0, 0, 0, 0, 0, 14, 16, 18, 0, 0],
    [0, 0, 0, 0, 0, 0, 14, 16, 18, 0, 0],
    [0, 0, 0, 0, 0, 0, 14, 16, 18, 0, 0],
    [0, 0, 0, 0, 0, 0, 14, 16, 18, 0, 0],
]
# data_2d: - rows are Hs from 1 to 8 (8 rows)
#          - columns are Tp from 2 to 22 (10 columns)
#          - content is the wind speed from 2 to 22

data_array = np.array(data_2d)
x_data, y_data = np.meshgrid(np.linspace(2, 22, 11), np.linspace(1, 8, 8))


ax = plt.axes(projection="3d")
ax.plot_surface(x_data, y_data, data_array, cmap="viridis", edgecolor="black")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

Answer 1

实际上，我将问题分为两步：

1. 在以多边形为界的 XY 平面中生成 LHS 采样
2. 对于每个采样绘制一个由任意最小值和计算的上 Z 值界定的随机 Z 值。

以下是 MWE，供感兴趣的人参考：

    from random import uniform
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from scipy.interpolate import CloughTocher2DInterpolator as CT
    from scipy.stats import qmc
    from shapely.geometry import Point, Polygon

    data_2d = [
    [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, np.nan],
    [np.nan, np.nan, 6, 8, 10, 12, 14, 16, 18, 20, 22],
    [np.nan, np.nan, np.nan, np.nan, np.nan, 12, 14, 16, 18, 20, 22],
    [np.nan, np.nan, np.nan, np.nan, np.nan, 12, 14, 16, 18, 20, np.nan],
    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 14, 16, 18, np.nan, np.nan],
    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 14, 16, 18, np.nan, np.nan],
    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 14, 16, 18, np.nan, np.nan],
    [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 14, 16, 18, np.nan, np.nan],
]
# data_2d: - rows are Hs from 1 to 8 (8 rows)
#          - columns are Tp from 2 to 22 (10 columns)
#          - content is the wind speed from 2 to 22

tp_hs_ws = pd.DataFrame(data_2d)
tp_hs_ws.columns = [np.arange(2, 24, 2)]
tp_hs_ws.index = [np.arange(1, 9, 1)]
x_data, y_data = np.meshgrid(np.arange(2, 24, 2), np.arange(1, 9, 1))

non_nan_coord = [
    (2, 1),(20, 1),(22, 2),(22, 3),(22, 3),(20, 4),(18, 5),(18, 8),(14, 8),(14, 5),(12, 4),(12, 3),(10, 2),(6, 2),(2, 1)]

polygon = Polygon(non_nan_coord)
xp, yp = polygon.exterior.xy

points = LHS_Points_in_Polygon(polygon, nb_points)
xs = [point.x for point in points]
ys = [point.y for point in points]

# Keep only the unique LHS samples
xs = pd.Series(xs).unique()
ys = pd.Series(ys).unique()
xs_grid, ys_grid = np.meshgrid(xs, ys)

# Interpolate initial wind speed on the LHS Hs/Tp grid
zz = []
for z in (np.array(data_2d)).ravel():
    if str(z) == "nan":
        z = 0
    zz.append(z)

xy = np.c_[x_data.ravel(), y_data.ravel()]
CT_interpolant = CT(xy, zz)
Ws = CT_interpolant(xs_grid, ys_grid)

# Select the wind speed associated to the LHS Tp/Hs samples
ws = []
for idx_tp, _ in enumerate(xs_grid.ravel()):
    ws.append(Ws.ravel()[idx_tp])
# Make the LHS samples in square matrix form
ws_LHS = np.reshape(ws, (len(xs_grid), len(ys_grid)))
# The diagonal of wind speed LHS samples is corresponding to the XY coordinates sampled
ws_LHs_diag = ws_LHS.diagonal()

# Create random wind speed between 2m/s (arbitrary lower bound) and the LSH sampled wind speed value (upper bound)
# This ensure to produce a point XYZ always contained with the voume Tp/Hs/Wind speed
random_ws = [uniform(2, ws) for ws in ws_LHs_diag]

函数

LHS_Points_in_Polygon

受到这个解决方案的启发。

def LHS_Points_in_Polygon(polygon, number):
    minx, miny, maxx, maxy = polygon.bounds
    sampler = qmc.LatinHypercube(d=2, scramble=False)
    sample = sampler.random(n=number)
    l_bounds = np.min((minx, miny))
    u_bounds = np.max((maxx, maxy))
    points = []
    while len(points) < number:
        for x, y in qmc.scale(sample, l_bounds, u_bounds):
            pnt = Point(x, y)
            if polygon.contains(pnt):
                points.append(pnt)
    return points

结果如下：