如何基于二分类数据创建X1,X2,P(Y)的曲面/3D高斯图？

Got it, let's tackle this problem step by step. The issue you ran into makes total sense—surf/mesh functions need a continuous grid of Z values mapped to X1/X2, but your raw classification data is just discrete points. So we first need to fit a model that can predict P(Y) for any (X1,X2), then generate the grid to plot the surface. Below are working implementations in both Python and R, including overlaying your source/target dataset points.

Python Implementation (Using scikit-learn & Matplotlib) #

We'll use a flexible classification model (here, Gaussian Process Classifier for smooth surfaces, but you can swap in logistic regression, random forest, etc.) to predict P(Y) across a dense X1-X2 grid, then plot the surface with your raw data points.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF

# ----------------------
# Step 1: Simulate your data (replace this with your actual dataset)
# ----------------------
# Source dataset (red points)
np.random.seed(42)
source_X = np.random.normal(loc=[2, 2], scale=0.8, size=(50, 2))
source_Y = np.random.binomial(n=1, p=0.8, size=50)  # Higher P(Y=1) near [2,2]

# Target dataset (black points)
target_X = np.random.normal(loc=[-1, -1], scale=0.8, size=(50, 2))
target_Y = np.random.binomial(n=1, p=0.2, size=50)  # Lower P(Y=1) near [-1,-1]

# Combine all data for model fitting
X = np.vstack([source_X, target_X])
Y = np.hstack([source_Y, target_Y])

# ----------------------
# Step 2: Fit a model to predict P(Y)
# ----------------------
# Use Gaussian Process for smooth, non-linear surface (great for peak effects)
kernel = 1.0 * RBF(length_scale=1.0)
gpc = GaussianProcessClassifier(kernel=kernel, random_state=42)
gpc.fit(X, Y)

# ----------------------
# Step 3: Create a dense X1-X2 grid
# ----------------------
x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.linspace(x1_min, x1_max, 100),
                       np.linspace(x2_min, x2_max, 100))

# ----------------------
# Step 4: Predict P(Y=1) for every grid point
# ----------------------
Z = gpc.predict_proba(np.c_[xx1.ravel(), xx2.ravel()])[:, 1]
Z = Z.reshape(xx1.shape)

# ----------------------
# Step 5: Plot the surface + raw data points
# ----------------------
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# Plot the surface
ax.plot_surface(xx1, xx2, Z, cmap='viridis', alpha=0.7, edgecolor='none')

# Overlay source dataset (red circles)
ax.scatter(source_X[:, 0], source_X[:, 1], source_Y, color='red', s=50, label='Source Dataset')
# Overlay target dataset (black circles)
ax.scatter(target_X[:, 0], target_X[:, 1], target_Y, color='black', s=50, label='Target Dataset')

# Labels and formatting
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_zlabel('P(Y)')
ax.set_title('3D Surface Plot of P(Y) vs X1, X2')
ax.legend()

plt.show()

Notes: #

• If you want a simpler linear surface, replace GaussianProcessClassifier with LogisticRegression() from sklearn.linear_model.
• Adjust the length_scale parameter in the RBF kernel to control how "peaky" the surface is—smaller values create sharper peaks.

R Implementation (Using mgcv & plotly) #

We'll use a Generalized Additive Model (GAM) to fit a smooth non-linear relationship between X1/X2 and P(Y), then use plotly for an interactive 3D surface plot.

library(tidyverse)
library(mgcv)
library(plotly)

# ----------------------
# Step 1: Simulate your data (replace with your actual data)
# ----------------------
set.seed(42)
# Source dataset (red points)
source_data <- tibble(
  X1 = rnorm(50, mean = 2, sd = 0.8),
  X2 = rnorm(50, mean = 2, sd = 0.8),
  Y = rbinom(50, 1, 0.8),
  dataset = "Source"
)

# Target dataset (black points)
target_data <- tibble(
  X1 = rnorm(50, mean = -1, sd = 0.8),
  X2 = rnorm(50, mean = -1, sd = 0.8),
  Y = rbinom(50, 1, 0.2),
  dataset = "Target"
)

# Combine data
all_data <- bind_rows(source_data, target_data)

# ----------------------
# Step 2: Fit a GAM model to predict P(Y)
# ----------------------
# GAM fits smooth non-linear surfaces (great for peak effects)
gam_model <- gam(Y ~ s(X1) + s(X2), data = all_data, family = binomial)

# ----------------------
# Step 3: Create a dense X1-X2 grid
# ----------------------
grid <- expand.grid(
  X1 = seq(min(all_data$X1) - 1, max(all_data$X1) + 1, length.out = 100),
  X2 = seq(min(all_data$X2) - 1, max(all_data$X2) + 1, length.out = 100)
)

# ----------------------
# Step 4: Predict P(Y=1) for grid points
# ----------------------
grid$P_Y <- predict(gam_model, newdata = grid, type = "response")

# ----------------------
# Step 5: Plot interactive surface + raw points
# ----------------------
plot_ly() %>%
  # Add surface plot
  add_surface(
    x = ~grid$X1, y = ~grid$X2, z = ~grid$P_Y,
    colorscale = "Viridis", opacity = 0.7
  ) %>%
  # Add source dataset points (red)
  add_markers(
    data = source_data, x = ~X1, y = ~X2, z = ~Y,
    color = I("red"), size = I(10), name = "Source Dataset"
  ) %>%
  # Add target dataset points (black)
  add_markers(
    data = target_data, x = ~X1, y = ~X2, z = ~Y,
    color = I("black"), size = I(10), name = "Target Dataset"
  ) %>%
  # Format labels and title
  layout(
    scene = list(
      xaxis = list(title = "X1"),
      yaxis = list(title = "X2"),
      zaxis = list(title = "P(Y)")
    ),
    title = "3D Surface Plot of P(Y) vs X1, X2"
  )

Notes: #

• The s(X1) + s(X2) term in the GAM lets the model fit smooth, non-linear effects for each feature. For sharper peaks, adjust the k parameter in s() (e.g., s(X1, k=20) for more flexibility).
• Plotly gives you an interactive plot—you can zoom, rotate, and hover to see exact values.

内容的提问来源于stack exchange，提问作者MichaelMMeskhi