def pdf(s1, s2, mu1, mu2, sig1, sig2, theta):
    val1 = -(((s1-mu1)*np.cos(theta) - (s2-mu2)*np.sin(theta))**2)/(2*sig1**2)
    val2 = -(((s1-mu1)*np.sin(theta) + (s2-mu2)*np.cos(theta))**2)/(2*sig2**2)
    val3 = - np.log(sig1*sig2)
    return np.exp(val1 + val2 + val3)