在保护生物学中，模拟物种的地理分布是一个重要的问题。本文以南美两种哺乳动物为例，基于过去的观察和14个环境变量来建模它们的地理分布。由于只有成功观察的案例（没有失败的观察），将这个问题视为一个密度估计问题，并使用OneClassSVM作为建模工具。数据集由Phillips等人（2006年）提供。如果可用，示例将使用basemap来绘制南美洲的海岸线和国界。
本文研究的两种物种是：
首先，使用OneClassSVM模型来拟合Bradypus variegatus的分布。模型训练完成后，绘制了南美洲的海岸线，并预测了物种的分布。ROC曲线下面积为0.868443，表明模型具有较好的区分能力。
接下来，对Microryzomys minutus进行了同样的建模过程。模型训练完成后，同样绘制了海岸线并预测了物种分布。ROC曲线下面积为0.993919，显示出极高的区分能力。
以下是使用Python和scikit-learn库实现上述建模过程的代码示例。代码中包含了数据加载、模型训练、预测分布以及绘制ROC曲线等步骤。
from time import time
import matplotlib.pyplot as plt
import numpy as np
from sklearn import metrics, svm
from sklearn.datasets import fetch_species_distributions
from sklearn.utils import Bunch
try:
from mpl_toolkits.basemap import Basemap
basemap = True
except ImportError:
basemap = False
def construct_grids(batch):
# 构建地图网格
xmin = batch.x_left_lower_corner + batch.grid_size
xmax = xmin + (batch.Nx * batch.grid_size)
ymin = batch.y_left_lower_corner + batch.grid_size
ymax = ymin + (batch.Ny * batch.grid_size)
xgrid = np.arange(xmin, xmax, batch.grid_size)
ygrid = np.arange(ymin, ymax, batch.grid_size)
return xgrid, ygrid
def create_species_bunch(species_name, train, test, coverages, xgrid, ygrid):
# 创建包含特定生物信息的bunch
bunch = Bunch(name=" ".join(species_name.split("_")[:2]))
species_name = species_name.encode("ascii")
points = dict(test=test, train=train)
for label, pts in points.items():
pts = pts[pts["species"] == species_name]
bunch["pts_%s" % label] = pts
ix = np.searchsorted(xgrid, pts["dd long"])
iy = np.searchsorted(ygrid, pts["dd lat"])
bunch["cov_%s" % label] = coverages[:, -iy, ix].T
return bunch
def plot_species_distribution(species=("bradypus_variegatus_0", "microryzomys_minutus_0")):
# 绘制物种分布图
t0 = time()
data = fetch_species_distributions()
xgrid, ygrid = construct_grids(data)
X, Y = np.meshgrid(xgrid, ygrid[::-1])
BV_bunch = create_species_bunch(species[0], data.train, data.test, data.coverages, xgrid, ygrid)
MM_bunch = create_species_bunch(species[1], data.train, data.test, data.coverages, xgrid, ygrid)
background_points = np.c_[np.random.randint(low=0, high=data.Ny, size=10000),
np.random.randint(low=0, high=data.Nx, size=10000), ].T
land_reference = data.coverages[6]
for i, species in enumerate([BV_bunch, MM_bunch]):
print("-" * 80)
print("Modeling distribution of species '%s'" % species.name)
mean = species.cov_train.mean(axis=0)
std = species.cov_train.std(axis=0)
train_cover_std = (species.cov_train - mean) / std
clf = svm.OneClassSVM(nu=0.1, kernel="rbf", gamma=0.5)
clf.fit(train_cover_std)
plt.subplot(1, 2, i + 1)
if basemap:
m = Basemap(projection="cyl", llcrnrlat=Y.min(), urcrnrlat=Y.max(), llcrnrlon=X.min(), urcrnrlon=X.max(), resolution="c")
m.drawcoastlines()
m.drawcountries()
else:
plt.contour(X, Y, land_reference, levels=[-9998], colors="k", linestyles="solid")
plt.xticks([])
plt.yticks([])
print("- predict species distribution")
Z = np.ones((data.Ny, data.Nx), dtype=np.float64)
idx = np.where(land_reference > -9999)
coverages_land = data.coverages[:, idx[0], idx[1]].T
pred = clf.decision_function((coverages_land - mean) / std)
Z *= pred.min()
Z[idx[0], idx[1]] = pred
levels = np.linspace(Z.min(), Z.max(), 25)
Z[land_reference == -9999] = -9999
plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)
plt.colorbar(format="%.2f")
plt.scatter(species.pts_train["dd long"], species.pts_train["dd lat"], s=2**2, c="black", marker="^", label="train")
plt.scatter(species.pts_test["dd long"], species.pts_test["dd lat"], s=2**2, c="black", marker="x", label="test")
plt.legend()
plt.title(species.name)
plt.axis("equal")
pred_background = Z[background_points[0], background_points[1]]
pred_test = clf.decision_function((species.cov_test - mean) / std)
scores = np.r_[pred_test, pred_background]
y = np.r_[np.ones(pred_test.shape), np.zeros(pred_background.shape)]
fpr, tpr, thresholds = metrics.roc_curve(y, scores)
roc_auc = metrics.auc(fpr, tpr)
plt.text(-35, -70, "AUC:%.3f" % roc_auc, ha="right")
print("\nArea under the ROC curve: %f" % roc_auc)
print("\ntime elapsed: %.2f s" % (time() - t0))
plot_species_distribution()
plt.show()