""" Example for Latin Hypercube Sampling
    100 samples are drawn using LHS and for standard MC
    The means and stdevs of x and y and the mean of x*y are closer to
    the theoretical means
    For MC you would need more samples to get comparable results
    To run this, you should have installed 
    numpy, scipy, pyDOE and matplotlib
"""
import matplotlib.pyplot as plt
import scipy.stats as ST
import numpy as NP
import pyDOE
 
NN = 100
points = pyDOE.lhs(2, samples=NN, criterion='center')
# points = pyDOE.lhs(2, samples=NN, criterion='maximin')
x = points[:, 0]
y = points[:, 1]
print("LHS")
print(NP.corrcoef(x, y))
print(NP.mean(x), NP.mean(y))  # should be 0.5, 0.5
print(NP.std(x), NP.std(y))  # should be 0.28867, 0.28867
print(NP.mean(x*y))  # should be 0.25
 
xa = ST.distributions.uniform.rvs(0, 1, size=NN)
ya = ST.distributions.uniform.rvs(0, 1, size=NN)
print("Standard")
print(NP.corrcoef(xa, ya))
print(NP.mean(xa), NP.mean(ya))
print(NP.std(xa), NP.std(ya))
print(NP.mean(xa*ya))
 
f, (xa1, xa2) = plt.subplots(1, 2, sharex=True, sharey=True)
xa1.scatter(x, y)
xa1.set_title("LHS")
 
xa2.scatter(xa, ya)
xa2.set_title("standard")
plt.show()

IiIiIEV4YW1wbGUgZm9yIExhdGluIEh5cGVyY3ViZSBTYW1wbGluZwogICAgMTAwIHNhbXBsZXMgYXJlIGRyYXduIHVzaW5nIExIUyBhbmQgZm9yIHN0YW5kYXJkIE1DCiAgICBUaGUgbWVhbnMgYW5kIHN0ZGV2cyBvZiB4IGFuZCB5IGFuZCB0aGUgbWVhbiBvZiB4KnkgYXJlIGNsb3NlciB0bwogICAgdGhlIHRoZW9yZXRpY2FsIG1lYW5zCiAgICBGb3IgTUMgeW91IHdvdWxkIG5lZWQgbW9yZSBzYW1wbGVzIHRvIGdldCBjb21wYXJhYmxlIHJlc3VsdHMKICAgIFRvIHJ1biB0aGlzLCB5b3Ugc2hvdWxkIGhhdmUgaW5zdGFsbGVkIAogICAgbnVtcHksIHNjaXB5LCBweURPRSBhbmQgbWF0cGxvdGxpYgoiIiIKaW1wb3J0IG1hdHBsb3RsaWIucHlwbG90IGFzIHBsdAppbXBvcnQgc2NpcHkuc3RhdHMgYXMgU1QKaW1wb3J0IG51bXB5IGFzIE5QCmltcG9ydCBweURPRQoKTk4gPSAxMDAKcG9pbnRzID0gcHlET0UubGhzKDIsIHNhbXBsZXM9Tk4sIGNyaXRlcmlvbj0nY2VudGVyJykKIyBwb2ludHMgPSBweURPRS5saHMoMiwgc2FtcGxlcz1OTiwgY3JpdGVyaW9uPSdtYXhpbWluJykKeCA9IHBvaW50c1s6LCAwXQp5ID0gcG9pbnRzWzosIDFdCnByaW50KCJMSFMiKQpwcmludChOUC5jb3JyY29lZih4LCB5KSkKcHJpbnQoTlAubWVhbih4KSwgTlAubWVhbih5KSkgICMgc2hvdWxkIGJlIDAuNSwgMC41CnByaW50KE5QLnN0ZCh4KSwgTlAuc3RkKHkpKSAgIyBzaG91bGQgYmUgMC4yODg2NywgMC4yODg2NwpwcmludChOUC5tZWFuKHgqeSkpICAjIHNob3VsZCBiZSAwLjI1Cgp4YSA9IFNULmRpc3RyaWJ1dGlvbnMudW5pZm9ybS5ydnMoMCwgMSwgc2l6ZT1OTikKeWEgPSBTVC5kaXN0cmlidXRpb25zLnVuaWZvcm0ucnZzKDAsIDEsIHNpemU9Tk4pCnByaW50KCJTdGFuZGFyZCIpCnByaW50KE5QLmNvcnJjb2VmKHhhLCB5YSkpCnByaW50KE5QLm1lYW4oeGEpLCBOUC5tZWFuKHlhKSkKcHJpbnQoTlAuc3RkKHhhKSwgTlAuc3RkKHlhKSkKcHJpbnQoTlAubWVhbih4YSp5YSkpCgpmLCAoeGExLCB4YTIpID0gcGx0LnN1YnBsb3RzKDEsIDIsIHNoYXJleD1UcnVlLCBzaGFyZXk9VHJ1ZSkKeGExLnNjYXR0ZXIoeCwgeSkKeGExLnNldF90aXRsZSgiTEhTIikKCnhhMi5zY2F0dGVyKHhhLCB5YSkKeGEyLnNldF90aXRsZSgic3RhbmRhcmQiKQpwbHQuc2hvdygp