/* StatsBox: self-defined statstical functions for Ch. 1-3 of * book https://lindnerdrwg.github.io/Statistics-Maxima-Companion.pdf * by Dr. W.G. Lindner, Leichlingen DE, 2026 */ fpprintprec : 5$ load(distrib)$ dim(X) := length(X)$ Sum(X) := apply("+", X)$ Mean(X) := Sum(X)/dim(X) $ sd(X) := float( sqrt(Sum((X-Mean(X))^2) /(dim(X)-1)) ) $ Var(M) := Sum((M-Mean(M))^2) / dim(M)$ Var1(X) := Sum((X-Mean(X))^2) /(dim(X)-1.) $ sem(X) := float( std1(X)/ sqrt(length(X)) )$ mad(L) := Mean(map(lambda([x], abs(x - Mean(L))), L)); rms(X) := sqrt(apply("+", map(lambda([x], x^2), X)) / length(X))$ load ("descriptive")$ /* mode by M. R. Riotorto */ mode(x) := block([], freq : discrete_freq(x), maxfreq : sublist_indices(freq[2], lambda([z], z= lmax(freq[2]))), makelist(freq[1][i], i, maxfreq))$ moment(X,A,r) := float( Mean (map (lambda ([x], (x-A)^r),X) ))$ skew(X) := moment(X,3)/ sqrt(moment(X,2)^3) /* compatibel with MatLab */ $ moment(X,r) := float( Sum((X-Mean(X))^r)/length(X) )$ kurtosis(X) := moment(X,4)/moment(X,2)^2$ /*** DISTRIBUTIONs ***/ /* binomial distribution */ choose(n,x) := binomial(n,x)$ /* binoPDF(x,n,p) := binomial(n,x)*p^x*(1-p)^(n-x) $ */ binoPDF(k,n,p) := choose(n,k)*p^k*(1-p)^(n-k) $ binoCDF(k,n,p) := sum( choose(n,j)*p^j*(1-p)^(n-j) , j,0,k) $ /* geometric distribution */ geoPDF(k,p) := float( p*(1-p)^k )$ geoCDF(n,p) := float( 1-(1-p)^(n+1) ) /*-- explicit formula */ $ geoCDF1(n,p) := sum( geoPDF(k,p), k,0,n) /*-- math. definition */ $ geoINV(u,p) := ceiling(( log(1-u)/log(1-p) - 1)) $ /* negative binomial distribution */ nbinPDF(x,r,p) := choose(x+r-1,x)*p^r*(1-p)^x $ nbinPDF1(x,n,p) := gamma(n+x)/(x!*gamma(n)) *p^n*(1-p)^x $ nbinCDF(x,n,p) := beta_incomplete_regularized(n,floor(x)+1,p)$ nbinCDF1(x,n,p) := p^n*sum( choose(n+j-1,n-1)*(1-p)^j ,j,0,x) $ nbinINV(alpha,n,p) := block( [x,k:0], do( k : k+1, x : nbinCDF(k,n,p), if x>alpha then return(k)))$ /* hypergeometric distribution */ hgPDF(N,R,n,r) := float(choose(R,r)*choose(N-R, n-r)/choose(N,n))$ hgCDF(N,R,n, x) := sum( hgPDF(N,R,n,k), k,0, floor(x)) $ hgINV(alpha, N,R,n) := block( [x,k:0], do( k : k+1, x : hgCDF(N,R,n, k), if x>alpha then return(k)))$ load(distrib)$ hgPDF1(N,R , n,r) := pdf_hypergeometric(r,R,N-R,n) /* compatibility with Maxima */ $ hgPDF2(x,n1,n2,n) := pdf_hypergeometric(x,n1,n2,n) /* original of distrib */ $ /* POISSON distribution */ poissonPDF(k, lambda) := float( lambda ^k*1/k! * exp(-lambda) )$ poissonCDF(x, lambda) := sum( poissonPDF(k, lambda) ,k, 0, floor(x))$ poissonINV(alpha,lambda) := block( [x,k:0], do( k : k+1, x : poissonCDF(k, lambda), if x>alpha then return(k)))$ /* normal distribution */ normPDF(x,mu,sigma) := float( 1/sqrt(2.*%pi)/sigma*exp(-0.5*((x-mu)/sigma)^2) )$ normCDF(x,m,s) := float( 1/2 + erf((x-m)/(s*sqrt(2)))/2 )$ stdnormCDF(x) := normCDF(x, 0,1)$ Phi(x) := stdnormCDF(x)$ normINV(q,m,s) := float( m + sqrt(2)*s*inverse_erf(2*q-1) )$ stdnormINV(q) := normINV(q,0,1) $ /* exponential distribution */ expPDF(x, m) := if x<0 then 0 else float( m*exp(-m*x) )$ expCDF(x, m) := if x<0 then 0 else float( 1 - exp(-m*x) )$ expINV(p,m) := - log(1-p)/m $ /* STUDENT t distribution */ Beta(a,b) := gamma(a)*gamma(b)/gamma(a+b)$ tPDF(x,f) := float( 1/(sqrt(f)*Beta(1/2,f/2))*(1+ x^2/f)^(-0.5*(f+1)) )$ betainc(a,b,c) := beta_incomplete_regularized(a,b,c)$ tCDF(x,n) := float( (1+signum(x))/2 - signum(x) * betainc(n/2,1/2,n/(n+x^2)) / 2 )$ load(distrib)$ tINV(t,df) := quantile_student_t(t,df) /* use function from distrib */ $ /* SNEDECOR F distribution */ /* the version by M.R. Riotorto in 'distrib' */ fPDF(x,m,n) := gamma((m+n)/2)*(m/n)^(m/2)*x^(m/2-1)*(1+m*x/n)^(-(m+n)/2) / (gamma(m/2)*gamma(n/2)) * unit_step(x) $ fCDF1(t,df1,df2) := betainc( 0.5*df1, 0.5*df2, t/(t+df2/df1))$ load(distrib)$ fINV(q,m,n) := quantile_f(q,m,n)$ /* CHI^2 distribution */ chi2PDF(x, k) := 1/(2^(k/2)*gamma(k/2))*x^(k/2-1)*exp(-x/2.) $ GammaP(x,k) := x^k*exp(-x) * sum( x^m/gamma(k+m+1), m,0,100)$ chi2CDF(x, k) := if x<=0 or k<=0 then return(0) else float(GammaP(x*0.5, k*0.5 )) $ /* build-in by M.R. Riotorto */ load(distrib)$ chi2cdf(x,n) := cdf_gamma(x,n/2,2)$ chi2INV(q,n) := quantile_gamma(q,n/2,2) /* alias */ $ /* PARETO distribution */ paretoPDF(x,m,a) := if x