CODEWITHSUNDEEP
pythonpython-3.xnumpymatplotlibrandom
Sara Krauss
Sara Krauss

Reputation: 398

Python uniform random number generation to a triangle shape

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)

enter image description here

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.

Upvotes: 2

Views: 443

Answers (1)

Lomtrur
Lomtrur

Reputation: 1798

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:

Example picture for 1000 points

Upvotes: 2

Related Questions