Python统一随机数生成为三角形

问题描述

我有三个数据点,我进行了线性拟合,并获得了1 sigma不确定性线。现在,我想生成均匀分布在1个sigma误差线(左侧的大三角形)之间的100k数据点,但是我不知道该怎么做。这是我的代码

import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.optimize import curve_fit

x = np.array([339.545772,339.545781,339.545803])
y = np.array([-0.430843,-0.43084,-0.430842])

def line(x,m,c):   
    return m*x + c

popt,pcov = curve_fit(line,x,y)
slope = popt[0]
intercept = popt[1]

xx = np.array([326.0,343.0])

fit  = line(xx,slope,intercept)
fit_plus1sigma = line(xx,slope + pcov[0,0]**0.5,intercept - pcov[1,1]**0.5)
fit_minus1sigma = line(xx,slope - pcov[0,intercept + pcov[1,1]**0.5)


plt.plot(xx,fit,"C4",label="Linear fit")
plt.plot(xx,fit_plus1sigma,'g--',label=r'One sigma uncertainty')
plt.plot(xx,fit_minus1sigma,'g--')
plt.fill_between(xx,facecolor="gray",alpha=0.15)

enter image description here

在NumPy中有一个Numpy random triangle function,但是我无法在我的情况下实现,而且我什至不确定这是否是正确的方法。感谢您的帮助。

解决方法

您可以使用answerSeverin Pappadeux来均匀采样为三角形。为此,您需要三角形的拐角点。

要找到线条相交的位置,可以按照answer跟随Norbu Tsering。然后,您只需要三角形形状的左上角和左下角坐标即可。

将所有这些放在一起,就可以解决这样的问题。

找到交叉点:

# Source: https://stackoverflow.com/a/42727584/5320601
def get_intersect(a1,a2,b1,b2):
    """
    Returns the point of intersection of the lines passing through a2,a1 and b2,b1.
    a1: [x,y] a point on the first line
    a2: [x,y] another point on the first line
    b1: [x,y] a point on the second line
    b2: [x,y] another point on the second line
    """
    s = np.vstack([a1,b2])  # s for stacked
    h = np.hstack((s,np.ones((4,1))))  # h for homogeneous
    l1 = np.cross(h[0],h[1])  # get first line
    l2 = np.cross(h[2],h[3])  # get second line
    x,y,z = np.cross(l1,l2)  # point of intersection
    if z == 0:  # lines are parallel
        return (float('inf'),float('inf'))
    return (x / z,y / z)

p1 = ((xx[0],fit_plus1sigma[0]),(xx[1],fit_plus1sigma[1]))
p2 = ((xx[0],fit_minus1sigma[0]),fit_minus1sigma[1]))
cross = get_intersect(p1[0],p1[1],p2[0],p2[1])

通过这种方式,您可以获得每条线上的两个点以及需要从该三角形形状内采样的交点。

然后您可以采样所需的点:

# Source: https://stackoverflow.com/a/47425047/5320601
def trisample(A,B,C):
    """
    Given three vertices A,C,sample point uniformly in the triangle
    """
    r1 = random.random()
    r2 = random.random()

    s1 = math.sqrt(r1)

    x = A[0] * (1.0 - s1) + B[0] * (1.0 - r2) * s1 + C[0] * r2 * s1
    y = A[1] * (1.0 - s1) + B[1] * (1.0 - r2) * s1 + C[1] * r2 * s1

    return (x,y)


points = []
for _ in range(100000):
    points.append(trisample(p1[0],cross))

示例图片为1000点:

Example picture for 1000 points