Python uniform random number generation to a triangle shape

Question

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.

Note that you get a different solution if you sample x uniformly over the horizontal interval and y in the corresponding vertical interval, than when you want x and y uniformly inside the triangle. — JohanC, Aug 13 '20 at 14:19

Lomtrur · Accepted Answer · 2020-09-06T14:24:53.830

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:

@SeverinPappadeux Is it possible to give weight to the part of the triangle easily? For instance, lets say I want to have the data between x=335 and x=337 have double the weight (# of data points) of the rest. I thought about generating a new set of uniform data for that range and that add them up but I was wondering if there is a more pythonic way to it! — Sara Krauss, Aug 13 '20 at 15:37
@SaraKrauss Sure. General approach is to triangulate area (any area), then for uniform density you pick up triangle `i` first with probability A_i/A_total (`A` denotes area), then sample in the triangle. That will give you uniform points density. Then it is easy to modify density to any of your liking - you pick up some subset of triangles with higher probability (and rest of triangles will have lower probability) and then sample points in selected triangle as before — Severin Pappadeux, Aug 13 '20 at 15:53

Python uniform random number generation to a triangle shape

1 Answers1