Probability of event using Central Limit Theorem + plotting results
You seem to know about str.format (since you use it), yet you are still doing string addition:

title1 = 'Probability of making loss: ' + '{0:.4g}'.format(P1*100) + '%'
title2 = 'Probability of earning more than $20k: ' + '{0:.4g}'.format(P2*100) + '%'


Just make it a single string:

title1 = 'Probability of making loss: {0:.4g}%'.format(P1*100)
title2 = 'Probability of earning more than $20k: {0:.4g}%'.format(P2*100)


You could also use the % formatting option for floats (described in the fifth table here), which automatically multiplies with 100 and adds the % sign. The number of digits after the decimal point work differently, though, compared to the g format:

"{:.4%}".format(0.5)
# '50.0000%'
"{:.4g}%".format(0.5*100.)
# '50%'
"{:.4g}%".format(0.5123123*100.)
# '51.23%'
"{:.4%}".format(0.5123123)
# '51.2312%'




You should also re-organize your code according to the following scheme:

import something

GLOBAL_CONSTANT = "value"

class Definitions

def functions()

if __name__ == '__main__':
    main()
    # or some small code using the stuff defined above


That last part is there to protect your code from being executed if you want to import some part of this script from another script. At the moment if you want to do from laplace import Z, your whole code would run. With the if __name__ == '__main__' guard, it will only run when you execute this script.



You should add docstrings to your functions, describing what they do, what arguments they take and what they return. Have a look at PEP257 for the guidelines regarding docstrings.



And, in order to answer my own question from the comments:

Yes, it is justified to use the normal distribution here, since the Central Limit Theorem guarantees that the sum of random variables tends towards a normal distribution. The approximation is quite good for \$N = 100\$, as can be seen with the following small script:

import numpy as np
import matplotlib.pyplot as plt

mu = 95.7
sigma = 1247.
n = 100    # how many random variables to sum for each value
N = 1000    # how many values to generate

x_l = np.random.laplace(mu, sigma/np.sqrt(2), (n, N)).sum(axis=0)
x_g = np.random.normal(n*mu, np.sqrt(n)*sigma, N)

plt.hist(x_g, label="gaus")
plt.hist(x_l, label="sum(laplace)", histtype='step')
plt.legend()
plt.show()