Calculate the autism rate of Lesswrong + create histogram of AQ scores

lesswrong_aq.py

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# https://web.archive.org/web/20130116004605/http://raikoth.net/Stuff/LessWrong/for_public.csv
df = pd.read_csv('lesswrong-2012-survey.csv')

# Clean scores
autism_scores = pd.to_numeric(df['AutismScore'], errors='coerce').dropna().to_numpy()
autism_scores = autism_scores.astype(int)
autism_scores = [score for score in autism_scores if 1 <= score <= 50]

# Distribution parameters (https://docs.autismresearchcentre.com/papers/2001_BCetal_AQ.pdf)
normie_mean, normie_sd = 16.4, 6.3 # Group 2
autist_mean, autist_sd = 35.8, 6.5 # Group 1

def draw_histogram():
    plt.figure(figsize=(10, 6))

    # Histogram
    plt.hist(
        autism_scores,
        bins=np.arange(1, 52) - 0.5,  # center bins on integers
        density=True,
        alpha=0.7,
        label="LessWrong 2012 respondents"
    )

    # X values for smooth curves
    x = np.linspace(0, 50, 1000)

    # Gaussian curves
    normie_pdf = norm.pdf(x, loc=normie_mean, scale=normie_sd)
    autist_pdf = norm.pdf(x, loc=autist_mean, scale=autist_sd)

    plt.plot(
        x, normie_pdf,
        linestyle=":",
        linewidth=2,
        label=f"Non-autistic distribution (μ={normie_mean}, σ={normie_sd})"
    )
    plt.plot(
        x, autist_pdf,
        linestyle=":",
        linewidth=2,
        label=f"Autistic distribution (μ={autist_mean}, σ={autist_sd})"
    )

    plt.xlabel("AQ Score")
    plt.ylabel("Density")
    plt.title("Distribution of Autism Quotient (AQ) Scores")
    plt.xlim(0, 50)
    plt.xticks(range(0, 51, 5))
    plt.legend()
    plt.tight_layout()

    plt.savefig('images/lw-2012-aq-histogram.png', dpi=300)
    plt.show()

def print_autism_stats(prior_autist):
    prior_normie = 1 - prior_autist

    # Likelihoods under each distribution
    p_score_given_autist = norm.pdf(autism_scores, autist_mean, autist_sd)
    p_score_given_normie = norm.pdf(autism_scores, normie_mean, normie_sd)

    # Bayes theorem
    posterior_autist = (
        p_score_given_autist * prior_autist
    ) / (
        p_score_given_autist * prior_autist
        + p_score_given_normie * prior_normie
    )

    # Expected number of autistic people
    expected_autistic_count = posterior_autist.sum()

    print(f"Proportion of autists: {round(100*(expected_autistic_count/len(autism_scores)))}%")

for prior_odds in [0.001, 0.005, 0.01]:
    print(f"Prior Odds: {round(100*prior_odds, 1)}%")
    print_autism_stats(prior_odds)
    print()

draw_histogram()