I have three data points which I performed a linear fit and obtained the 1 sigma uncertainty lines. Now I would like to generate 100k data point uniformly distributed between the 1 sigma error bars (the big triangle on the left side) but I do not have any idea how am I able to do that. Here is my code
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,0]**0.5, 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, fit_plus1sigma, fit_minus1sigma, facecolor="gray", alpha=0.15)
In NumPy there is a Numpy random triangle function, however, I was not able to implement that in my case and I am not even sure if that is the right approach. I appreciate any help.
Advertisement
Answer
You can use this answer by Severin Pappadeux to sample uniformly in a triangle shape. You need the corner points of your triangle for that.
To find where your lines intersect you can follow this answer by Norbu Tsering. Then you just need the top-left corner and bottom-left corner coordinates of your triangle shape.
Putting all of this together, you can solve your problem like this.
Find the intersection:
# 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, a2, b1, 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]), (xx[1], fit_minus1sigma[1])) cross = get_intersect(p1[0], p1[1], p2[0], p2[1])
This way you get your two points per line that are on it, and the intersection point, which you need to sample from within this triangle shape.
Then you can sample the points you need:
# Source: https://stackoverflow.com/a/47425047/5320601 def trisample(A, B, C): """ Given three vertices A, B, 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], p2[0], cross))
Example picture for 1000 points: