8  MCMC

import { createSlider, injectStyle } from "./_slider.js"
d3 = require("d3@7")

8.1 Thuật toán Metropolis - Hastings

viewof mh = {
  // ═══════════════════ MATH ═══════════════════
  function normalPDF(x,mu,s){return Math.exp(-0.5*((x-mu)/s)**2)/(s*Math.sqrt(2*Math.PI));}
  function lnGamma(z){
    if(z<0.5)return Math.log(Math.PI/Math.sin(Math.PI*z))-lnGamma(1-z);
    z--;const c=[0.99999999999980993,676.5203681218851,-1259.1392167224028,771.32342877765313,
      -176.61502916214059,12.507343278686905,-0.13857109526572012,9.9843695780195716e-6,1.5056327351493116e-7];
    let x=c[0];for(let i=1;i<9;i++)x+=c[i]/(z+i);
    const t=z+7.5;return 0.5*Math.log(2*Math.PI)+(z+0.5)*Math.log(t)-t+Math.log(x);
  }
  function lnBeta(a,b){return lnGamma(a)+lnGamma(b)-lnGamma(a+b);}
  function betaPDF(x,a,b){return(x<=0||x>=1)?0:Math.exp((a-1)*Math.log(x)+(b-1)*Math.log(1-x)-lnBeta(a,b));}
  function randn(){let u=0,v=0;while(!u)u=Math.random();while(!v)v=Math.random();return Math.sqrt(-2*Math.log(u))*Math.cos(2*Math.PI*v);}

  // ═══════════════════ DYNAMIC MODEL ═══════════════════
  let priorA=2,priorB=5,dataK=7,dataN=20;
  let rawLikMax=1,postA_=9,postB_=18;
  function recompute(){
    const mle=dataK/Math.max(dataN,1);
    rawLikMax=(mle>0&&mle<1)?Math.pow(mle,dataK)*Math.pow(1-mle,dataN-dataK):1e-20;
    postA_=priorA+dataK; postB_=priorB+dataN-dataK;
  }
  recompute();
  function prior(t){return betaPDF(t,priorA,priorB);}
  function likRaw(t){if(t<=0||t>=1)return 0;return Math.pow(t,dataK)*Math.pow(1-t,dataN-dataK);}
  function likNorm(t){return rawLikMax>0?likRaw(t)/rawLikMax:0;}
  function target(t){return prior(t)*likRaw(t);}
  function truePosterior(t){return betaPDF(t,postA_,postB_);}

  // ═══════════════════ PALETTE ═══════════════════
  const C={current:"#1565C0",propose:"#E65100",accept:"#2E7D32",reject:"#C62828",
    prior:"#7B1FA2",lik:"#00897B",post:"#E65100",
    hist:"#1565C0",histStroke:"#0D47A1",histTop:"#64B5F6",
    histNew:"#E65100",histNewStroke:"#BF360C",histNewTop:"#FF8A65",
    burned:"#EF9A9A",burnLine:"#E53935",
    filt:"#43A047",filtStroke:"#2E7D32",filtTop:"#81C784",
    ghost:"#CFD8DC",
    trace:"#1565C0",traceLight:"#90CAF9",proposal:"#FF6F00",
    txt:"#212121",sub:"#9E9E9E",border:"#E0E0E0"};
  const mono="'SF Mono',SFMono-Regular,Menlo,monospace";

  // ═══════════════════ STATE ═══════════════════
  let chain=[0.3],proposed=null,stepPhase="idle",accepted=null;
  let acceptRatio=0,uRand=0,samples=[],proposalSigma=0.08,stepCount=0;
  let autoRunning=false,autoTimer=null;
  let propStepIdx=0,lastAcceptedBin=-1;

  // ═══════════════════ WRAPPER ═══════════════════
  const outer=document.createElement("div");
  outer.style.cssText=`display:flex;flex-direction:column;align-items:center;gap:12px;
    font-family:'Inter',-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,sans-serif;
    width:100%;max-width:900px;margin:0 auto;`;
  if(!document.querySelector('link[href*="Inter"]')){
    const lk=document.createElement('link');lk.rel='stylesheet';
    lk.href='https://fonts.googleapis.com/css2?family=Inter:wght@400;500;600;700;800&display=swap';
    document.head.appendChild(lk);
  }
  outer.appendChild(injectStyle());

  // ── BUTTONS ──
  const ctrlRow=document.createElement("div");
  ctrlRow.style.cssText="display:flex;gap:10px;width:100%;align-items:center;flex-wrap:wrap;";
  function mkBtn(text,color,big){
    const b=document.createElement("button");
    b.style.cssText=`padding:${big?"10px 28px":"7px 16px"};border-radius:${big?12:10}px;border:none;
      background:linear-gradient(135deg,${color},${color}dd);color:#fff;font-size:${big?15:13}px;font-weight:700;cursor:pointer;
      font-family:inherit;transition:all 0.2s ease;letter-spacing:0.3px;
      ${big?"box-shadow:0 4px 14px "+color+"44;":"box-shadow:0 2px 8px rgba(0,0,0,0.12);"}`;
    b.textContent=text;
    b.addEventListener('mouseenter',()=>{b.style.transform='translateY(-1px)';b.style.boxShadow=`0 6px 20px ${color}55`;});
    b.addEventListener('mouseleave',()=>{b.style.transform='translateY(0)';b.style.boxShadow=big?`0 4px 14px ${color}44`:'0 2px 8px rgba(0,0,0,0.12)';});
    return b;
  }
  const btnStep=mkBtn("▶  Bước tiếp theo",C.current,true);
  const btnAuto=mkBtn("⏩ Chạy tự động","#6A1B9A",false);
  const btnReset=mkBtn("↺ Đặt lại","#78909C",false);
  const stepLabel=document.createElement("span");
  stepLabel.style.cssText=`font-size:13px;font-weight:700;color:${C.sub};font-family:${mono};margin-left:auto;`;
  ctrlRow.append(btnStep,btnAuto,btnReset,stepLabel);
  outer.appendChild(ctrlRow);

  // ── MODEL SLIDERS ──
  const SL={};
  SL.sigma=createSlider("Proposal σ",0.01,0.30,0.01,0.08,C.proposal,"amber");
  SL.pA=createSlider("Prior α",0.5,10,0.5,2,C.prior,"purple");
  SL.pB=createSlider("Prior β",0.5,10,0.5,5,C.prior,"purple");
  SL.dK=createSlider("Quan sát k",0,50,1,7,C.lik,"green");
  SL.dN=createSlider("Số thử n",1,50,1,20,C.lik,"green");

  const sr1=document.createElement("div");sr1.style.cssText="display:flex;gap:14px;width:100%;";
  sr1.append(SL.sigma.el,SL.pA.el,SL.pB.el);
  const sr2=document.createElement("div");sr2.style.cssText="display:flex;gap:14px;width:100%;";
  sr2.append(SL.dK.el,SL.dN.el);
  outer.append(sr1,sr2);

  // ── MCMC DIAGNOSTICS ──
  SL.burnin=createSlider("Burn-in",0,300,5,0,C.burnLine,"red");
  SL.thin=createSlider("Thinning (every k-th)",1,20,1,1,C.filt,"green");
  function mkToggle(label,color){
    const w=document.createElement("label");
    w.style.cssText=`display:inline-flex;align-items:center;gap:6px;cursor:pointer;
      padding:6px 14px;border-radius:10px;background:#f5f5f5;
      border:1.5px solid #E0E0E0;transition:all 0.2s;font-size:12px;font-weight:600;
      color:#9E9E9E;white-space:nowrap;user-select:none;`;
    const cb=document.createElement("input");cb.type="checkbox";
    cb.style.cssText=`accent-color:${color};width:15px;height:15px;cursor:pointer;`;
    w.append(cb,document.createTextNode(label));
    function s_(){w.style.background=cb.checked?color+'15':'#f5f5f5';w.style.borderColor=cb.checked?color:'#E0E0E0';w.style.color=cb.checked?color:'#9E9E9E';}
    cb.addEventListener("change",s_);
    return{el:w,cb};
  }
  const togBurn=mkToggle("Burn-in",C.burnLine);
  const togThin=mkToggle("Thinning",C.filt);
  const diagRow=document.createElement("div");
  diagRow.style.cssText="display:flex;gap:12px;width:100%;align-items:end;";
  diagRow.append(SL.burnin.el,SL.thin.el,togBurn.el,togThin.el);
  outer.appendChild(diagRow);

  // ── STATUS BOX ──
  const statusBox=document.createElement("div");
  statusBox.style.cssText=`width:100%;padding:14px 20px;border-radius:14px;
    background:linear-gradient(135deg,#FAFAFA,#F5F5F5);backdrop-filter:blur(8px);
    border:2px solid ${C.border};font-size:13px;font-family:${mono};color:${C.txt};
    line-height:1.7;min-height:44px;transition:all 0.3s ease;
    box-shadow:0 2px 12px rgba(0,0,0,0.04);`;

  // ═══════════════════ MAIN SVG ═══════════════════
  const W=900,totalH=260;
  const traceW=640,histW_=W-traceW;
  const mg={l:50,r:8,t:28,b:30};
  const tIW=traceW-mg.l-mg.r,tIH=totalH-mg.t-mg.b;
  const hMg={l:4,r:24,t:28,b:30};
  const hIW=histW_-hMg.l-hMg.r,hIH=totalH-hMg.t-hMg.b;

  const svg=d3.create("svg").attr("viewBox",[0,0,W,totalH])
    .style("width","100%").style("max-width",W+"px")
    .style("border-radius","14px").style("border",`1px solid ${C.border}`)
    .style("background","#fff").style("box-shadow","0 2px 16px rgba(0,0,0,0.06)");
  outer.appendChild(svg.node());

  // Trace panel
  const traceG=svg.append("g").attr("transform",`translate(${mg.l},${mg.t})`);
  traceG.append("rect").attr("x",-mg.l).attr("y",-mg.t).attr("width",traceW).attr("height",totalH).attr("fill","#FAFAFA");
  traceG.append("line").attr("x2",tIW).attr("y1",tIH).attr("y2",tIH).attr("stroke","#E0E0E0");
  traceG.append("line").attr("y2",tIH).attr("stroke","#E0E0E0");
  svg.append("text").attr("x",mg.l+tIW/2).attr("y",14).attr("text-anchor","middle")
    .attr("font-size",13).attr("font-weight",700).attr("fill",C.txt).text("Biểu đồ vết");
  svg.append("text").attr("x",14).attr("y",mg.t+tIH/2).attr("text-anchor","middle")
    .attr("font-size",11).attr("font-weight",600).attr("fill",C.sub)
    .attr("transform",`rotate(-90,14,${mg.t+tIH/2})`).text("θ");
  svg.append("text").attr("x",mg.l+tIW/2).attr("y",totalH-6).attr("text-anchor","middle")
    .attr("font-size",10).attr("fill",C.sub).attr("font-weight",600).text("bước");
  const traceBurnG=traceG.append("g");
  const traceLineG=traceG.append("g");
  const traceProposalG=traceG.append("g");
  const tracePointsG=traceG.append("g");

  // Histogram panel
  const histG=svg.append("g").attr("transform",`translate(${traceW+hMg.l},${hMg.t})`);
  histG.append("line").attr("y2",hIH).attr("stroke","#E0E0E0");
  const histTitle=svg.append("text").attr("x",traceW+histW_/2).attr("y",14).attr("text-anchor","middle")
    .attr("font-size",13).attr("font-weight",700).attr("fill",C.txt);
  svg.append("line").attr("x1",traceW).attr("x2",traceW).attr("y2",totalH).attr("stroke","#CFD8DC");
  const histGhostG=histG.append("g");
  const histBarsG=histG.append("g");
  const histCurveG=histG.append("g");

  // ── MERGED Prior+Likelihood ──
  const BH=200;
  const bsvg=d3.create("svg").attr("viewBox",[0,0,W,BH])
    .style("width","100%").style("max-width",W+"px")
    .style("border-radius","14px").style("border",`1px solid ${C.border}`)
    .style("background","#fff").style("box-shadow","0 2px 16px rgba(0,0,0,0.06)");
  outer.appendChild(bsvg.node());
  outer.appendChild(statusBox);

  const bMg={t:26,b:28,l:44,r:16};
  const bw=W-bMg.l-bMg.r,halfBH=(BH-bMg.t-bMg.b)/2;
  const bxS=d3.scaleLinear().domain([0,1]).range([0,bw]);
  const midY=bMg.t+halfBH;

  const mergedTitleEl=bsvg.append("text").attr("x",bMg.l+bw/2).attr("y",14).attr("text-anchor","middle")
    .attr("font-size",12).attr("font-weight",700).attr("fill",C.txt);
  const mergedG=bsvg.append("g").attr("transform",`translate(${bMg.l},0)`);
  mergedG.append("line").attr("x2",bw).attr("y1",midY).attr("y2",midY).attr("stroke","#BDBDBD");
  mergedG.append("line").attr("y1",bMg.t).attr("y2",BH-bMg.b).attr("stroke","#E0E0E0");
  bsvg.append("text").attr("x",bMg.l+bw/2).attr("y",BH-6).attr("text-anchor","middle")
    .attr("font-size",10).attr("fill",C.sub).attr("font-weight",600).text("θ");
  [0,0.25,0.5,0.75,1].forEach(v=>{
    mergedG.append("text").attr("x",bxS(v)).attr("y",BH-bMg.b+14).attr("text-anchor","middle")
      .attr("font-size",9).attr("fill",C.sub).attr("font-family",mono).text(v.toFixed(2));
  });
  bsvg.append("text").attr("x",12).attr("y",midY-halfBH/2+4)
    .attr("text-anchor","middle").attr("font-size",9).attr("font-weight",600).attr("fill",C.prior).text("Prior");
  bsvg.append("text").attr("x",12).attr("y",midY+halfBH/2+4)
    .attr("text-anchor","middle").attr("font-size",9).attr("font-weight",600).attr("fill",C.lik).text("Lik.");

  const staticG=mergedG.append("g");
  const mergedMarkers=mergedG.append("g");

  let pyS,lyS;
  const curvePts=d3.range(0,1.005,0.005);

  function drawStatic(){
    staticG.selectAll("*").remove();
    const pV=curvePts.map(t=>({x:t,y:prior(t)}));
    const lV=curvePts.map(t=>({x:t,y:likNorm(t)}));
    const mxP=Math.max(...pV.map(d=>d.y))||1;
    pyS=d3.scaleLinear().domain([0,mxP*1.1]).range([midY,bMg.t]);
    lyS=d3.scaleLinear().domain([0,1.1]).range([midY,BH-bMg.b]);
    const pArea=d3.area().x(d=>bxS(d.x)).y0(midY).y1(d=>pyS(d.y)).curve(d3.curveBasis);
    const pLine=d3.line().x(d=>bxS(d.x)).y(d=>pyS(d.y)).curve(d3.curveBasis);
    staticG.append("path").attr("d",pArea(pV)).attr("fill",C.prior).attr("opacity",0.13);
    staticG.append("path").attr("d",pLine(pV)).attr("fill","none").attr("stroke",C.prior).attr("stroke-width",2);
    const lArea=d3.area().x(d=>bxS(d.x)).y0(midY).y1(d=>lyS(d.y)).curve(d3.curveBasis);
    const lLine=d3.line().x(d=>bxS(d.x)).y(d=>lyS(d.y)).curve(d3.curveBasis);
    staticG.append("path").attr("d",lArea(lV)).attr("fill",C.lik).attr("opacity",0.13);
    staticG.append("path").attr("d",lLine(lV)).attr("fill","none").attr("stroke",C.lik).attr("stroke-width",2);
    mergedTitleEl.text(`Prior: Beta(${priorA},${priorB})  ·  Lik: Bin(${dataK}/${dataN},θ)  ·  Post: Beta(${postA_},${postB_})`);
  }
  drawStatic();

  // ═══════════════════ RENDER ═══════════════════
  function render(){
    proposalSigma=SL.sigma.val();
    const useBurn=togBurn.cb.checked,useThin=togThin.cb.checked;
    const burnin=SL.burnin.val(),thin=SL.thin.val();
    // Only actually filter when the toggle is on AND the value changes something
    const burnActive=useBurn&&burnin>0;
    const thinActive=useThin&&thin>1;
    const useFilter=burnActive||thinActive;
    const effBurn=burnActive?burnin:0,effThin=thinActive?thin:1;
    const curVal=chain[chain.length-1];
    const thetaS=d3.scaleLinear().domain([0,1]).range([tIH,0]);
    const histThetaS=d3.scaleLinear().domain([0,1]).range([hIH,0]);
    const maxVis=Math.max(chain.length+5,60);
    const stepS=d3.scaleLinear().domain([Math.max(0,chain.length-maxVis),chain.length+2]).range([0,tIW]);

    // ── TRACE ──
    traceLineG.selectAll("*").remove();
    traceProposalG.selectAll("*").remove();
    tracePointsG.selectAll("*").remove();
    traceBurnG.selectAll("*").remove();
    traceG.selectAll(".yt").remove();

    [0,0.25,0.5,0.75,1].forEach(v=>{
      traceG.append("text").attr("class","yt").attr("x",-8).attr("y",thetaS(v)+3)
        .attr("text-anchor","end").attr("font-size",9).attr("fill",C.sub).attr("font-family",mono).text(v.toFixed(2));
      traceG.append("line").attr("class","yt").attr("x1",0).attr("x2",tIW)
        .attr("y1",thetaS(v)).attr("y2",thetaS(v)).attr("stroke","#F0F0F0");
    });

    // Burn-in overlay
    if(burnActive&&burnin<chain.length){
      const bx=stepS(burnin);
      traceBurnG.append("rect").attr("x",0).attr("y",0).attr("width",Math.max(0,bx)).attr("height",tIH)
        .attr("fill",C.burnLine).attr("opacity",0.06);
      traceBurnG.append("line").attr("x1",bx).attr("x2",bx).attr("y1",0).attr("y2",tIH)
        .attr("stroke",C.burnLine).attr("stroke-width",2).attr("stroke-dasharray","6,4");
      traceBurnG.append("text").attr("x",bx+4).attr("y",12)
        .attr("font-size",10).attr("font-weight",700).attr("fill",C.burnLine).attr("font-family",mono).text(`burn-in = ${burnin}`);
    }

    // Trace line
    if(chain.length>1){
      if(burnActive&&burnin<chain.length){
        // Burned portion faded
        const bp=chain.slice(0,Math.min(burnin+1,chain.length));
        if(bp.length>1) traceLineG.append("path")
          .attr("d",bp.map((v,i)=>`${i?'L':'M'}${stepS(i).toFixed(1)},${thetaS(v).toFixed(1)}`).join(""))
          .attr("fill","none").attr("stroke",C.burned).attr("stroke-width",1.5).attr("opacity",0.4);
        const pp=chain.slice(burnin);
        if(pp.length>1) traceLineG.append("path")
          .attr("d",pp.map((v,i)=>`${i?'L':'M'}${stepS(burnin+i).toFixed(1)},${thetaS(v).toFixed(1)}`).join(""))
          .attr("fill","none").attr("stroke",C.trace).attr("stroke-width",1.3).attr("opacity",0.7);
      } else {
        traceLineG.append("path")
          .attr("d",chain.map((v,i)=>`${i?'L':'M'}${stepS(i).toFixed(1)},${thetaS(v).toFixed(1)}`).join(""))
          .attr("fill","none").attr("stroke",C.trace).attr("stroke-width",1.2).attr("opacity",0.6);
      }
    }

    // Points
    const visStart=Math.max(0,chain.length-Math.floor(maxVis));
    for(let i=visStart;i<chain.length;i++){
      const isLast=i===chain.length-1;
      const isBurned=burnActive&&i<burnin;
      const isKept=!isBurned&&(effThin<=1||(i-effBurn)%effThin===0);
      tracePointsG.append("circle")
        .attr("cx",stepS(i)).attr("cy",thetaS(chain[i]))
        .attr("r",isLast?5: (thinActive&&!isBurned&&isKept)?3:1.5)
        .attr("fill",isLast?C.current: isBurned?C.burned: isKept?C.trace:C.traceLight)
        .attr("opacity",isBurned?0.25: isKept?1:0.2)
        .attr("stroke",isLast?"#fff":"none").attr("stroke-width",isLast?2:0);
    }

    tracePointsG.append("text").attr("x",stepS(chain.length-1)).attr("y",thetaS(curVal)-10)
      .attr("text-anchor","middle").attr("font-size",11).attr("font-weight",700)
      .attr("fill",C.current).attr("font-family",mono).text(`θ = ${curVal.toFixed(3)}`);

    // Proposal
    if(stepPhase==="proposed"||stepPhase==="decided"){
      const pc=stepPhase==="decided"?(accepted?C.accept:C.reject):C.propose;
      const cIdx=propStepIdx-1,cT=chain[cIdx],cx=stepS(cIdx),px=stepS(propStepIdx);
      const nP=60,pp=[];
      for(let i=0;i<=nP;i++){const z=-3+6*i/nP;
        pp.push({theta:Math.max(0,Math.min(1,cT+z*proposalSigma)),density:normalPDF(z*proposalSigma,0,proposalSigma)});}
      const mD=Math.max(...pp.map(d=>d.density));
      let pD=`M${cx},${thetaS(pp[0].theta)}`;
      for(const d of pp) pD+=`L${(cx+(d.density/mD)*40).toFixed(1)},${thetaS(d.theta).toFixed(1)}`;
      pD+=`L${cx},${thetaS(pp[pp.length-1].theta)}Z`;
      traceProposalG.append("path").attr("d",pD).attr("fill",C.proposal).attr("opacity",0.18)
        .attr("stroke",C.proposal).attr("stroke-width",1).attr("stroke-opacity",0.3);
      traceProposalG.append("line").attr("x1",cx).attr("y1",thetaS(cT)).attr("x2",px).attr("y2",thetaS(proposed))
        .attr("stroke",pc).attr("stroke-width",1.5).attr("stroke-dasharray","4,3");
      traceProposalG.append("circle").attr("cx",px).attr("cy",thetaS(proposed)).attr("r",5)
        .attr("fill",pc).attr("stroke","#fff").attr("stroke-width",2);
      const sym=stepPhase==="decided"?(accepted?"✓":"✗"):"?";
      traceProposalG.append("text").attr("x",px).attr("y",thetaS(proposed)+(proposed>cT?-10:14))
        .attr("text-anchor","middle").attr("font-size",11).attr("font-weight",700)
        .attr("fill",pc).attr("font-family",mono).text(`${sym} θ*=${proposed.toFixed(3)}`);
    }

    // ── HISTOGRAM ──
    histGhostG.selectAll("*").remove();
    histBarsG.selectAll("*").remove();
    histCurveG.selectAll("*").remove();
    histG.selectAll(".yt").remove();

    const nBins=30,binH_=1/nBins,barH_=hIH/nBins;
    const binsRaw=Array(nBins).fill(0);
    samples.forEach(s=>{binsRaw[Math.min(Math.floor(s/binH_),nBins-1)]++;});

    // Filtered samples from full chain
    const filtSamp=[];
    for(let i=effBurn;i<chain.length;i++) if(effThin<=1||(i-effBurn)%effThin===0) filtSamp.push(chain[i]);
    const binsFilt=Array(nBins).fill(0);
    filtSamp.forEach(s=>{binsFilt[Math.min(Math.floor(s/binH_),nBins-1)]++;});

    const binsShow=useFilter?binsFilt:binsRaw;
    const maxBin=Math.max(...binsShow,1);
    const barSc=d3.scaleLinear().domain([0,maxBin]).range([0,hIW-8]);

    histTitle.text(useFilter?`Đã lọc (${filtSamp.length})`:"Biểu đồ tần suất");

    [0,0.25,0.5,0.75,1].forEach(v=>{
      histG.append("text").attr("class","yt").attr("x",hIW+6).attr("y",histThetaS(v)+3)
        .attr("text-anchor","start").attr("font-size",9).attr("fill",C.sub).attr("font-family",mono).text(v.toFixed(2));
    });

    // Ghost bars (raw) if filtering active
    if(useFilter){
      const mxR=Math.max(...binsRaw,1);
      const gSc=d3.scaleLinear().domain([0,mxR]).range([0,hIW-8]);
      binsRaw.forEach((c,bi)=>{if(!c)return;
        histGhostG.append("rect").attr("x",1).attr("y",histThetaS((bi+1)*binH_)+1)
          .attr("width",gSc(c)).attr("height",barH_-2).attr("rx",2).attr("fill",C.ghost).attr("opacity",0.4);
      });
    }

    // Active bars
    binsShow.forEach((count,bi)=>{
      if(!count)return;
      const by=histThetaS((bi+1)*binH_),bHt=barH_-2,tw=barSc(count);
      const isNew=bi===lastAcceptedBin&&!useFilter;
      const fc=useFilter?C.filt:isNew?C.histNew:C.hist;
      const sc=useFilter?C.filtStroke:isNew?C.histNewStroke:C.histStroke;
      const tc=useFilter?C.filtTop:isNew?C.histNewTop:C.histTop;
      histBarsG.append("rect").attr("x",1).attr("y",by+1).attr("width",tw).attr("height",bHt)
        .attr("rx",2).attr("fill",fc).attr("stroke",sc).attr("stroke-width",0.6);
      if(bHt>3) histBarsG.append("rect").attr("x",1).attr("y",by+1).attr("width",tw)
        .attr("height",Math.min(2,bHt*0.25)).attr("rx",2).attr("fill",tc).attr("opacity",0.5);
      const cw=Math.min(10,tw/count);
      if(cw>3) for(let j=1;j<count;j++){const lx=j*cw+1;
        if(lx<tw) histBarsG.append("line").attr("x1",lx).attr("x2",lx).attr("y1",by+1).attr("y2",by+1+bHt)
          .attr("stroke",sc).attr("stroke-width",0.4).attr("opacity",0.5);}
    });

    // True posterior curve
    const trueV=curvePts.map(t=>({theta:t,density:truePosterior(t)}));
    const mxT=Math.max(...trueV.map(d=>d.density));
    histCurveG.append("path")
      .attr("d",d3.line().x(d=>d.density*((hIW-8)/mxT*0.85)).y(d=>histThetaS(d.theta)).curve(d3.curveBasis)(trueV))
      .attr("fill","none").attr("stroke",C.post).attr("stroke-width",2).attr("stroke-dasharray","5,3").attr("opacity",0.6);

    // ── MERGED MARKERS ──
    mergedMarkers.selectAll("*").remove();
    function drawMM(val,color,label){
      const px=bxS(val),pP=pyS(prior(val)),pL=lyS(likNorm(val));
      mergedMarkers.append("line").attr("x1",px).attr("x2",px).attr("y1",pP).attr("y2",pL)
        .attr("stroke",color).attr("stroke-width",1.5).attr("stroke-dasharray","3,2");
      mergedMarkers.append("circle").attr("cx",px).attr("cy",pP).attr("r",4).attr("fill",color).attr("stroke","#fff").attr("stroke-width",1.5);
      mergedMarkers.append("circle").attr("cx",px).attr("cy",pL).attr("r",4).attr("fill",color).attr("stroke","#fff").attr("stroke-width",1.5);
      mergedMarkers.append("text").attr("x",px).attr("y",Math.max(bMg.t+6,pP-8)).attr("text-anchor","middle")
        .attr("font-size",9).attr("font-weight",700).attr("fill",color).attr("font-family",mono).text(`π=${prior(val).toFixed(2)}`);
      mergedMarkers.append("text").attr("x",px).attr("y",Math.min(BH-bMg.b-2,pL+14)).attr("text-anchor","middle")
        .attr("font-size",9).attr("font-weight",700).attr("fill",color).attr("font-family",mono).text(`L=${likNorm(val).toFixed(3)}`);
      mergedMarkers.append("text").attr("x",px).attr("y",Math.max(bMg.t+6,pP-18)).attr("text-anchor","middle")
        .attr("font-size",8).attr("font-weight",600).attr("fill",color).attr("font-family",mono).text(label);
    }
    drawMM(curVal,C.current,`θ=${curVal.toFixed(3)}`);
    if(proposed!==null&&(stepPhase==="proposed"||stepPhase==="decided"))
      drawMM(proposed,stepPhase==="decided"?(accepted?C.accept:C.reject):C.propose,`θ*=${proposed.toFixed(3)}`);

    // ── STATUS ──
    const nAcc=samples.length,rate=stepCount>0?(nAcc/stepCount*100).toFixed(0):"—";
    const pill=(t,bg)=>`<span style="display:inline-block;padding:2px 10px;border-radius:20px;background:${bg};color:#fff;font-weight:700;font-size:12px;letter-spacing:0.3px;">${t}</span>`;
    const vl=(v,c,big)=>`<span style="color:${c};font-weight:${big?800:700};font-size:${big?16:13}px;">${v}</span>`;
    const lb=(t,c)=>`<span style="color:${c||C.sub};font-size:11px;font-weight:600;">${t}</span>`;
    const stats=`<div style="display:flex;gap:16px;align-items:center;margin-top:6px;font-size:11px;color:${C.sub};">`+
      `${lb('Bước')} <b>${stepCount}</b> &nbsp;${lb('Đã nhận')} <b>${nAcc}</b> &nbsp;${lb('Tỷ lệ')} <b>${rate}%</b>`+
      (useFilter?` &nbsp;${lb('Đã lọc',C.filt)} <b style="color:${C.filt}">${filtSamp.length}</b>`:"")+`</div>`;

    if(stepPhase==="idle"){
      statusBox.style.borderColor=C.border;statusBox.style.background='linear-gradient(135deg,#FAFAFA,#F5F5F5)';
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:12px;">`+pill('SẴN SÀNG','#546E7A')+
        `<span style="font-size:14px;">Vị trí hiện tại: ${vl('θ = '+curVal.toFixed(4),C.current,true)}</span>`+
        `<span style="margin-left:auto;color:${C.sub};font-size:12px;">Nhấn <b style="color:${C.current}">Bước tiếp theo</b> ▶</span></div>`+stats;
    } else if(stepPhase==="proposed"){
      const pC=prior(curVal),pP=prior(proposed),lC=likNorm(curVal),lP=likNorm(proposed);
      statusBox.style.borderColor=C.propose;statusBox.style.background='linear-gradient(135deg,#FFF8E1,#FFF3E0)';
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:10px;margin-bottom:8px;">`+
        pill('ĐỀ XUẤT',C.propose)+`<span style="font-size:13px;">θ* = ${vl(proposed.toFixed(4),C.propose,true)} từ N(${curVal.toFixed(3)}, ${proposalSigma}²)</span></div>`+
        `<table style="width:100%;border-collapse:separate;border-spacing:0;border-radius:10px;overflow:hidden;border:1px solid #E0E0E0;font-size:12px;">`+
        `<thead><tr style="background:#F5F5F5;"><th style="padding:6px 12px;text-align:left;font-weight:600;color:${C.sub};border-bottom:1px solid #E0E0E0;"></th>`+
        `<th style="padding:6px 12px;text-align:center;font-weight:700;color:${C.current};border-bottom:1px solid #E0E0E0;">Hiện tại θ = ${curVal.toFixed(3)}</th>`+
        `<th style="padding:6px 12px;text-align:center;font-weight:700;color:${C.propose};border-bottom:1px solid #E0E0E0;">Đề xuất θ* = ${proposed.toFixed(3)}</th></tr></thead><tbody>`+
        `<tr style="background:#fff;"><td style="padding:5px 12px;font-weight:600;color:${C.prior};border-bottom:1px solid #F0F0F0;">π(·) Prior</td>`+
          `<td style="padding:5px 12px;text-align:center;border-bottom:1px solid #F0F0F0;">${vl(pC.toFixed(4),C.prior)}</td>`+
          `<td style="padding:5px 12px;text-align:center;border-bottom:1px solid #F0F0F0;">${vl(pP.toFixed(4),C.prior)}</td></tr>`+
        `<tr style="background:#FAFAFA;"><td style="padding:5px 12px;font-weight:600;color:${C.lik};border-bottom:1px solid #F0F0F0;">L(·) Likelihood</td>`+
          `<td style="padding:5px 12px;text-align:center;border-bottom:1px solid #F0F0F0;">${vl(lC.toFixed(4),C.lik)}</td>`+
          `<td style="padding:5px 12px;text-align:center;border-bottom:1px solid #F0F0F0;">${vl(lP.toFixed(4),C.lik)}</td></tr>`+
        `<tr style="background:#fff;"><td style="padding:5px 12px;font-weight:600;color:${C.txt};">π·L Mục tiêu</td>`+
          `<td style="padding:5px 12px;text-align:center;">${vl((pC*lC).toFixed(6),C.txt)}</td>`+
          `<td style="padding:5px 12px;text-align:center;">${vl((pP*lP).toFixed(6),C.txt)}</td></tr></tbody></table>`+
        `<div style="display:flex;align-items:center;gap:20px;margin-top:8px;padding:8px 12px;background:#fff;border-radius:8px;border:1px solid #E0E0E0;">`+
          `<span>${lb('Tỷ lệ chấp nhận','#333')} &nbsp;α = ${vl(acceptRatio.toFixed(4),'#333',true)}</span>`+
          `<span style="color:#BDBDBD;">│</span>`+
          `<span>${lb('Rút ngẫu nhiên','#333')} &nbsp;u = ${vl(uRand.toFixed(4),'#E65100',true)}</span>`+
          `<span style="margin-left:auto;font-size:12px;color:${C.sub};">Nhấn <b>Bước tiếp theo</b> để quyết định</span></div>`;
    } else {
      const dc=accepted?C.accept:C.reject;
      statusBox.style.borderColor=dc;
      statusBox.style.background=accepted?'linear-gradient(135deg,#E8F5E9,#F1F8E9)':'linear-gradient(135deg,#FFEBEE,#FBE9E7)';
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:12px;">`+
        pill(accepted?'✓ CHẤP NHẬN':'✗ TỪ CHỐI',dc)+
        `<span style="font-size:14px;font-weight:600;">α = ${vl(acceptRatio.toFixed(4),C.txt)} `+
        `<span style="color:${dc};font-weight:800;font-size:15px;">${accepted?'≥':'<'}</span> `+
        `u = ${vl(uRand.toFixed(4),C.propose)}</span>`+
        `<span style="margin-left:auto;font-size:13px;color:${dc};font-weight:600;">`+
        (accepted?`→ Di chuyển đến θ* = ${proposed.toFixed(4)}`:`→ Giữ tại θ = ${curVal.toFixed(4)}`)+`</span></div>`+stats;
    }
    stepLabel.textContent=`Bước ${stepCount} · ${nAcc} mẫu`;
  }

  // ═══════════════════ STEP LOGIC ═══════════════════
  function resetChain(){chain=[0.3];proposed=null;stepPhase="idle";accepted=null;samples=[];stepCount=0;propStepIdx=0;lastAcceptedBin=-1;}
  function doStep(){
    proposalSigma=SL.sigma.val();
    if(stepPhase==="idle"){
      const cur=chain[chain.length-1];
      let prop=cur+randn()*proposalSigma;
      prop=Math.max(0.001,Math.min(0.999,prop));
      proposed=prop;propStepIdx=chain.length;
      const tC=target(cur),tP=target(prop);
      acceptRatio=tC>0?Math.min(1,tP/tC):1;
      uRand=Math.random();stepPhase="proposed";render();
    } else if(stepPhase==="proposed"){
      accepted=uRand<=acceptRatio;stepCount++;
      stepPhase="decided";render();
    } else {
      if(accepted){chain.push(proposed);samples.push(proposed);lastAcceptedBin=Math.min(Math.floor(proposed/(1/30)),29);}
      else{chain.push(chain[chain.length-1]);lastAcceptedBin=-1;}
      proposed=null;stepPhase="idle";render();
      if(autoRunning)setTimeout(doStep,50);
    }
  }
  function fullStep(){
    proposalSigma=SL.sigma.val();
    const cur=chain[chain.length-1];
    let prop=cur+randn()*proposalSigma;
    prop=Math.max(0.001,Math.min(0.999,prop));
    propStepIdx=chain.length;
    const tC=target(cur),tP=target(prop);
    acceptRatio=tC>0?Math.min(1,tP/tC):1;
    uRand=Math.random();accepted=uRand<=acceptRatio;stepCount++;
    if(accepted){chain.push(prop);samples.push(prop);lastAcceptedBin=Math.min(Math.floor(prop/(1/30)),29);}
    else{chain.push(cur);lastAcceptedBin=-1;}
    if(chain.length>500)chain=chain.slice(-400);
    proposed=prop;stepPhase="decided";render();stepPhase="idle";
  }

  // ═══════════════════ EVENTS ═══════════════════
  function stopAuto(){autoRunning=false;clearInterval(autoTimer);btnAuto.textContent="⏩ Chạy tự động";btnAuto.style.background="linear-gradient(135deg,#6A1B9A,#6A1B9Add)";}
  btnStep.addEventListener("click",()=>{if(autoRunning)stopAuto();doStep();});
  btnAuto.addEventListener("click",()=>{
    if(autoRunning)stopAuto();
    else{autoRunning=true;btnAuto.textContent="⏸ Tạm dừng";btnAuto.style.background=`linear-gradient(135deg,${C.reject},${C.reject}dd)`;
      stepPhase="idle";autoTimer=setInterval(()=>{if(autoRunning)fullStep();else clearInterval(autoTimer);},80);}
  });
  btnReset.addEventListener("click",()=>{stopAuto();resetChain();render();});
  SL.sigma.input.addEventListener("input",()=>SL.sigma.sync());

  // Model changes → rebuild curves + reset chain
  function onModel(){
    for(const s of Object.values(SL))s.sync();
    const nA=SL.pA.val(),nB=SL.pB.val(),nK=SL.dK.val(),nN=SL.dN.val();
    if(nK>nN){SL.dK.input.value=nN;SL.dK.sync();}
    const dk=Math.min(SL.dK.val(),nN);
    if(nA!==priorA||nB!==priorB||dk!==dataK||nN!==dataN){
      priorA=nA;priorB=nB;dataK=dk;dataN=nN;
      recompute();drawStatic();resetChain();
    }
    render();
  }
  [SL.pA,SL.pB,SL.dK,SL.dN].forEach(s=>s.input.addEventListener("input",onModel));

  // Diagnostics → just re-render (no reset)
  SL.burnin.input.addEventListener("input",()=>{SL.burnin.sync();render();});
  SL.thin.input.addEventListener("input",()=>{SL.thin.sync();render();});
  togBurn.cb.addEventListener("change",render);
  togThin.cb.addEventListener("change",render);

  render();
  invalidation.then(()=>{autoRunning=false;clearInterval(autoTimer);});
  outer.value={};return outer;
}

8.2 Fit mô hình

viewof sir = {
  // ═══════════════════ MATH UTILITIES ═══════════════════
  function lnGamma(z){
    if(z<0.5)return Math.log(Math.PI/Math.sin(Math.PI*z))-lnGamma(1-z);
    z--;const c=[0.99999999999980993,676.5203681218851,-1259.1392167224028,771.32342877765313,
      -176.61502916214059,12.507343278686905,-0.13857109526572012,9.9843695780195716e-6,1.5056327351493116e-7];
    let x=c[0];for(let i=1;i<9;i++)x+=c[i]/(z+i);
    const t=z+7.5;return 0.5*Math.log(2*Math.PI)+(z+0.5)*Math.log(t)-t+Math.log(x);
  }
  function randn(){let u=0,v=0;while(!u)u=Math.random();while(!v)v=Math.random();return Math.sqrt(-2*Math.log(u))*Math.cos(2*Math.PI*v);}
  function mulberry32(seed){return function(){seed|=0;seed=seed+0x6D2B79F5|0;let t=Math.imul(seed^seed>>>15,1|seed);t=t+Math.imul(t^t>>>7,61|t)^t;return((t^t>>>14)>>>0)/4294967296;};}
  function randnSeeded(rng){let u=0,v=0;while(!u)u=rng();while(!v)v=rng();return Math.sqrt(-2*Math.log(u))*Math.cos(2*Math.PI*v);}
  function poissonRng(lam,rng){
    if(lam<=0)return 0;
    if(lam<30){const L=Math.exp(-lam);let k=0,p=1;do{k++;p*=rng();}while(p>L);return k-1;}
    return Math.max(0,Math.round(lam+Math.sqrt(lam)*randnSeeded(rng)));
  }
  function normalPDF(x,mu,s){return Math.exp(-0.5*((x-mu)/s)**2)/(s*Math.sqrt(2*Math.PI));}

  // ═══════════════════ SIR MODEL ═══════════════════
  function solveSIR(R0v,gam,s0,i0,r0){
    const N=s0+i0+r0,beta=R0v*gam/N,dt=1,tMax=364;
    const sol=[{t:0,S:s0,I:i0,R:r0}];
    let S=s0,I=i0,R=r0;
    for(let t=0;t<tMax;t++){
      const fS=(s,i)=>-beta*s*i, fI=(s,i)=>beta*s*i-gam*i;
      const k1s=fS(S,I),k1i=fI(S,I);
      const k2s=fS(S+.5*k1s,I+.5*k1i),k2i=fI(S+.5*k1s,I+.5*k1i);
      const k3s=fS(S+.5*k2s,I+.5*k2i),k3i=fI(S+.5*k2s,I+.5*k2i);
      const k4s=fS(S+k3s,I+k3i),k4i=fI(S+k3s,I+k3i);
      S+=(k1s+2*k2s+2*k3s+k4s)/6; I+=(k1i+2*k2i+2*k3i+k4i)/6;
      S=Math.max(0,Math.min(N,S)); I=Math.max(0,Math.min(N,I)); R=N-S-I;
      sol.push({t:t+1,S,I,R});
    }
    return sol;
  }
  function weeklyInc(sol){
    const w=[];for(let i=0;i<52;i++){let s=0;for(let d=0;d<7;d++){const j=i*7+d;s+=Math.max(0,sol[j].S-sol[j+1].S);}w.push(s);}return w;
  }
  function logLik(obs,pred){
    let ll=0;for(let i=0;i<obs.length;i++){const lam=Math.max(pred[i],1e-10);ll+=obs[i]*Math.log(lam)-lam-lnGamma(obs[i]+1);}return ll;
  }

  // ═══════════════════ DATA GENERATION ═══════════════════
  const TRUE_R0=1.5, TRUE_GAM=0.2, S0=999, I0=1, R0P=0;
  const R0_LO=1,R0_HI=5, GAM_LO=0.05,GAM_HI=0.8;
  const trueSol=solveSIR(TRUE_R0,TRUE_GAM,S0,I0,R0P);
  const trueWeekly=weeklyInc(trueSol);
  const rng=mulberry32(42);
  const obsData=trueWeekly.map(l=>poissonRng(l,rng));

  function targetLL(r0,gam){
    if(r0<R0_LO||r0>R0_HI||gam<GAM_LO||gam>GAM_HI)return-Infinity;
    try{const s=solveSIR(r0,gam,S0,I0,R0P);const p=weeklyInc(s);const ll=logLik(obsData,p);return isFinite(ll)?ll:-Infinity;}
    catch(e){return-Infinity;}
  }

  // ═══════════════════ PALETTE ═══════════════════
  const C={current:"#1565C0",propose:"#E65100",accept:"#2E7D32",reject:"#C62828",
    data:"#FF9800",dataStroke:"#F57C00",model:"#1565C0",modelLight:"#90CAF9",
    traceR0:"#1565C0",traceGam:"#00897B",trueLine:"#E53935",
    hist:"#1565C0",histStroke:"#0D47A1",histGam:"#00897B",histGamStroke:"#00695C",
    scatter:"#42A5F5",scatterAccept:"#66BB6A",
    txt:"#212121",sub:"#9E9E9E",border:"#E0E0E0",bg:"#FAFAFA"};
  const mono="'SF Mono',SFMono-Regular,Menlo,monospace";

  // ═══════════════════ STATE ═══════════════════
  let chain=[{R0:2.5,gamma:0.15}], stepPhase="idle", stepCount=0, nAcc=0;
  let proposedR0=null, proposedGam=null, curLL=null, propLL=null;
  let logAlpha=0, uRand=0, accepted=false;
  let autoRunning=false, autoTimer=null;
  let sigR0=0.15, sigGam=0.03, autoSpeed=200;

  // ═══════════════════ WRAPPER ═══════════════════
  const outer=document.createElement("div");
  outer.style.cssText=`display:flex;flex-direction:column;align-items:center;gap:12px;
    font-family:'Inter',-apple-system,BlinkMacSystemFont,"Segoe UI",Roboto,sans-serif;
    width:100%;max-width:920px;margin:0 auto;`;
  if(!document.querySelector('link[href*="Inter"]')){
    const lk=document.createElement('link');lk.rel='stylesheet';
    lk.href='https://fonts.googleapis.com/css2?family=Inter:wght@400;500;600;700;800&display=swap';
    document.head.appendChild(lk);
  }
  outer.appendChild(injectStyle());

  // ═══════════════════ CONTROLS ═══════════════════
  const ctrlRow=document.createElement("div");
  ctrlRow.style.cssText="display:flex;gap:10px;width:100%;align-items:center;flex-wrap:wrap;";
  function mkBtn(text,color,big){
    const b=document.createElement("button");
    b.style.cssText=`padding:${big?"10px 28px":"7px 16px"};border-radius:${big?12:10}px;border:none;
      background:linear-gradient(135deg,${color},${color}dd);color:#fff;font-size:${big?15:13}px;font-weight:700;cursor:pointer;
      font-family:inherit;transition:all 0.2s ease;letter-spacing:0.3px;
      ${big?"box-shadow:0 4px 14px "+color+"44;":"box-shadow:0 2px 8px rgba(0,0,0,0.12);"}`;
    b.textContent=text;
    b.addEventListener('mouseenter',()=>{b.style.transform='translateY(-1px)';b.style.boxShadow=`0 6px 20px ${color}55`;});
    b.addEventListener('mouseleave',()=>{b.style.transform='translateY(0)';b.style.boxShadow=big?`0 4px 14px ${color}44`:'0 2px 8px rgba(0,0,0,0.12)';});
    return b;
  }
  const btnStep=mkBtn("▶  Bước tiếp theo",C.current,true);
  const btnAuto=mkBtn("⏩ Chạy tự động","#6A1B9A",false);
  const btnReset=mkBtn("↺ Đặt lại","#78909C",false);
  const stepLabel=document.createElement("span");
  stepLabel.style.cssText=`font-size:13px;font-weight:700;color:${C.sub};font-family:${mono};margin-left:auto;`;
  ctrlRow.append(btnStep,btnAuto,btnReset,stepLabel);
  outer.appendChild(ctrlRow);

  // ── SLIDERS ──
  const SL={};
  SL.sigR0=createSlider("Proposal σ(R₀)",0.01,0.50,0.01,0.15,C.propose,"amber");
  SL.sigGam=createSlider("Proposal σ(γ)",0.005,0.10,0.005,0.030,C.propose,"amber");
  SL.speed=createSlider("Tốc độ (ms/bước)",50,500,10,200,"#6A1B9A","purple");
  const sr1=document.createElement("div");sr1.style.cssText="display:flex;gap:14px;width:100%;";
  sr1.append(SL.sigR0.el,SL.sigGam.el,SL.speed.el);
  outer.appendChild(sr1);

  // ═══════════════════ MAIN SVG ═══════════════════
  const W=920,mainH=360,gap=40;
  const pW=Math.floor((W-gap)/2);
  const M={t:30,b:38,l:50,r:15};
  const iW=pW-M.l-M.r, iH=mainH-M.t-M.b;
  // Trace area: left part = trace line, right part = density histogram
  const densW=60; // width reserved for density
  const tGap=42,tH=Math.floor((mainH-tGap)/2);
  const tM={t:24,b:8},tMb={t:8,b:30};
  const tIH=tH-tM.t-tM.b, tIHb=tH-tMb.t-tMb.b;
  const tIW=pW-M.l-M.r-densW-8; // trace inner width (leave room for density)

  const mainSvg=d3.create("svg").attr("viewBox",`0 0 ${W} ${mainH}`)
    .style("width","100%").style("max-width",W+"px");

  // ── Clip path for incidence plot ──
  const clipId="inc-clip-"+Math.random().toString(36).slice(2,8);
  mainSvg.append("defs").append("clipPath").attr("id",clipId)
    .append("rect").attr("x",0).attr("y",0).attr("width",iW).attr("height",iH);

  // ── LEFT: Incidence plot ──
  const gInc=mainSvg.append("g").attr("transform",`translate(${M.l},${M.t})`);
  const xInc=d3.scaleLinear().domain([1,52]).range([0,iW]);
  let yMax=d3.max(obsData)*1.3;
  let yInc=d3.scaleLinear().domain([0,yMax]).range([iH,0]);
  const incAxisX=gInc.append("g").attr("transform",`translate(0,${iH})`);
  const incAxisY=gInc.append("g");
  function drawIncAxes(){
    incAxisX.call(d3.axisBottom(xInc).ticks(8).tickSize(-iH).tickFormat(d3.format("d")))
      .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
    incAxisY.call(d3.axisLeft(yInc).ticks(6).tickSize(-iW))
      .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
  }
  drawIncAxes();
  gInc.append("text").attr("x",iW/2).attr("y",-12).attr("text-anchor","middle")
    .attr("font-size",13).attr("font-weight",700).attr("fill",C.txt).text("Dữ liệu & Dự đoán mô hình");
  gInc.append("text").attr("x",iW/2).attr("y",iH+32).attr("text-anchor","middle")
    .attr("font-size",10).attr("fill",C.sub).text("Tuần");
  gInc.append("text").attr("transform","rotate(-90)").attr("x",-iH/2).attr("y",-38)
    .attr("text-anchor","middle").attr("font-size",10).attr("fill",C.sub).text("Số ca mới / tuần");

  // clipped content group
  const gIncClip=gInc.append("g").attr("clip-path",`url(#${clipId})`);
  gIncClip.selectAll(".obs-dot").data(obsData).join("circle").attr("class","obs-dot")
    .attr("cx",(d,i)=>xInc(i+1)).attr("cy",d=>yInc(d)).attr("r",3.5)
    .attr("fill",C.data).attr("stroke",C.dataStroke).attr("stroke-width",1).attr("opacity",0.9);
  const modelLine=gIncClip.append("path").attr("fill","none").attr("stroke",C.model).attr("stroke-width",2.5).attr("opacity",0.85);
  const propLine=gIncClip.append("path").attr("fill","none").attr("stroke",C.propose).attr("stroke-width",2)
    .attr("stroke-dasharray","6,4").attr("opacity",0);

  // ── RIGHT: Trace plots + density ──
  function makeTracePanel(g,color,histColor,histStroke,label,trueVal,innerH,margins){
    const xTr=d3.scaleLinear().domain([0,10]).range([0,tIW]);
    const yTr=d3.scaleLinear().domain([trueVal-0.5,trueVal+0.5]).range([innerH,0]);
    const axX=g.append("g").attr("transform",`translate(0,${innerH})`);
    const axY=g.append("g");
    g.append("text").attr("x",tIW/2).attr("y",-8).attr("text-anchor","middle")
      .attr("font-size",11).attr("font-weight",700).attr("fill",C.txt).text(label);
    const trueLn=g.append("line").attr("x1",0).attr("x2",tIW)
      .attr("stroke",C.trueLine).attr("stroke-width",1.5).attr("stroke-dasharray","4,3").attr("opacity",0.7);
    const tracePath=g.append("path").attr("fill","none").attr("stroke",color).attr("stroke-width",1.5);
    const curDot=g.append("circle").attr("r",4).attr("fill",color).attr("stroke","#fff").attr("stroke-width",2);
    const propDot=g.append("circle").attr("r",4).attr("fill","none").attr("stroke",C.propose).attr("stroke-width",2).attr("opacity",0);
    // density histogram group (to the right of trace)
    const gDens=g.append("g").attr("transform",`translate(${tIW+8},0)`);
    const densPath=gDens.append("path").attr("fill",histColor).attr("fill-opacity",0.35).attr("stroke",histStroke).attr("stroke-width",1.2);
    const densTrueLn=gDens.append("line").attr("stroke",C.trueLine).attr("stroke-width",1.5).attr("stroke-dasharray","4,3").attr("opacity",0.7);
    return {xTr,yTr,axX,axY,trueLn,tracePath,curDot,propDot,innerH,margins,trueVal,densPath,densTrueLn,gDens};
  }
  const gTrR0=mainSvg.append("g").attr("transform",`translate(${pW+gap+M.l},${tM.t})`);
  const gTrGam=mainSvg.append("g").attr("transform",`translate(${pW+gap+M.l},${tH+tGap+tMb.t})`);
  const trR0=makeTracePanel(gTrR0,C.traceR0,C.hist,C.histStroke,"Chuỗi R₀",TRUE_R0,tIH,tM);
  const trGam=makeTracePanel(gTrGam,C.traceGam,C.histGam,C.histGamStroke,"Chuỗi γ",TRUE_GAM,tIHb,tMb);
  gTrGam.append("text").attr("x",tIW/2).attr("y",tIHb+25).attr("text-anchor","middle")
    .attr("font-size",10).attr("fill",C.sub).text("Bước");

  outer.appendChild(mainSvg.node());

  // ═══════════════════ STATUS BOX ═══════════════════
  const statusBox=document.createElement("div");
  statusBox.style.cssText=`width:100%;padding:14px 18px;border-radius:14px;border:1.5px solid ${C.border};
    background:#fff;font-family:${mono};font-size:13px;color:${C.txt};min-height:40px;
    box-shadow:0 2px 10px rgba(0,0,0,0.04);transition:border-color 0.3s;`;
  outer.appendChild(statusBox);

  // ═══════════════════ SCATTER SVG ═══════════════════
  const scW=420,scH=280,scM={t:20,b:35,l:50,r:15};
  const scIW=scW-scM.l-scM.r,scIH=scH-scM.t-scM.b;
  const scatSvg=d3.create("svg").attr("viewBox",`0 0 ${scW} ${scH}`)
    .style("width","100%").style("max-width",scW+"px");
  const gScat=scatSvg.append("g").attr("transform",`translate(${scM.l},${scM.t})`);
  const xSc=d3.scaleLinear().domain([GAM_LO,GAM_HI]).range([0,scIW]);
  const ySc=d3.scaleLinear().domain([R0_LO,R0_HI]).range([scIH,0]);
  gScat.append("g").attr("transform",`translate(0,${scIH})`).call(d3.axisBottom(xSc).ticks(6).tickSize(-scIH))
    .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
  gScat.append("g").call(d3.axisLeft(ySc).ticks(6).tickSize(-scIW))
    .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
  gScat.append("text").attr("x",scIW/2).attr("y",-6).attr("text-anchor","middle")
    .attr("font-size",12).attr("font-weight",700).attr("fill",C.txt).text("Không gian tham số");
  gScat.append("text").attr("x",scIW/2).attr("y",scIH+30).attr("text-anchor","middle")
    .attr("font-size",10).attr("fill",C.sub).text("γ");
  gScat.append("text").attr("transform","rotate(-90)").attr("x",-scIH/2).attr("y",-38)
    .attr("text-anchor","middle").attr("font-size",10).attr("fill",C.sub).text("R₀");
  gScat.append("line").attr("x1",xSc(TRUE_GAM)-8).attr("x2",xSc(TRUE_GAM)+8)
    .attr("y1",ySc(TRUE_R0)).attr("y2",ySc(TRUE_R0)).attr("stroke",C.trueLine).attr("stroke-width",2.5);
  gScat.append("line").attr("x1",xSc(TRUE_GAM)).attr("x2",xSc(TRUE_GAM))
    .attr("y1",ySc(TRUE_R0)-8).attr("y2",ySc(TRUE_R0)+8).attr("stroke",C.trueLine).attr("stroke-width",2.5);
  const scatDots=gScat.append("g");
  const scatCur=gScat.append("circle").attr("r",5).attr("fill",C.current).attr("stroke","#fff").attr("stroke-width",2);
  const scatProp=gScat.append("circle").attr("r",5).attr("fill","none").attr("stroke",C.propose).attr("stroke-width",2).attr("opacity",0);

  outer.appendChild(scatSvg.node());

  // ═══════════════════ RENDER FUNCTIONS ═══════════════════
  const lineGen=d3.line().x((_,i)=>xInc(i+1)).curve(d3.curveMonotoneX);

  function renderIncidence(){
    const cur=chain[chain.length-1];
    const sol=solveSIR(cur.R0,cur.gamma,S0,I0,R0P);
    const pred=weeklyInc(sol);
    // dynamically adjust y-axis to fit both data and model
    const allVals=[...obsData,...pred];
    let propPred=null;
    if(stepPhase==="proposed"&&proposedR0!=null){
      const pSol=solveSIR(proposedR0,proposedGam,S0,I0,R0P);
      propPred=weeklyInc(pSol);
      allVals.push(...propPred);
    }
    const newYMax=d3.max(allVals)*1.2;
    yInc.domain([0,Math.max(newYMax,10)]);
    drawIncAxes();
    // re-position data dots with new scale
    gIncClip.selectAll(".obs-dot").attr("cy",d=>yInc(d));
    modelLine.attr("d",lineGen.y(d=>yInc(d))(pred));
    if(propPred){
      propLine.attr("d",lineGen.y(d=>yInc(d))(propPred)).attr("opacity",0.8);
    } else propLine.attr("opacity",0);
  }

  function updateTrace(tr,vals,propVal,key){
    const n=vals.length;
    tr.xTr.domain([0,Math.max(10,n)]);
    const ext=d3.extent(vals.map(v=>v[key]));
    const pad=Math.max((ext[1]-ext[0])*0.2,0.05);
    const lo=Math.min(ext[0],tr.trueVal)-pad, hi=Math.max(ext[1],tr.trueVal)+pad;
    tr.yTr.domain([lo,hi]);
    tr.axX.call(d3.axisBottom(tr.xTr).ticks(5).tickFormat(d3.format("d")).tickSize(-tr.innerH))
      .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
    tr.axY.call(d3.axisLeft(tr.yTr).ticks(4).tickSize(-tIW))
      .call(g=>{g.selectAll(".tick line").attr("stroke","#eee");g.select(".domain").attr("stroke",C.border);});
    tr.trueLn.attr("y1",tr.yTr(tr.trueVal)).attr("y2",tr.yTr(tr.trueVal));
    const line=d3.line().x((_,i)=>tr.xTr(i)).y(d=>tr.yTr(d[key]));
    tr.tracePath.attr("d",line(vals));
    const last=vals[vals.length-1];
    tr.curDot.attr("cx",tr.xTr(n-1)).attr("cy",tr.yTr(last[key]));
    if(stepPhase==="proposed"&&propVal!=null){
      tr.propDot.attr("cx",tr.xTr(n)).attr("cy",tr.yTr(propVal)).attr("opacity",1);
    } else tr.propDot.attr("opacity",0);
    // ── Density histogram (rotated 90°, drawn to the right of trace) ──
    if(n>=3){
      const data=vals.map(v=>v[key]);
      const mean=d3.mean(data), sd=d3.deviation(data)||0.01;
      const nBins=30;
      const yLo=tr.yTr.domain()[0], yHi=tr.yTr.domain()[1];
      const step=(yHi-yLo)/nBins;
      const pts=[];
      for(let i=0;i<=nBins;i++){
        const v=yLo+i*step;
        pts.push({v,d:normalPDF(v,mean,sd)});
      }
      const xDens=d3.scaleLinear().domain([0,d3.max(pts,p=>p.d)]).range([0,densW]);
      const area=d3.area().x0(0).x1(p=>xDens(p.d)).y(p=>tr.yTr(p.v)).curve(d3.curveBasis);
      tr.densPath.attr("d",area(pts));
      tr.densTrueLn.attr("x1",0).attr("x2",densW)
        .attr("y1",tr.yTr(tr.trueVal)).attr("y2",tr.yTr(tr.trueVal));
    } else {
      tr.densPath.attr("d","");
    }
  }

  function renderTraces(){
    updateTrace(trR0,chain,proposedR0,"R0");
    updateTrace(trGam,chain,proposedGam,"gamma");
  }

  function renderScatter(){
    scatDots.selectAll("circle").data(chain).join("circle")
      .attr("cx",d=>xSc(d.gamma)).attr("cy",d=>ySc(d.R0))
      .attr("r",2.5).attr("fill",C.scatter).attr("opacity",0.4);
    const last=chain[chain.length-1];
    scatCur.attr("cx",xSc(last.gamma)).attr("cy",ySc(last.R0));
    if(stepPhase==="proposed"&&proposedR0!=null){
      scatProp.attr("cx",xSc(proposedGam)).attr("cy",ySc(proposedR0)).attr("opacity",1);
    } else scatProp.attr("opacity",0);
  }

  function pill(label,color){return `<span style="display:inline-block;padding:3px 12px;border-radius:6px;background:${color};color:#fff;font-weight:700;font-size:12px;letter-spacing:0.5px;">${label}</span>`;}
  function vl(t,c,b){return `<span style="color:${c};font-weight:${b?800:600};">${t}</span>`;}
  function lb(t,c){return `<span style="color:${c};font-size:11px;font-weight:600;opacity:0.7;">${t}</span>`;}

  function updateStatus(){
    const cur=chain[chain.length-1];
    const stats=`<div style="margin-top:6px;font-size:11px;color:${C.sub};">Bước ${stepCount} · ${nAcc} chấp nhận</div>`;
    if(stepPhase==="idle"){
      statusBox.style.borderColor=C.border;
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:12px;">`+pill('SẴN SÀNG','#546E7A')+
        `<span style="font-size:14px;">Vị trí hiện tại: ${vl('R₀ = '+cur.R0.toFixed(3),C.current,true)}, ${vl('γ = '+cur.gamma.toFixed(4),C.traceGam,true)}</span>`+
        `<span style="margin-left:auto;color:${C.sub};font-size:12px;">Nhấn <b style="color:${C.current}">Bước tiếp theo</b> ▶</span></div>`+stats;
    } else if(stepPhase==="proposed"){
      statusBox.style.borderColor=C.propose;
      const clampedAlpha=Math.min(1,Math.exp(logAlpha));
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:10px;margin-bottom:8px;">`+
        pill('ĐỀ XUẤT',C.propose)+`<span style="font-size:13px;">R₀* = ${vl(proposedR0.toFixed(3),C.propose,true)}, γ* = ${vl(proposedGam.toFixed(4),C.propose,true)}</span></div>`+
        `<table style="width:100%;border-collapse:separate;border-spacing:0;border-radius:10px;overflow:hidden;border:1px solid #E0E0E0;font-size:12px;">`+
        `<thead><tr style="background:#F5F5F5;"><th style="padding:6px 12px;text-align:left;border-bottom:1px solid #E0E0E0;"></th>`+
        `<th style="padding:6px 12px;text-align:center;font-weight:700;color:${C.current};border-bottom:1px solid #E0E0E0;">Hiện tại</th>`+
        `<th style="padding:6px 12px;text-align:center;font-weight:700;color:${C.propose};border-bottom:1px solid #E0E0E0;">Đề xuất</th></tr></thead><tbody>`+
        `<tr><td style="padding:5px 12px;font-weight:600;">Log-Likelihood</td>`+
        `<td style="padding:5px 12px;text-align:center;color:${C.current};">${curLL.toFixed(2)}</td>`+
        `<td style="padding:5px 12px;text-align:center;color:${C.propose};">${(propLL===-Infinity?"-∞":propLL.toFixed(2))}</td></tr></tbody></table>`+
        `<div style="display:flex;align-items:center;gap:20px;margin-top:8px;padding:8px 12px;background:#fff;border-radius:8px;border:1px solid #E0E0E0;">`+
        `<span>${lb('Tỷ lệ chấp nhận','#333')} α = ${vl(clampedAlpha.toFixed(4),'#333',true)}</span>`+
        `<span style="color:#BDBDBD;">│</span>`+
        `<span>${lb('Rút ngẫu nhiên','#333')} u = ${vl(uRand.toFixed(4),'#E65100',true)}</span>`+
        `<span style="margin-left:auto;font-size:12px;color:${C.sub};">Nhấn <b>Bước tiếp theo</b> để quyết định</span></div>`;
    } else {
      const dc=accepted?C.accept:C.reject;
      statusBox.style.borderColor=dc;
      const clampedAlpha=Math.min(1,Math.exp(logAlpha));
      statusBox.innerHTML=`<div style="display:flex;align-items:center;gap:12px;">`+
        pill(accepted?'✓ CHẤP NHẬN':'✗ TỪ CHỐI',dc)+
        `<span style="font-size:14px;font-weight:600;">α = ${vl(clampedAlpha.toFixed(4),C.txt)} `+
        `<span style="color:${dc};font-weight:800;font-size:15px;">${accepted?'≥':'<'}</span> `+
        `u = ${vl(uRand.toFixed(4),C.propose)}</span>`+
        `<span style="margin-left:auto;font-size:13px;color:${dc};font-weight:600;">`+
        (accepted?`→ Di chuyển đến R₀*=${proposedR0.toFixed(3)}, γ*=${proposedGam.toFixed(4)}`
                 :`→ Ở lại R₀=${cur.R0.toFixed(3)}, γ=${cur.gamma.toFixed(4)}`)+`</span></div>`+stats;
    }
    stepLabel.textContent=`Bước ${stepCount} · ${nAcc} mẫu`;
  }

  function renderAll(){renderIncidence();renderTraces();renderScatter();updateStatus();}

  // ═══════════════════ MCMC LOGIC ═══════════════════
  function doStep(){
    if(stepPhase==="idle"){
      sigR0=SL.sigR0.val(); sigGam=SL.sigGam.val();
      const cur=chain[chain.length-1];
      proposedR0=cur.R0+sigR0*randn();
      proposedGam=cur.gamma+sigGam*randn();
      curLL=targetLL(cur.R0,cur.gamma);
      propLL=targetLL(proposedR0,proposedGam);
      logAlpha=propLL-curLL;
      uRand=Math.random();
      stepPhase="proposed";
    } else if(stepPhase==="proposed"){
      accepted=Math.log(uRand)<logAlpha;
      if(accepted) chain.push({R0:proposedR0,gamma:proposedGam});
      else chain.push({R0:chain[chain.length-1].R0,gamma:chain[chain.length-1].gamma});
      stepCount++; if(accepted)nAcc++;
      stepPhase="decided";
    } else {
      stepPhase="idle";
    }
    renderAll();
  }

  function resetAll(){
    stopAuto();
    chain=[{R0:2.5,gamma:0.15}]; stepPhase="idle"; stepCount=0; nAcc=0;
    proposedR0=null; proposedGam=null; curLL=null; propLL=null;
    renderAll();
  }

  function stopAuto(){autoRunning=false;clearInterval(autoTimer);btnAuto.textContent="⏩ Chạy tự động";btnAuto.style.background="linear-gradient(135deg,#6A1B9A,#6A1B9Add)";}

  // ═══════════════════ EVENT HANDLERS ═══════════════════
  btnStep.addEventListener("click",()=>{if(autoRunning)stopAuto();doStep();});
  btnReset.addEventListener("click",resetAll);
  btnAuto.addEventListener("click",()=>{
    if(autoRunning){stopAuto();return;}
    autoRunning=true;btnAuto.textContent="⏸ Tạm dừng";btnAuto.style.background=`linear-gradient(135deg,${C.reject},${C.reject}dd)`;
    autoSpeed=SL.speed.val();
    autoTimer=setInterval(doStep,autoSpeed);
  });
  // When speed slider changes during auto-run, restart timer with new speed
  SL.speed.input.addEventListener("input",()=>{
    SL.speed.sync();
    if(autoRunning){clearInterval(autoTimer);autoSpeed=SL.speed.val();autoTimer=setInterval(doStep,autoSpeed);}
  });
  for(const k of Object.keys(SL)){if(k!=="speed")SL[k].input.addEventListener("input",()=>SL[k].sync());}

  // ═══════════════════ INITIAL RENDER ═══════════════════
  renderAll();
  return outer;
}