如何解决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)
在NumPy中有一个Numpy random triangle function,但是我无法在我的情况下实现,而且我什至不确定这是否是正确的方法。感谢您的帮助。
解决方法
您可以使用answer的Severin 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点:
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。