This one requires geopandas but it's a quick and standard solution for sampling within odd shapes (called Monte Carlo Sampling ). Most of the comments below question outline the same concept.

Solution

# grab shape within which to sample
url = ""
us = gpd.read_file(url).explode()
## filter out parts of the US that are far away from mainland, I have no idea what they are (Guam islands?)
us = us.loc[us.geometry.apply(lambda x: x.exterior.bounds[2])<-60]

# grab bounding box within which to generate random numbers
x_min,y_min,x_max,y_max = us.geometry.unary_union.bounds 

# the sampling
np.random.seed(2) # set seed (needed for reproducible results
N = 10000
rndn_sample = pd.DataFrame({'x':np.random.uniform(x_min,x_max,N),'y':np.random.uniform(y_min,y_max,N)}) # actual generation
# re-save results in a geodataframe
rndn_sample = gpd.GeoDataFrame(rndn_sample, geometry = gpd.points_from_xy(x=rndn_sample.x, y=rndn_sample.y),crs = us.crs)

# filtering
inUS = rndn_sample['geometry'].apply(lambda s: s.within(us.geometry.unary_union)) # check if within the U.S. bounds
rndn_sample.loc[inUS,:].plot() # plot for visual inspection of results

Explanation

Grab US outline within which we want to randomly sample

# grab shapefile of the US from an official source
url = ""
us = gpd.read_file(url).explode() 

Note, with explode(), I expand the multi-part polygon into separate rows. This enables for easier filtering of the area we are interested in because we can grab bounds for each part of the multi-part polygon as below. Note that -60 is just an approximate longitude of the most eastern part of mainland US (Puerto Rico). Feel free to decrease it to exclude PR

## filter out parts of the US that are far away from mainland, I have no idea what they are (Guam islands?)
us = us.loc[us.geometry.apply(lambda x: x.exterior.bounds[2])<-60]