30:T37d3,class stats { /** Basic Statistical Operations * sum * mean * median * mad * sumSquares * variance * stDeviation * prod */ // -----------------------------------------------------------------------------------------------||| static sum (arr) { return arr.reduce( (acc,val) => acc + val, 0); } // -----------------------------------------------------------------------------------------------||| static mean (arr) { let len = arr.length; return stats.sum(arr) / len; } // -----------------------------------------------------------------------------------------------||| static median(arr) { const n = arr.length; const vals = [...arr]; vals.sort((a,b) => a-b); const middleIndex = parseInt(n/2); const isEven = !(n % 2); return isEven ? stats.mean(vals.splice(middleIndex-1,2)) : vals[middleIndex]; } // -----------------------------------------------------------------------------------------------||| static mad(arr,mode="mean",deviation="median") { const outerFx = mode === "mean"? stats.mean : stats.median; const innerFx = deviation === "median"? stats.median : stats.mean; return outerFx(arr.map( v => {return Math.abs(v-innerFx(arr))} )); } // -----------------------------------------------------------------------------------------------||| static sumSquares(arr,internal_stat = stats.mean) { if (typeof internal_stat !== "function") { throw new Error("Internal Statistic must be a function which accepts a single array arg."); } const sub = internal_stat(arr); return arr.map( v => (v-sub)**2 ).reduce( (sum,dif) => sum + dif, 0); } // -----------------------------------------------------------------------------------------------||| static variance(arr,bias="unbiased") { // bias can be one of: // unbiased, biased, uncorrected // which correspond to the scaling by: // N-1, N+1, N const ssq = stats.sumSquares(arr,stats.mean); let div = arr.length - 1; switch (bias) { case "biased": div = arr.length + 1; break; case "uncorrected": div = arr.length; break; } return !div? 0 : ssq/div; } // -----------------------------------------------------------------------------------------------||| static stDeviation(arr,...varArgs) { // standard deviation return Math.sqrt(stats.variance(arr,...varArgs)); } // -----------------------------------------------------------------------------------------------||| static prod(arr) { let newArray = arr.filter( v => v !== 0 ); return newArray.reduce( (a,v) => a * v, 1 ); } // -----------------------------------------------------------------------------------------------||| /** distribution based operations and sampling * randInt * randNorm (TODO) * cdfNormal * ppfNormal * ecdf * quantile */ //Pseudo-random.. use Math.random() though I'd rather use a pcg or xoroshiro... static randInt (size=1,{min=0,max=1,inclusive=true}={}) { // validate range min = Math.floor(min); max = Math.ceil(max); // setup array const randNums = []; for (let i=0; i= 0 ? r : 2 - r; }; v = typeof v === "number"? [v] : v; return v.map(x => 0.5 * erfc(-(x - mu) / (Math.sqrt(nu) * Math.sqrt(2)))); } // -----------------------------------------------------------------------------------------------||| static ppfNormal(q,mu,nu){ // Complementary error function // From Numerical Recipes in C 2e p221 function erfc(x) { var z = Math.abs(x); var t = 1 / (1 + z / 2); var r = t * Math.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277))))))))) return x >= 0 ? r : 2 - r; }; // Inverse complementary error function // From Numerical Recipes 3e p265 function ierfc(x) { if (x >= 2) { return -100; } if (x <= 0) { return 100; } var xx = (x < 1) ? x : 2 - x; var t = Math.sqrt(-2 * Math.log(xx / 2)); var r = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t); for (var j = 0; j < 2; j++) { var err = erfc(r) - xx; r += err / (1.12837916709551257 * Math.exp(-(r * r)) - r * err); } return (x < 1) ? r : -r; }; q = typeof q === "number"? [q] : q; return q.map(x => mu - Math.sqrt(nu) * Math.sqrt(2) * ierfc(2 * x) ); } // -----------------------------------------------------------------------------------------------||| static sample(arr,size,replace=false,prob=null){ //validate if (!replace && size > arr.length) { throw new Error("Cannot permute more than available when 'replace = false'."); } const n = arr.length; if (prob) { // we need to weight out values prob = typeof prob === "number"? Array.from({length:n}, d => prob ): prob; if (prob.length != n) { throw new Error("Probability weights must be scalar or same length as arr."); } // convert to fractions let pSum = prob.reduce( (acc,val) => acc+val, 0); prob = prob.map( v => v/pSum ); } else { prob = Array.from({length:n}, d => 1/n); } // map probs to counts for each element, fractions get rounded up // which may alter probability. let probNums = prob.map( v => v*n ); // make sure we get some accuracy for small n and wild probs let minProbs = Math.min(...probNums); let fact = 1; if (minProbs < 1) { while(minProbs < 5){ fact *= 10**0.25; minProbs *= 10**0.25; } } else { fact = 1; } probNums = probNums.map( v => Math.ceil(v*fact)); // create expanded values array const values = arr.map( (v,i) => { const newValues = []; for (let j=1; j<=probNums[i]; j++){ newValues.push(v); } return newValues; }).flat(2); let n_possible = values.length; const allInds = stats.seq(n_possible,0,{by:1}); const indexes = []; while (indexes.length < size) { let testIndex = stats.randInt(1,{min:0,max:n_possible,inclusive:false}); const val = allInds[testIndex]; if (!replace) { allInds.splice(testIndex,1); n_possible -= 1; } indexes.push(val); } const output = []; for (let i = 0; i < size; i++) { output.push(values[indexes[i]]); } return output; } // -----------------------------------------------------------------------------------------------||| static ecdf(arr) { // returns a function which accepts an array of data and returns the empirical percentile const x = [...arr].filter( v => !isNaN(v)).sort((a,b) => a-b); const vals = stats.unique(x,false); const n = x.length; // get percentiles from data // here, duplicates are accounted for by the match (reduce) const frequency = vals.map( v => x.reduce((count,u) => count+(u===v),0)); const cumFreq = stats.cumsum(frequency).map(v=>v/n); cumFreq.unshift(0.0); vals.unshift(-Infinity); // define a function that uses "constant" interpolation const fx = (x,y) => (v) => { if (isNaN(v)) return NaN; const n = x.length; var i = 0; var j = n - 1 ; if (v < x[i]) return 0; if (v > x[j]) return 1; // search the array for the two elements closest to the value while (i < j-1) { let ij = Math.floor((i + j)/2); if(v < x[ij]) { j = ij; } else { i = ij; } } if(v === x[j]) return y[j]; return y[i]; }; // return the function return fx(vals,cumFreq); } // -----------------------------------------------------------------------------------------------||| static quantile(arr,probs,type=8) { // approxmate probs using linear interpolation between points const vals = [...arr].sort((a,b)=>a-b); const n = vals.length; let m,gCal; // define convenience functions. const k = (p,m) => n*p + m(p); let g = (p,m) => k(p,m)-Math.floor(k(p,m)); // Parameters based on type switch (type) { case 1: m = (p) => 0; gCal = (p,m) => (g(p,m) > 0)? 1 : 0; break; case 2: m = (p) => 0; gCal = (p,m) => (g(p,m) > 0)? 1 : 1/2; break; case 3: m = (p) => -1/2; gCal = (p,m) => ((g(p,m) == 0) && !(Math.floor(k(p,m))%2))? 0 : 1; break; case 4: m = (p) => 0; break; case 5: m = (p) => 1/2; break; case 6: m = (p) => p; break; case 7: m = (p) => (1-p); break; case 9: m = (p) => (p+3/8)/4; break; default: // 8 m = (p) => (p+1)/3; break; } return [...probs].flat().map( p => { let j = Math.floor(k(p,m)); let gamma = type >= 4? g(p,m): gCal(p,m); return (1-gamma)*vals[j-1] + gamma*vals[j]; }); } // -----------------------------------------------------------------------------------------------||| /** Methods that return arrays (copies from inputs) after modification * unique * cumsum * diff * movingAverage * trapezoid (:= discrete integral) * seq */ // array operations static unique(arr,sort=false) { let output = {}; let vals = arr.filter(v => !(output[v]=v in output) ); vals = sort? vals.sort((a,b) => a-b) : vals; return vals; } // -----------------------------------------------------------------------------------------------||| static cumsum(arr,sort=false) { let newArray = [...arr];//copy if (sort) { newArray.sort(); } let output = []; let summed = newArray.reduce( (acc,val) => { output.push( acc ); return acc + val; }, 0); return output.slice(1).concat(summed); } // -----------------------------------------------------------------------------------------------||| static diff(arr) { let i = 0; let result = []; while (i < arr.length-1) { result.push(arr[++i] - arr[i-1]);//arr[1]-arr[1-1] for i=0 } return result; } // -----------------------------------------------------------------------------------------------||| static movingAverage(arr,n=2) { // very simple moving average calculation // no interpolation is applied at the array extremes let result = []; let i = 0; while (i < arr.length) { let end = i+n; if (end > arr.length) { break; } let currentArray = arr.slice(i++,end); result.push(stats.mean(currentArray)); } return result; } // -----------------------------------------------------------------------------------------------||| static trapezoid(X,Y) { if (arguments.length < 2) { Y = X; X = stats.seq(Y.length);// use interval == 1 } else { // trapezoid(X,Y) if (typeof X === "number") { //X is scalar X = stats.seq(Y.length*X,X,{length:Y.length,by:X});//interval == X } if (X.length !== Y.length) { throw new Error("X must be scalar or same length as Y"); } } // Cacl sigfig let dec = Math.max(...Y.map(v=> { return v.toString().split(".").filter( (v,i) => !!i ) }).map( v => v[0].length )); // Calculate Trapz let dx = stats.diff(X); let sy = stats.movingAverage(Y,2); // sum of change in x * avg change in y return parseFloat(sy.reduce( (acc,v,i) => acc + v * dx[i],0).toFixed(dec+1)); } // -----------------------------------------------------------------------------------------------||| static seq(to,from=0,{length=null,by=null,inclusive=false}={}) { // linear spaced vector from to let range = to-from; let direction = Math.sign(range); range = Math.abs(range);//get rid of negatives if (length === null & by === null) { by = parseFloat(range / Math.floor(range)); by = by === Infinity ? range : by; } if (length === null) { length = Math.ceil(range/by)+inclusive; length = length < 2? 2 : length; // can't be zero } if (by === null) { by = range / (inclusive? length-1:length); } let vec = []; for (let i = 0; i < length; i++) { vec.push(from + direction * i * by); } return vec } // Distance metrics static levenshtein(a,b) { if (a.length === 0) return b.length if (b.length === 0) return a.length let tmp, i, j, prev, val, row // swap to save some memory O(min(a,b)) instead of O(a) if (a.length > b.length) { tmp = a a = b b = tmp } row = Array(a.length + 1) // init the row for (i = 0; i <= a.length; i++) { row[i] = i } // fill in the rest for (i = 1; i <= b.length; i++) { prev = i for (j = 1; j <= a.length; j++) { if (b[i-1] === a[j-1]) { val = row[j-1] // match } else { val = Math.min(row[j-1] + 1, // substitution Math.min(prev + 1, // insertion row[j] + 1)) // deletion } row[j - 1] = prev prev = val } row[a.length] = prev } return row[a.length] } }31:T8d2,ecdfDescription = { // ECDF example for thumbnail const x = [0.79, 1.15, 1.61, -0.73, -0.44, -0.55, 1.98, 1.3, 0.03, -1.45, 0.2]; x.sort((a,b)=> a-b ); const fx = stats.ecdf(x); const xECDF = x.map( v => fx(v) ); const x_bar = stats.mean(x); const var_x = stats.variance(x); const pNorm = stats.ppfNormal(xECDF,0,1); const tIQR = stats.quantile(pNorm,[0.25,0.75],7); const rIQR = stats.quantile(x,[0.25,0.75],7); const dx = stats.diff(tIQR); const dy = stats.diff(rIQR); const qqSlope = dy/dx; //y = mx+b ... b = y-mx = int const itcp = ((rIQR[0] - qqSlope*tIQR[0])+(rIQR[1] - qqSlope*tIQR[1]))/2; const qqline = utils.BivariateDatum( d3.extent(pNorm.slice(0,-1)), d3.extent(pNorm.slice(0,-1)).map( v=> v*qqSlope + itcp), "qqline", "lines", {width:2,style: "dotted",color:"rgb(166,166,166)"} ); // plot ecdf const modelX = stats.seq(2.1,-1.6,{length:1000,inclusive:true}); const modelY = modelX.map( v => fx(v)); const ax = linePlot( modelX, modelY, {ylab:"Empircal Distribution",xlab:"Values",dims:[width/2-20,400],mode:"lines"}, {width:2} ); ax.XZero = false; ax.YZero = false; ax.XGrid = true; ax.YGrid = true; // append markers const markers = utils.BivariateDatum(x,x.map(v=>fx(v)),"mk","markers",{},{size:8,color:"rgb(211, 84, 0)"}); ax.appendPlot(markers); // qqplot const qqax = linePlot( pNorm.slice(0,-1), x.slice(0,-1), {ylab:"Sample Quantiles",xlab:"Theoretical Quantiles",dims:[width/2-20,400],mode:"markers"}, {}, {size:8,color:"rgb(10,10,10)"} ); qqax.appendPlot(qqline); qqax.XZero = false; qqax.YZero = false; qqax.XGrid = true; qqax.YGrid = true; // place in box const par1 = ax.here(); par1.style.setProperty("margin", "0 10 0 10"); par1.style.setProperty("float", "left"); const par2 = qqax.here(); par2.style.setProperty("margin", "0 10 0 10"); par2.style.setProperty("float", "left"); const box = DOM.element("div"); box.append(par1); box.append(par2); d3.select(box) .attr("width", width) .attr("height",400) .style("margin", 0) .style("width", width) .style("height",400) .style("overflow", "hidden"); return box; }32:T549,tab ={ function conditionEntry (ret,condition) { const output = [`\\left\\{ \\begin{array}{ll} `]; for (let i = 0; i < ret.length; i++) { output.push(`${ret[i]} & ${condition[i]} \\\\`); } output.push(`\\end{array}\\right. `); return output.join('\n'); } const types = [ { Type:1, m:tex`m=0`.outerHTML, gamma:tex`\gamma = ${conditionEntry([0,1],["g = 0","g \\neq 0"])}`.outerHTML }, { Type:2, m:tex`m=0`.outerHTML, gamma:tex`\gamma = ${conditionEntry(["1/2",1],["g = 0","g \\neq 0"])}`.outerHTML }, { Type:3, m:tex`m=-\frac{1}{2}`.outerHTML, gamma:tex`\gamma = ${conditionEntry([0,1],["g=0\\land j_{even}","g\\neq 0\\lor j_{odd}"])}`.outerHTML }, {Type:4,m:tex`m=0`.outerHTML,gamma:tex`\gamma = g`.outerHTML}, {Type:5,m:tex`m=\frac{1}{2}`.outerHTML,gamma:tex`\gamma = g`.outerHTML}, {Type:6,m:tex`m=p`.outerHTML,gamma:tex`\gamma = g`.outerHTML}, {Type:7,m:tex`m=1-p`.outerHTML,gamma:tex`\gamma = g`.outerHTML}, {Type:8,m:tex`m=\frac{p+1}{3}`.outerHTML,gamma:tex`\gamma = g`.outerHTML}, {Type:9,m:tex`m=\frac{p}{4}+\frac{3}{8}`.outerHTML,gamma:tex`\gamma = g`.outerHTML} ]; const tab = table( types, { enableFilter: false, enableCSVDownload: false, maxHeight: 1000 }); return tab }33:T62d,{ tab.style.setProperty("height","360px","important"); // center rows const td = tab.getElementsByTagName("td"); const n = td.length; for (let i=0; i < n; i++) { td[i].style.setProperty('vertical-align', 'middle', 'important'); } // add title and caption const iDiv = tab.querySelector("div:nth-child(1)") iDiv.style.setProperty("height","330px","important"); let titleDiv = tab.getElementsByClassName("legend"); if(titleDiv.length === 0){ titleDiv = document.createElement("div"); titleDiv.className = "legend"; iDiv.after(titleDiv) } titleDiv.innerHTML = ` Table 1: Definitions for ${tex`m`.outerHTML} and ${tex`\gamma`.outerHTML} from Hyndman and Fan. `; return md` The quantile distribution estimators are a bit of mess. In R, the default implementation is type 7, the case where ${tex`p_k= \mathrm{mode}\left[ F(x_k)\right] `}, however Hyndman and Fan recommend the median-unbiased estimator such that ${tex`p_k\approx \mathrm{median}\left[ F(x_k)\right] .`} The benefit of type 8, they argue, is that estimates are distribution free and robust. A comparison of quantile types for the set ${tex.block`\left\{ ${qData.join(", ")}\right\} `} is shown below. D3.js also includes a \`\`quantile()\`\` method which, without having looked at any of the code, appears to match the output of type 7, similar to the default case in R. One caveat to the D3.js approach is that the input data must be sorted (ascending: \`\`x.sort( (a,b) => a - b );\`\`) beforehand otherwise out-of-order values are returned. ` }34:Ta8b,quantileDescription = { function arrow(color) { const x = d3.scaleLinear().domain([0,1]).range([0.5,14-0.5]); const y = d3.scaleLinear().domain([0,1]).range([17-0.5,0.5]); const ctx = DOM.context2d(14,17); const line = d3.line().x((d,i) => x(i)).y(y).context(ctx); ctx.beginPath(); ctx.lineWidth = 4; line([0.5,0.5]); ctx.strokeStyle = color; ctx.shadowOffsetX = 3/4; ctx.shadowOffsetY = 3/2; ctx.shadowBlur = 3; ctx.shadowColor = "rgba(104, 104, 104, 0.12)"; ctx.stroke(); return ctx.canvas; } function hexToRgb(hex) { let result = /^#?([a-f\d]{2})([a-f\d]{2})([a-f\d]{2})$/i.exec(hex); return result ? { r: parseInt(result[1], 16), g: parseInt(result[2], 16), b: parseInt(result[3], 16) } : null; } // const x = [1.36, 0.6, 0.72, -1.13, -0.25]; const x = [...qData]; // copy q x.sort((a,b) => a-b); const probs = stats.seq(1,0,{length:50,inclusive:true}); const types = stats.seq(10,1); const colors = [ "#e6194B","#f58231","#ffe119","#bfef45","#3cb44b", "#42d4f4","#4363d8","#911eb4","#f032e6" ]; const data = types.map( (t,i) => { let estimates = stats.quantile(x,probs,t); let col = hexToRgb(colors[i]); return utils.BivariateDatum( estimates, probs, "type"+t, "lines", { //lineOpts width:2, color:`rgb(${col.r},${col.g},${col.b})`, opacity:0.8, style: t < 4? "dashed" : "solid" } ); }); data.push( //d3.quantile variant utils.BivariateDatum( probs.map( q => d3.quantile(x,q)), probs, "d3Quantile", "lines", {width:2,color:`rgb(20,20,20)`,opacity:0.8,style:"dotted"} ) ); // ecdf for original quantiles const cdf = stats.ecdf(x); const y = x.map( v => 0); const ax = linePlot( x, y, {ylab:"Quantiles",xlab:"Data (Estimates)",dims:[width-20,600],mode:"markers"}, {}, {size:12,color:"rgb(0, 0, 0)"} ); ax.XZero = false; ax.YZero = true; ax.XGrid = true; ax.YGrid = true; ax.appendPlot(data); const legend = html`

${data.map( o => html`${o.name}: ${arrow(o.line.color)} `)}

`; legend.style.setProperty("background", "#FFFFFF"); legend.style.setProperty("display","block","important"); legend.style.setProperty("position","absolute","important"); legend.style.setProperty("margin","-12px 0px 5px 65px","important"); legend.style.setProperty("width", width-20); legend.style.setProperty("height", 45); const parent = ax.here(); parent.append(legend); parent.style.setProperty("margin","auto"); return parent; }35:Tbb1,md` ## Quantile ~ Quantile Plot Now, imagine you wanted to to visualize how your data compare to the expected quantiles of another distribution, how might you go about this? One handy method is to compare your data against the corresponding quantiles of another distribution. Often times, we want to know how our data compare to the standard normal distribution, i.e., when we want to ask if our data are normally distributed, or some other distribution. We can turn to the Quantile\~Quantile plot, or QQ plot for short. The QQ plot is a simple mechanism by which we can correlated the quantiles of one dataset with another. In the case of sample $A$ versus distribution $B$, we can simply sort the data, and assign each point a weight of ${tex`\frac{1}{N}`}, assuming the data have equal length. If the lengths are not the same, we then make some linear interpolation of the datasets to make them the same length. This effectively gives us data points that correspond to ${tex`\frac{0:(N-1)}{N-1}`} quantiles. We then plot these these data, $A$ against $B$, in a scatter plot. To assess the relationship between the sample and the distribution, we plot a diagonal line, i.e., the identity line, and inspect where the points are in relation to the line. _How do we make a QQ plot?_ This is a question for which the answer is not all that straight forward. Generally speaking, a QQ plot may be constructed from a sample and the corresponding quantiles of the distribution of interest. In the common case of testing whether the sample may be treated as being sampled from a normal distribution, the distribution is the standard normal distribution. In this case of testing data against the standard normal distribution, we need to consider that the standard normal distribution has tails that stretch into ${tex`(-\infty,\, \infty)`}. Thus, if our data are not at either of those extremes, we must modify their quantiles such that we achieve a closeness to the standard normal distribution. For small sample sizes (${tex`N\lt 10`}), it is common to adjust the percentiles by a factor of 0.375, or ${tex`\frac{3}{8}`}, such that the data quantiles are evenly spaced between 0 and 1, e.g., for ${tex`N=7`} we would have ${tex`\frac{5}{2(4N+1)}=0.0862`} for the lowest element, and ${tex`\frac{2n-1}{2n}=0.9138`} for the highest element. In other word, the quantile fraction is adjust by an offset of ${tex`\alpha`} such that probability points along the empirical cumulative distribution function (ECDF) would be set at ${tex`\frac{\{1,\, \ldots,\, N\} - \alpha}{N+1-2\alpha}`}. So let's build a QNorm or (QQ plot of our quantiles against a normal distribution). A common approch is to use the standard normal distribution where ${tex`\mathcal{N}(\mu=0,\sigma^2=1)`}. I have read arguments for and against this approach, suffice to say, it doesn't seem to make much difference whether you use the standard normal (below, right), or simulated normal from your own data (below, left). `36:Taa2,qqplot = { const x = [...qData].sort((a,b) => a-b); const N = x.length; const fx = stats.ecdf(x); const xECDF = x.map( v => fx(v) ); const x_bar = stats.mean(x); const var_x = stats.variance(x); // FOR QQ STANDARD NORM const alpha = 3/8; const ppoints = stats.seq(N,1,{length:N,inclusive:true}).map(v => (v-alpha)/(N+1-2*alpha)); const pNorm = stats.ppfNormal(ppoints,0,1); const tIQR = stats.quantile(pNorm,[0.25,0.75],9); const rIQR = stats.quantile(x,[0.25,0.75],9); const dx = stats.diff(tIQR); const dy = stats.diff(rIQR); const qqSlope = dy/dx; //y = mx+b ... b = y-mx = int const itcp = ((rIQR[0] - qqSlope*tIQR[0])+(rIQR[1] - qqSlope*tIQR[1]))/2; const qqline = utils.BivariateDatum( d3.extent(pNorm.slice(0,-1)), d3.extent(pNorm.slice(0,-1)).map( v=> v*qqSlope + itcp), "qqline", "lines", {width:2,style: "dotted",color:"rgb(166,166,166)"} ); const qqax = linePlot( pNorm.slice(0,-1), x.slice(0,-1), {ylab:"Sample Quantiles",xlab:"Standard Quantiles",dims:[width/2-20,400],mode:"markers"}, {}, {size:8,color:"rgb(10,10,10)"} ); qqax.appendPlot(qqline); qqax.XZero = false; qqax.YZero = false; qqax.XGrid = true; qqax.YGrid = true; // FOR QQPlot const qrNorm = stats.ppfNormal(xECDF,x_bar,var_x); const trIQR = stats.quantile(qrNorm,[0.25,0.75],9); const rrIQR = stats.quantile(x,[0.25,0.75],9); const drx = stats.diff(trIQR); const dry = stats.diff(rrIQR); const qrSlope = dry/drx; //y = mx+b ... b = y-mx = int const qritcp = ((rrIQR[0] - qrSlope*trIQR[0])+(rrIQR[1] - qrSlope*trIQR[1]))/2; const qrline = utils.BivariateDatum( d3.extent(qrNorm.slice(0,-1)), d3.extent(qrNorm.slice(0,-1)).map( v=> v*qrSlope + qritcp), "qqline", "lines", {width:2,style: "dotted",color:"rgb(166,166,166)"} ); const qrax = linePlot( qrNorm.slice(0,-1), x.slice(0,-1), {ylab:"Sample Quantiles",xlab:"Theoretical Quantiles",dims:[width/2-20,400],mode:"markers"}, {}, {size:8,color:"rgb(10,10,10)"} ); qrax.appendPlot(qrline); qrax.XZero = false; qrax.YZero = false; qrax.XGrid = true; qrax.YGrid = true; // place in box const par1 = qrax.here(); par1.style.setProperty("margin", "0 10 0 10"); par1.style.setProperty("float", "left"); const par2 = qqax.here(); par2.style.setProperty("margin", "0 10 0 10"); par2.style.setProperty("float", "left"); const box = DOM.element("div"); box.append(par1); box.append(par2); d3.select(box) .attr("width", width) .attr("height",400) .style("margin", 0) .style("width", width) .style("height",400) .style("overflow", "hidden"); return box; }37:T53f,{ const ecdf = (arr) => { // returns a function which accepts an array of data and returns the empirical percentile const x = [...arr].filter( v => !isNaN(v)).sort((a,b) => a-b); const vals = stats.unique(x,false); const n = x.length; // get percentiles from data // here, duplicates are accounted for by the match (reduce) const frequency = vals.map( v => x.reduce((count,u) => count+(u===v),0)); const cumFreq = stats.cumsum(frequency).map(v=>v/n); cumFreq.unshift(0.0); vals.unshift(-Infinity); // define a function that uses "constant" interpolation const fx = (x,y) => (v) => { if (isNaN(v)) return NaN; const n = x.length; var i = 0; var j = n - 1 ; if (v < x[i]) return 0; if (v > x[j]) return 1; // search the array for the two elements closest to the value while (i < j-1) { let ij = Math.floor((i + j)/2); if(v < x[ij]) { j = ij; } else { i = ij; } } if(v === x[j]) return y[j]; return y[i]; }; // return the function return fx(vals,cumFreq); } // const arr = [0.79, 1.15, 1.61, -0.73, -0.44, -0.55, 1.98, 1.3, 0.03, -1.45, -0.44, 0.2]; const efx = ecdf(arr); arr.sort((a,b) => a-b); return {"arr":arr,"ecdf":arr.map( v => efx(v) )} }38:[] 2e:{"title":"Statistics","cells":[{"id":0,"value":"md`# Statistics\n\nThe following stats class is a collection of (mostly) vanilla JS stats methods for use in other notebooks. Sorry the code is a bit of a mess... someday I'll clean it up. Feel free to play around with the functions. Code sources are noted otherwise you're looking at my feeble attempts to wrangle JavaScript according to my will. Scroll to the bottom for usage notes and other information.\n\nThis library can be called by:\n~~~JS\nimport {stats} from \"${document.referrer.replace('https://observablehq.com/','')}\"\n~~~\n\nAs with most of my library-type notebooks... This notebook is under development *in perpetuum*. Enjoy!\n

-Khris

\n\n${(navigator.userAgent.match('HeadlessChrome')) ? html`

` : ''}\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1383,"value":"md`\nThe \\`\\`stats\\`\\` class contains the following static methods:\n${Object.getOwnPropertyNames(stats).filter( p => typeof stats[p] === \"function\" ).map( p => \"- \\`\\`stats.\"+p+\"()\\`\\`\").join(\"\\n\")}\n\n---\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":143,"value":"md`\n# Notes\n\nNotes, ideas or expanded comments on some of the above functions.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":147,"value":"md`\n## \\`\\`trapezoid\\`\\`\n\\`\\`trapezoid(X,Y)\\`\\` or \\`\\`trapezoid(Y)\\`\\`\nIntegrate Y with respect to X. If only a vector is supplied, \\`\\`trapezoid(Y)\\`\\`, then X is considered 1. If X is scalar then \\`\\`Ig = trapezoid(X,Y) = X*trapezoid(Y)\\`\\`. If X is non-scalar, then A) it must be the same length as Y and B) each [X,Y] pair is used to compute the \"integration\". This implementation is more like the \"chained-trapezoid\" rule, specific information (as far as one can trust good 'ol wikipedia) can be found on the wiki: [Trapezoidal Rule](https://en.wikipedia.org/wiki/Trapezoidal_rule).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1106,"value":"md`\n## \\`\\`ecdf\\`\\`\n\nThe [Empirical Cumulative Density Function](https://en.wikipedia.org/wiki/Empirical_distribution_function) (ECDF) is an estimate of the underlying cumulative distribution function. Here we use the ${tex`1/n`} approach for the indicator function: ${tex.block`\\hat{\\mathrm{F}}_n\\left( t\\right) = \\frac{1}{n}\\sum_{i=1}^{n}1_{X_i\\leq t}`}\n\nThe ECDF has many uses in statistics, from data exploration to distribution fitting. There are also a number of use cases for null hypothesis testing... (⚠️ opinion coming in 3,2,1...) though if you can't concisely define what a *p-value* is, you shouldn't be performing any statistical hypothesis testing.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":738,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":816,"value":"md`\n## \\`\\`quantile\\`\\`\n\nSee Hyndman and Fan ([1996](https://www.jstor.org/stable/2684934)) for more info. Here I implement all 9 types using 8 by default as recommended in the text. This implementation is similar to the [R function](https://www.rdocumentation.org/packages/stats/versions/3.6.0/topics/quantile) of the same name.\n\nIn general, the quantile function is defined as: \n${tex.block`Q_{i}(p)=(1- \\gamma)x_{j}+\\gamma x_{j+1}`}\nwhere ${tex`j=\\lfloor k\\rfloor`} for ${tex`k=np+m`} and ${tex`g=k-j`} and ${tex`i`} indicates the type. Parameters ${tex`m`} and ${tex`\\gamma`} are a function of quantile type. By defualt, I implement type 8 such that,\n\n${tex.block`Q_{8}(p)=(1-(k-\\left\\lfloor k \\right\\rfloor ))x_{\\left( \\left\\lfloor k\\right\\rfloor \\right) } + (k-\\left\\lfloor k \\right\\rfloor )x_{\\left( \\left\\lfloor k\\right\\rfloor +1\\right) }`}\n\n${tex.block`\\mathrm{where}\\, k=\\frac{p(3n+1)+1}{3}`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":958,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":969,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1589,"value":"qData = [1.36, 0.6, 0.72, -1.13, -0.25, -0.45, 1.62, 0.01,-0.8,1.01];","pinned":false,"hidden":false,"mode":"ojs"},{"id":1147,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1606,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":1631,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":731,"value":"md`\n---\n\n# Dependencies and Other Libraries\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1160,"value":"d3=require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":735,"value":"import { linePlot, utils } from \"@khlick/reusable-graphics-and-grobs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":917,"value":"import { table } from \"@gampleman/table\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1492,"value":"$37","pinned":true,"hidden":false,"mode":"ojs"},{"id":1505,"value":"{\n const v = [NaN,0.1,0.5,10,NaN,-20];\n return v.filter( val => !isNaN(val));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1546,"value":"---\n\n## Test Bed","pinned":false,"hidden":false,"mode":"md"},{"id":1549,"value":"expandGrid([[\"a\"],[\"b\"],[\"c\"]], [[\"a\"],[\"b\"],[\"c\"],[\"d\"],[\"e\"],[\"f\"]], [[\"a\"],[\"b\"],[\"c\"],[\"d\"]], [[\"enter\"]])","pinned":true,"hidden":false,"mode":"ojs"},{"id":1553,"value":"function expandGrid(...arrays){\n let N = arrays.length;\n let L = arrays.map(a => a.length);\n let T = L.reduce((a,b) => a*b);\n let m = arrays.map((a,i) => {\n // let nCopies = T/L[i]; // each item in array appears nCopies Times\n let nRepeats = i === 0 ? 1 : L.slice(0,i).reduce((a,b)=>a*b);\n // return a.map(v => Array(nRepeats).fill(Array.isArray(v) ? v[0] : v)).flat();\n let base = a.map(v => Array(nRepeats).fill(Array.isArray(v) ? v[0] : v)).flat();\n let nCopies = T/base.length\n return Array(nCopies).fill(base).flat()\n });\n // const zip = (...arr) => Array(Math.max(...arr.map(a => a.length))).fill().map((_,i) => arr.map(a => a[i]));\n return Array(T).fill().map((_,i) => m.map(v => v[i]))\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W38","stdlib":"1","titling":"auto"}