我需要您的帮助来创建线性或更佳的截断拉普拉斯分布,该分布线性地依赖于另一个拉普拉斯分布变量。
不幸的是,这是到目前为止,我在派生的y分布中使用NaN实现的结果:
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde, truncnorm
slope = 0.2237
intercept = 1.066
spread = 4.8719
def dependency(x):
y_lin = slope * x + intercept
lower = slope / spread * 3 * x
upper = slope * spread / 3 * x + 2 * intercept
y_lin_noise = np.random.laplace(loc=0, scale=spread, size=len(y_lin)) + y_lin
y_lin_noise[y_lin_noise < lower] = np.nan # This is the desperate solution where
y_lin_noise[y_lin_noise > upper] = np.nan # NaNs are introduced
return y_lin_noise
max = 100
min = 1
mean = 40
sigma = 25
x = truncnorm((min-mean)/sigma, (max-mean)/sigma, loc=mean, scale=sigma).rvs(5000)
y = dependency(x)
# Plotting
xx = np.linspace(np.nanmin(x), np.nanmax(x), 100)
yy = slope * xx + intercept
lower = slope/spread*3*xx
upper = slope*spread/3*xx + 2*intercept
mask = ~np.isnan(y) & ~np.isnan(x)
x = x[mask]
y = y[mask]
xy = np.vstack([x, y])
z = gaussian_kde(xy)(xy)
idz = z.argsort()
x, y, z = x[idz], y[idz], z[idz]
fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(xx, upper, 'r-.', label='upper constraint')
plt.plot(xx, lower, 'r--', label='lower constraint')
ax.scatter(x, y, c=z, s=3)
plt.xlabel(r'$\bf X_{laplace}$')
plt.ylabel(r'$\bf Y_{{derived}}$')
plt.plot(xx, yy, 'r', label='regression model')
plt.legend()
plt.tight_layout()
plt.show()
最后,我要得到的是没有NaN的y分布,因此对于每个x,在上限/下限阈值范围内都有一个对应的y。可以这么说,回归线上的下/上截断分布。
我期待创意!
感谢您和最诚挚的问候。
您要做的基本上是重试超出范围的值。您可以定义一个针对每个值执行此操作的函数,也可以定义一个初始值并遍历答案以“更正”这样超出范围的值:
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kde, truncnorm
slope = 0.2237
intercept = 1.066
spread = 4.8719
#slow but effective
def truncated_noise(y, lower, upper):
while True:
y_noise = np.random.laplace(loc=0, scale=spread, size=1) + y
if upper > y_noise > lower:
return y_noise
def refine_noise(y, y_lin_noise, lower, upper):
for i in range(len(y_lin_noise)):
if upper[i] < y_lin_noise[i] or lower[i] > y_lin_noise[i]:
y_lin_noise[i] = truncated_noise(y[i], lower[i], upper[i])
return y_lin_noise
def dependency(x):
y_lin = slope * x + intercept
lower = slope / spread * 3 * x
upper = slope * spread / 3 * x + 2 * intercept
y_lin_noise = np.random.laplace(loc=0, scale=spread, size=len(y_lin)) + y_lin
y_lin_noise = refine_noise(y_lin, y_lin_noise, lower, upper)
return y_lin_noise
max = 100
min = 1
mean = 40
sigma = 25
x = truncnorm((min-mean)/sigma, (max-mean)/sigma, loc=mean, scale=sigma).rvs(5000)
y = dependency(x)
# Plotting
xx = np.linspace(np.nanmin(x), np.nanmax(x), 100)
yy = slope * xx + intercept
lower = slope/spread*3*xx
upper = slope*spread/3*xx + 2*intercept
xy = np.vstack([x, y])
z = gaussian_kde(xy)(xy)
idz = z.argsort()
x, y, z = x[idz], y[idz], z[idz]
fig, ax = plt.subplots(figsize=(5, 5))
plt.plot(xx, upper, 'r-.', label='upper constraint')
plt.plot(xx, lower, 'r--', label='lower constraint')
ax.scatter(x, y, c=z, s=3)
plt.xlabel(r'$\bf X_{laplace}$')
plt.ylabel(r'$\bf Y_{{derived}}$')
plt.plot(xx, yy, 'r', label='regression model')
plt.legend()
plt.tight_layout()
plt.show()
我不知道您想将其用作什么,但知道它不再是随机的拉普拉斯分布,因此,如果您的函数需要它,请不要这样做。