calcsd volume

import matplotlib.pyplot as plt
import numpy as np

# Average values are from old CalcSD site, not sure if it's still accurate
mean_length = 13.98
sd_length = 1.72

mean_girth = 11.64
sd_girth = 1.31

# Correlation between erect girth and bone-pressed erect length
corr = 0.55 #https://sci-hub.se/10.1111/jsm.12894

# Create Correlation Matrix
cor_mat = np.array([[1, corr], [corr, 1]])

# Create standard deviation matrix
std_mat = np.array([[sd_length**2, sd_length * sd_girth],
                    [sd_length * sd_girth, sd_girth**2]])

# Create covariance matrix
cov_mat = np.multiply(cor_mat, std_mat)

# Create mean matrix
avg_mat = np.array([[mean_length], [mean_girth]])

# Calculate
L = np.linalg.cholesky(cov_mat)
X = np.random.normal(size=(2, 1_000_000)) #calculate 1MIL values
Y = L.dot(X) + avg_mat

# Split into separate arrays
lengths, girths = Y[0, :], Y[1, :]

# Calculate volumes
volumes = 0.9 * lengths * np.pi * (girths / (2 * np.pi))**2

print(np.mean(volumes), np.std(volumes))

# Plot
volumes = np.sort(volumes)

fig, axs = plt.subplots(1)

axs.hist(volumes, bins = 1000)

plt.show()