rice_ilness_botec.js

function to(low, high) {
    if (low >= high) throw new Error("Upper bound must be greater than lower bound");
    
    const p = 0.9;
    const a = 0.147;
    const y = Math.log(1 - (2*p-1) ** 2);
    const z = 2/(Math.PI * a) + y/2;
    const zScore = Math.sqrt(2) * Math.sqrt(Math.sqrt(z*z - y/a) - z);
    
    if (low > 0) {
        const mu = (Math.log(low) + Math.log(high)) / 2;
        const sigma = (Math.log(high) - mu) / zScore;
        return () => Math.exp(mu + sigma * (
            Math.sqrt(-2 * Math.log(Math.random())) * 
            Math.cos(2 * Math.PI * Math.random())
        ));
    }
    const mu = (low + high) / 2;
    const sigma = (high - mu) / zScore;
    return () => mu + sigma * (
        Math.sqrt(-2 * Math.log(Math.random())) * 
        Math.cos(2 * Math.PI * Math.random())
    );
}

function simulateCase(region) {
    // Initial contamination
    const rawContaminationRate = region === 'US' ? to(0.4, 0.6)() : to(0.5, 0.7)();
    if (Math.random() > rawContaminationRate) return 0;
    
    // Spore survival through cooking
    const sporesSurvival = to(0.1, 0.3)();
    if (Math.random() > sporesSurvival) return 0;
    
    // Initial post-cooking load of surviving spores
    const initialLoad = to(10, 500)();  // CFU/g when present
    
    // Growth over 48h
    const logIncrease = region === 'US' ? to(2, 5)() : to(3, 6)();
    const finalLoad = initialLoad * Math.pow(10, logIncrease);
    
    // Illness threshold varies by strain and host
    const illnessThreshold = to(1e4, 1e6)();
    
    // Individual susceptibility
    const susceptibilityFactor = Math.random() < to(0.2, 0.4)();
    
    return finalLoad > illnessThreshold && susceptibilityFactor ? 1 : 0;
}

const trials = 50000;
let usIllnesses = 0;
let seaIllnesses = 0;

for (let i = 0; i < trials; i++) {
    usIllnesses += simulateCase('US');
    seaIllnesses += simulateCase('SEA');
}

console.log(`US illness probability: ${(usIllnesses/trials*100).toFixed(2)}% (80% CI: ${(usIllnesses/trials*80).toFixed(2)}%-${(usIllnesses/trials*120).toFixed(2)}%)`);
console.log(`SEA illness probability: ${(seaIllnesses/trials*100).toFixed(2)}% (80% CI: ${(seaIllnesses/trials*80).toFixed(2)}%-${(seaIllnesses/trials*120).toFixed(2)}%)`);

// Calculate component probabilities for one run
const sampleSize = 10000;
let usInitialContam = 0;
let usSporesSurvived = 0;
let usGrowthToThreshold = 0;

for (let i = 0; i < sampleSize; i++) {
    const rawContam = Math.random() < to(0.4, 0.6)();
    usInitialContam += rawContam;
    
    if (rawContam) {
        const sporesSurvived = Math.random() < to(0.1, 0.3)();
        usSporesSurvived += sporesSurvived;
    }
}

console.log(`\nComponent probabilities (US):
Raw contamination: ${(usInitialContam/sampleSize*100).toFixed(2)}%
Spores survival given contamination: ${(usSporesSurvived/usInitialContam*100).toFixed(2)}%`);