macro calc_correlation (y(),x())
local n, sumx, sumy
local sumyy,sumxx,sumyx
n = first(y) * second(y)
sumy = sum(y)
sumx = sum(x)
sumyx = z_dotprod (y,x)
sumxx = z_dotprod (x,x)
sumyy = z_dotprod (y,y)
= (n * sumyx - sumy * sumx)/(sqrt((n * sumxx - sumx * sumx)*(n * sumyy - sumy * sumy)))
endmacro

macro calc_multiple_regression (y(), x())
local n, c, res

n = first(y) * second(y)
if first (x) = n then 
c = second (x)
else
c = first (x)
endif


res = multiple_regression (y,x,n,c)


= "finished"
endmacro

macro calc_pdf_binomial (no_successes, no_trials, prob_success)
= 100 * factorial (no_trials) / (factorial (no_successes) * factorial (no_trials - no_successes)) * prob_success ^ no_successes * (1- prob_success) ^ (no_trials - no_successes)
endmacro

macro calc_pdf_poisson (x,mean)
= 100 * exp (-mean) * mean ^ x / z_factorial (x)
endmacro

macro calc_pp_normal (p)
local const1,const2,const3,const4,const5,const6,q,weight,weight1,weight2,z
const1 = 2.515517
const2 = 0.802853
const3 = 0.010328
const4 = 1.432788
const5 = 0.189269
const6 = 0.001308

q = 0.5 - abs (p - 0.5)

weight = (-2 * ln (q)) ^ .5
weight1 = const1 + weight * (const2 + const3 * weight)
weight2 = 1 + weight * (const4 + weight * (const5 + const6 * weight))

z = weight - weight1/ weight2
z = z * sign (p - 0.5)

= z
endmacro

macro calc_pp_Students_t (p,df)
local acc,alpha,evar,evar2,fvar1,fvar2,tvar3,tvar4,w,z,t,res

res = angle (1)

acc = 0.0001
evar = 0.005
evar2 = evar + evar

alpha = 2 * p - 1

if df = 1 then 
t = tan (1.57079633 * alpha)
endif
 
if df = 2 then
t = alpha * (2 / ( 1 - alpha * alpha)) ^ .5
endif

if df > 2 then

z = calc_pp_normal (p)

w = z * (1 + 2 / (1 + 8 * df))
tvar3 = df * (exp (w * w / df) -1)
tvar3 = sign (z) * tvar3 ^ .5

aaaa:
fvar1 = df_Students_t ( tvar3, df) - p
fvar2 = df_Students_t ( tvar3 + evar, df)
fvar2 = fva2 - df_Students_t ( (tvar3 - evar), df)

fvar2 = fvar2 / evar2
tvar4 = tvar3 - fvar1 / fvar2

if abs(tvar4 - tvar3) > acc then 
tvar3 = tvar4
goto aaaa
endif

t = tvar4

endif
=t
endmacro


macro calc_regression (y(), x())
rem results() is matrix for output
local n, i, j, k, beta, alpha
local sumxy,sumxx,sumyy,sumres,variance,df
local sst, ssr, sse
local name, res

n = first(x) * second(x)

sumxy = (n * z_dotprod(y,x) - sum(y) * sum(x)) / n
sumxx = (n * z_dotprod(x,x) - sum(x) ^ 2) / n
sumyy = (n * z_dotprod(y,y) - sum(y) ^ 2) / n

beta = sumxy / sumxx
alpha = sum(y)/n - beta * sum(x)/n

sst = sumyy
ssr = beta * sumxy
sse = sst - ssr 

df = n-2
variance = sse / df

rem i = cell (results(0,0) >< results(0,0), "x")
rem j = cell (results(0,0) >< results(0,0), "y")

name = create()
res = selectsheet("name")
res = changewidth(a1><d10,250)
res = changewidth(f1><h10,200)

a1(0,0) = "Simple linear regression"

i = 1
j = 4
a1 (j-2,i-1) = "Coefficient"
a1 (j-2,i) = "Estimate"
a1 (j-2,i+1) = "Standard error"
a1 (j-2,i+2) = "T statistic"
a1 (j-1,i-1) = "Intercept = "
a1 (j,i-1) = "Gradient = "
a1 (j-1,i) = alpha
a1 (j,i) = beta
a1 (j-1,i+1) = sqrt (variance * (1/n + average(x)^2 / sumxx))
a1 (j,i+1) = sqrt (variance / sumxx)
a1 (j-1,i+2) = alpha / sqrt (variance * (1/n + average(x)^2 / sumxx))
a1 (j,i+2) = beta / sqrt (variance / sumxx)

a1 (j+9,i-1) = "R2 = "
a1 (j+9,i) = ssr / sst

a1 (j+2,i-1) = "Source"
a1 (j+2,i) =   "Sums of squares"
a1 (j+2,i+1) = "Deg. of freedom"
a1 (j+2,i+2) = "Mean squares"
a1 (j+3,i-1) = "Regression"
a1 (j+3,i) = ssr
a1 (j+3,i+1) = 1
a1 (j+3,i+2) = ssr / 1

a1 (j+4,i-1) = "Residual"
a1 (j+4,i) = sse
a1 (j+4,i+1) = df
a1 (j+4,i+2) = variance

a1 (j+5,i-1) = "Total"
a1 (j+5,i) = sst
a1 (j+5,i+1) = n-1

a1 (j+7,i-1) = "F(1," cat df cat ") = "
a1 (j+7,i) = (ssr / 1) / variance

i = 6
j = 4
a1 (j-2,i-1) = "Observed y"
a1 (j-2,i) = "Predicted y"
a1 (j-2,i+1) = "Residual"

for k = 0 to n-1
a1 (j+k-1,i-1) = y(k)
a1 (j+k-1,i) = alpha + beta * x(k)
a1 (j+k-1,i+1) = y(k) - (alpha + beta * x(k))
next k 

 = "finished"
endmacro

macro calc_t_1_sample (y(),mean_h0)
local m,data,res,obs

if first(y) > 1 and second(y) > 1 then
 data = "freq"
 if first(y) >= second (y) then
  m = first(y)
  obs = "row"
 else
  m = second(y)
  obs = "col"
 endif
else
 data = "raw"
 m = first(y) * second(y)
endif

res = z_t_single_sample (y,m,mean_h0,data,obs)

= "finished"
endmacro

macro calc_t_2_sample (y(),x())
local n,m,data,res,obs

if first(y) > 1 and second(y) > 1 then
 data = "freq"
 else
 data = "raw"
endif

if first(y) >= second (y) then
  m = first(y)
  n = first(x)
  obs = "row"
 else
  m = second(y)
  n = second(x)
  obs = "col"
endif

res = z_t_two_sample (y,x,m,n,data,obs)

= "finished"
endmacro

macro create_freq_d (y(), classes)
local n, lowerval, upperval, classwidth, a, w, i, j, k, l, m, data
local name,res
local f(classes),midpoint(classes)
n = first(y) * second(y)

data="null"
if first (y) = 1 then 
data = "col"
endif
if second (y) = 1 then
data = "row"
endif
if data = "null" then i = reporterror ("Error in macro", "Invalid parameter", 0x11)

lowerval= min (y)
upperval= max (y)

a= z_scale (lowerval, upperval, upperval-lowerval)
classwidth = z_extract (a, "i=")
lowerval = z_extract (a, "l=")
upperval = z_extract (a, "h=")

rem lowerval= input ("Lower end point", "lowerval", lowerval, "OK", "ABORT")
rem upperval= input ("Upper end point", "upperval", upperval, "OK", "ABORT")

rem classes = int ((upperval-lowerval)/ classwidth)

rem classwidth =  int ((upperval-lowerval)/ (classes-1))
rem classwidth =  input ("Width of classes", "classwidth", classwidth, "OK", "ABORT")

for i = 0 to classes-1
f(i) = 0
next i

for i = 0 to n-1 
if data = "col" then 
w = int ((y(0,i)- lowerval + classwidth) / classwidth) -1
else
w = int ((y(i)- lowerval + classwidth) / classwidth) -1
endif
f(w) = f(w) + 1
next i

name = create()
res = selectsheet("name")
res = changewidth(a1><a10,360)
res = changewidth(b1><b10,200)
res = changewidth(c1><c10,200)

a1(0,0) = "Frequency distribution"
a1(2,0) = "Class interval"
a1(2,1) = "Mid-point"
a1(2,2) = "Frequency"

for i=0 to classes-1
j = lowerval + classwidth * (i)
k = j + classwidth

if i=classes-1 then 
k=upperval
endif

midpoint(i) = j + .5 * (k - j)

a1(i+3,0) = j cat" <= x < " cat k
a1(i+3,1) = j + .5 * (k - j)
a1(i+3,2) = f(i)

next i

a1(classes+4,0) = "mean (freq data)"
a1(classes+4,1) = freq_d_mean(midpoint,f)
a1(classes+5,0) = "variance (freq data)"
a1(classes+5,1) = freq_d_variance(midpoint,f)

a1(classes+4,2) = average(y)
a1(classes+5,2) = stdev(y)^2

= "finished"
endmacro

macro df_Chisquare (X, df)
local a,g,kvar3,l,s,t,xvar3, p
xvar3 = 0.5 * X
kvar3 = 0.5 * df
l = z_lngamma_2 (kvar3 + 1)

s = 0
a = exp (kvar3 * ln (xvar3) - l - xvar3)
if a = 0 then goto aaaa
t = 1
s = 1
g = kvar3

bbbb:
g = g + 1
t = t * xvar3 / g
s = s + t
if t / s > 1.0E-6 then goto bbbb

aaaa:
p = s * a
= p
endmacro

macro df_Fishers_F (F, m, n)
local i,acc,hvar1,hvar2,hvar3,l,lvar1,p,tvar1,tvar2,wvar1,wvar2,wvar3,wvar4,x

acc = 0.001
if m = 1 then acc = 0.00001

hvar2 = n / 2
hvar1 = m / 2
hvar3 = hvar1 + hvar2

lvar1 = 0

lvar1 = lvar1 + z_lngamma_2 (hvar1)
lvar1 = lvar1 + z_lngamma_2 (hvar2)
lvar1 = lvar1 - z_lngamma_2 (hvar3)

wvar1 = n / (n + m * F)
x = 0
if hvar2 >= hvar3 * wvar1 then
wvar2 = 1 - wvar1
x = 1
else
wvar2 = wvar1
wvar1 = 1 - wvar1
wvar3 = hvar2
hvar2 = hvar1
hvar1 = wvar3
endif

tvar1 = 1
wvar3 = 1
p = 1
i = int (hvar1 + wvar2 * hvar3)
wvar4 = wvar1 / wvar2

xxxx:
tvar2 = hvar1 - wvar3
if i <> 0 then goto yyyy
wvar4 = wvar1

yyyy:
tvar1 = tvar1 * tvar2 * wvar4 / (hvar2 + wvar3)
p = p + tvar1
tvar2 = abs (tvar1)
if (tvar2 <= acc) and (tvar2 <= acc * p) then goto zzzz

wvar3 = wvar3 + 1
i = i -1
if i>= 0 then goto xxxx

tvar2 = hvar3
hvar3 = hvar3 + 1
goto yyyy

zzzz:
wvar1 = exp (hvar2 * ln (wvar1) + (hvar1 - 1) * ln (wvar2) - lvar1)
p = p * wvar1 / hvar2
if x = 0 then goto wwww
p = 1 - p

wwww:
= p
endmacro

macro df_normal (x, mean, variance)
local const1,const2,const3,const4,const5,weight,weight1, z, prob
const1 = 0.319381530
const2 =-0.356563782
const3 = 1.781477937
const4 =-1.821255978
const5 = 1.330274429

z= (x-mean) / (variance ^ .5)

weight1 = abs (z)
weight = 1 / (1 + 0.2316419 * weight1)
weight1 = 0.39894228 * exp(-.5 * weight1 * weight1)
prob = weight * (const1 + weight * (const2 + weight * (const3 + weight * (const4 + const5 * weight))))
prob = weight1 * prob
if z > 0 then prob = 1 - prob

= prob
endmacro

macro df_Students_t (t, df)
local pi, c, s, i, sumx, term, theta,res

res = angle (1)

pi = 3.14159265
term = df
theta = atan (t/ term^(1/2))
s = sin (theta)
c = cos (theta)
sumx = 0
if df > 1 then
  if df % 2 = 1 then
   i = 3
   term = c
   else
   i = 2
   term = 1
  endif
 sumx = term
 while i < df
  term = term * c * c * (i-1) / i
  sumx = sumx + term
  i = i + 2
 endwhile
 sumx = s * sumx
endif
if df % 2 = 1 then 
 sumx = 2/pi * (sumx + theta)
endif
= 0.5 * (1 + sumx)
endmacro


macro freq_d_mean (x(),f())
local n, i
n = first(x) * second(x)
= z_dotprod (x,f) / sum (f)
endmacro

macro freq_d_variance (x(),f())
local n
n = first(x) * second(x)
= (z_dotprod_y2x (x,f)  - sum (f) * freq_d_mean (x, f) ^ 2) / (sum (f) - 1)
endmacro

macro linear_solver (s(), c, b(), n_1)
local n
local name, res
local w(c)
local i,i_1,i_2,i_3,j,j_1,j_2,j_3,n_2,n_3,c_1
local error,v,t

rem inv(c)

error = 1.0E-6

c_1=c+1

rem factorisation
n_1=0
j_1=1
j_2=0
for j = 1 to c
 i_1 = 0
 for i = 1 to j
  j_2=j_2+1
  v=s(j_2-1)
  i_2=j_1
  for i_3 = 1 to i
   i_1 = i_1+1
   if i_3=i then goto xxxx
   v=v-s(i_1-1)*s(i_2-1)

   i_2=i_2+1
  next i_3

  xxxx:
  if i=j then goto yyyy
  if s(i_1-1)=0 then goto wwww
  s(j_2-1)=v/s(i_1-1)
  goto zzzz

  wwww:
  s(j_2-1)=0
  
 zzzz:
 next i

 yyyy:
 if abs(v)<error then goto aaaa
 if v<0 then res = reporterror ("Negative diagonal element", "S(" cat j_2 cat ")",0x11)
 s(j_2-1)=sqrt(v)
 goto bbbb

 aaaa:
 s(j_2-1)=0
 n_1=n_1+1
 res = reporterror ("Zero diagonal element", "S(" cat j_2 cat ")",0x11)

 bbbb:
 j_1=j_1+j
next j


rem inversion of s
n_2=c*c_1/2
i=c
n_3=n_2

iiii:
if s(n_3-1)=0 then goto cccc
i_1=n_3
for i_3=i to c
 w(i_3-1)=s(i_1-1)
 i_1=i_1+i_3
next i_3
j=c
j_3=n_2
i_2=n_2

hhhh:
i_1=j_3
v=0
if i<>j then goto dddd
v=1.0/w(i-1)

dddd:
j_2=c

ffff:
if j_2=i then goto eeee
v=v-w(j_2-1)*s(i_1-1)
j_2=j_2-1
i_1=i_1-1
if i_1<=i_2 then goto ffff
i_1=i_1-j_2+1
goto ffff

eeee:
s(i_1-1)=v/w(i-1)
if i=j then goto gggg
i_2=i_2-j
j=j-1
j_3=j_3-1
goto hhhh

cccc:
i_1=n_3
for j_1=i to c
 s(i_1-1)=0
 i_1=i_1+j_1
next j_1

gggg:
n_3=n_3-i
i=i-1
if i<>0 then goto iiii
 
rem form solution vector b()
for i=1 to c
 w(i-1)=s(n_2+i-1)
next i


i_1=0
for i=1 to c
 i_1=i_1+i-1
 t=0
 for j=1 to i
  i_2=i_1+j
  t=t+s(i_2-1)*w(j-1)
 next j
 if i=c then goto jjjj
 for j=i+1 to c
  i_2=i_2+j-1
  t=t+s(i_2-1)*w(j-1)
 next j

jjjj:
b(i-1)=t
next i


= "finished"
endmacro

macro multiple_regression (y(), z(), n, c)
local name, res
local i,i_1,j,k, c_1, n_1,n_2
local t,t_1,t_2,w,s_2,a,s_3, w_1, w_2,s_0,sea
local x(n,c+1), m(c+1), s((c+1)*(c+1+1)/2)
local f(n),r(n)
local b(c)

c_1=c+1

for i=1 to n
for j=1 to c
x(i-1,j-1)=z(i-1,j-1)
next j
next i

for i = 1 to n
x(i-1,c_1-1)=y(i-1)
next i

for j=1 to c_1
 t=0
 for i=1 to n
  t=t+x(i-1,j-1)
 next i
 m(j-1)=t/n
next j

for i=1 to c_1
 t_1=m(i-1)
 i_1=i*(i-1)/2
 for j=1 to i
  t=0
  t_2=m(j-1)
  for k=1 to n
   t=t+(x(k-1,i-1)-t_1)*(x(k-1,j-1)-t_2)
  next k
  s(i_1+j-1)=t
 next j
next i


res = linear_solver (s,c,b,n_1)

rem calculation of fitted values and residuals

w=0
for i=1 to c
 w=w+b(i-1)*m(i-1)
next i

a=m(c_1-1)-w
s_2=0
for i=1 to n
 t=a
 for j=1 to c
  t=t+b(j-1)*x(i-1,j-1)
 next j
 f(i-1)=t
 w=y(i-1)-t
 r(i-1)=w
 s_2=s_2+w*w
next i
s_3=s(c_1*(c_1+1)/2-1)


rem output multiple regression
rem input intercept a, regression coefficints b(), sum of squares s_2 and s_3
rem fitted values f(), residuals r(), dependent variables y()

name = create()
res = selectsheet("name")
res = changewidth(a1><d10,250)
res = changewidth(f1><h10,200)

a1(0,0) = "Multiple linear regression"

w_1=s_2/(n-c_1)

rem standard error of a
s_0=s(0)*m(0)*m(0)
if c>1 then
for i=2 to c
 for j=1 to i-1
  k= (i*(i-1))/2+j
  s_0=s_0+2*s(k-1)*m(i-1)*m(j-1)
 next j
 s_0=s_0+s(k+1-1)*m(i-1)*m(i-1)
next i
endif
sea = sqrt(w_1 * (s_0 + 1/n))

rem ------------

k = 3
j = 0
a1 (k-1,j) = "Coefficient"
a1 (k-1,j+1) = "Estimate"
a1 (k-1,j+2) = "Standard error"
a1 (k-1,j+3) = "T statistic"

a1 (k,j) = "a = "
a1 (k,j+1) = a
a1 (k,j+2) = sea

for i=1 to c
w_2=w_1*s(i*(i+1)/2-1)
a1 (k+i,j) = "b" cat i cat " ="
a1 (k+i,j+1) = b(i-1)
a1 (k+i,j+2) = sqrt (w_2)
a1 (k+i,j+3) = b(i-1)/sqrt (w_2)
next i

k=k+c+3

a1 (k-1,j) = "Source"
a1 (k-1,j+1) = "Sums of squares"
a1 (k-1,j+2) = "Deg. of freedom"
a1 (k-1,j+3) = "Mean squares"

w_2=(s_3-s_2)/c
n_2=n-c-1

a1 (k,j) = "Regression"
a1 (k,j+1) = s_3-s_2
a1 (k,j+2) = c
a1 (k,j+3) = w_2

a1 (k+1,j) = "Residual"
a1 (k+1,j+1) = s_2
a1 (k+1,j+2) = n_2
a1 (k+1,j+3) = w_1

a1 (k+2,j) = "Total"
a1 (k+2,j+1) = s_3
a1 (k+2,j+2) = n-1

a1 (k+4,j) = "F(" cat c cat "," cat n_2 cat ") = "
a1 (k+4,j+1) = w_2/w_1

a1 (k+6,j)= "R2 = "
a1 (k+6,j+1) = (s_3-s_2)/s_3

k=3
j=5

a1 (k-1,j) = "Observed y"
a1 (k-1,j+1) = "Predicted y"
a1 (k-1,j+2) = "Residual"

for i=1 to n
a1 (k+i-1,j) = y(i-1)
a1 (k+i-1,j+1) = f(i-1)
a1 (k+i-1,j+2) = r(i-1)
next i

= "finished"
endmacro

macro pdf_Fishers_F (F, m, n)
local pi, a, b, c, d, e
pi = 3.14159265
a = z_factorial ( (m+n-2) / 2 ) 
b = z_factorial ((m-2)/2) * z_factorial ((n-2)/2) 
c = (m/n)^(m/2) 
d = F^((m-2)/2)
e = (1 + m*F/n)^((m+n)/2)
= (a / b ) * c * (d / e)
endmacro

macro pdf_Fishers_F_1 (F, m, n)
local pi, a, b, c, d, e
pi = 3.14159265
a = gamma ((m+n)/2)
b = gamma (m/2) * gamma (n/2)
c = (m/n)^(m/2)
d = F ^ ((m-2)/2)
e = (1+(m/n)*F)^((m+n/2))
= (a / b) * c * (d / e)
endmacro

macro pdf_normal (x, mean, variance)
local standarddev, pi
pi = 3.14159265
standarddev = variance ^ .5
= 1 / (standarddev * (2 * pi) ^ .5) * exp (-.5 * ((x - mean) ^ 2) / variance)
endmacro

macro pdf_Students_t (t, df)
local pi, a, b, c
pi = 3.14159265
a =  z_gamma ((df+1)/2 ) 
b = (df * pi) ^ 0.5 * z_gamma (df/2) 
c = (1 + t*t/df) ^ ( (df+1)/2 )
= (a/b) * (1/c)
endmacro

macro pdf_Students_t_1 (t, df)
local pi, a, b, c
pi = 3.14159265
a = z_factorial ((df - 1)/2)
b = (df*pi)^(0.5) * z_factorial ((df-2)/2)
c = (1 + (t^2/df))^((df+1)/2)
= a/b*1/c
endmacro


macro stat_t_single_mean (mean, variance, sample_size, mean_h0)
local t,df,st_error,name, cv_5_1t,cv_5_2t,prob,res,cint

st_error = (variance/sample_size)^.5
t = (mean - mean_h0) / st_error
df = sample_size - 1
cv_5_2t = calc_pp_Students_t(.975,df)
cv_5_1t = calc_pp_Students_t(.95,df)
prob = 1- df_Students_t(t,df)

cint = cv_5_2t * st_error

name = create()
res = selectsheet ("name")
res = changewidth (a1><a10,750)

a1(0,0) = "T test (mean of population)"
a1(2,0) = "Mean under null hypothesis"
a1(2,1) = mean_h0
a1(3,0) = "Sample mean"
a1(3,1) = mean
a1(4,0) = "Standard error of mean"
a1(4,1) = st_error
a1(5,0) = "t-statistic"
a1(5,1) = t
a1(6,0) = "Degrees of freedom"
a1(6,1) = df
a1(7,0) = "Critical value for two-tail test (95% level)"
a1(7,1) = cv_5_2t
a1(8,0) = "Critical value for one-tail test (95% level)"
a1(8,1) = cv_5_1t
a1(9,0) = "Probability mean <> sample mean (two tail test)"
a1(9,1) = 2 * prob
a1(10,0) = "Probability mean <> sample mean (one tail test)"
a1(10,1) = prob

a1(11,0)= "95% central confidence interval: lower limit"
a1(11,1)= mean - cint
a1(12,0)= "95% central confidence interval: upper limit"
a1(12,1)= mean + cint

= "finished"
endmacro

macro stat_t_two_means (mean_x,variance_x,sample_size_x,mean_y,variance_y,sample_size_y)
local t,df,variance,st_error,cint,name,cv_5_1t,cv_5_2t,prob,res
local f,degf1,degf2

df = sample_size_x + sample_size_y -2

variance = ((sample_size_x -1) * variance_x + (sample_size_y -1) * variance_y) / df
st_error = (variance / sample_size_x + variance / sample_size_y) ^ .5
t = (mean_x - mean_y) / st_error

cv_5_2t = calc_pp_Students_t(.975,df)
cv_5_1t = calc_pp_Students_t(.95,df)
prob = 1- df_Students_t(t,df)

if variance_y >= variance_x then
 f = variance_y / variance_x
 degf1 = sample_size_y -1
 degf2 = sample_size_x -1
else
f = variance_x / variance_y
 degf1 = sample_size_x -1
 degf2 = sample_size_y -1
endif

cint = cv_5_2t * st_error

name = create()
res = selectsheet ("name")
res = changewidth (a1><a11,750)

a1(0,0) = "T test (means of unpaired samples)"
a1(2,0) = "Difference in two means under null hypothesis"
a1(2,1) = 0
a1(3,0) = "Difference in sample means"
a1(3,1) = mean_x - mean_y
a1(4,0) = "Estimated variance of x"
a1(4,1) = variance_x
a1(5,0) = "Estimated variance of y"
a1(5,1) = variance_y
a1(6,0) = "F test of variances: F (" cat degf1 cat "," cat degf2 cat ")"
a1(6,1) = f
a1(7,0) = "Standard error of difference in means"
a1(7,1) = st_error
a1(8,0) = "t-statistic"
a1(8,1) = t
a1(9,0) = "Degrees of freedom"
a1(9,1) = df
a1(10,0) = "Critical value for two-tail test (95% level)"
a1(10,1) = cv_5_2t
a1(11,0) = "Critical value for one-tail test (95% level)"
a1(11,1) = cv_5_1t
a1(12,0) = "Probability mean <> sample mean (two-tail test)"
a1(12,1) = 2 * prob
a1(13,0) = "Probability mean <> sample mean (one-tail test)"
a1(13,1) = prob
a1(14,0) = "95% confidence interval: lower limit"
a1(14,1) = (mean_x - mean_y) - cint
a1(15,0) = "95% confidence interval: upper limit"
a1(15,1) = (mean_x - mean_y) + cint
= "finished"
endmacro

macro s_df_normal (a, b, mean, variance)
local h, n, i, integral, lastint, zsum1, zsum2, zsum4, error
error = .01
n = 2
h = 0.5 * (b-a)
zsum1 = ( pdf_normal (a, mean, variance) + pdf_normal (b, mean, variance) )
zsum2 = 0
zsum4 = pdf_normal ( (a + h), mean, variance)
integral = h * (zsum1 + 4 * zsum4)
repeat
 lastint = integral
 n = n * 2
 h = 0.5 * h
 zsum2 = zsum2 + zsum4
 zsum4 = 0
 i = 1
  repeat
   zsum4 = zsum4 + pdf_normal (a+i*h, mean, variance)
   i = i + 2
  until i > n
 integral = h * (zsum1 + 2 * zsum2 + 4 * zsum4)
until abs(integral - lastint) < error
= integral / 3
endmacro

macro z_dotprod (y(),x())
local rows, cols, dot, j
rows = first(y)
cols = second (y)
dot = 0
if rows = 1 then
for j = 0 to cols-1
dot = dot + x(0,j) * y(0,j)
next j
endif
if cols = 1 then
for j = 0 to rows-1
dot = dot + x(j) * y(j)
next j
endif
= dot
endmacro

macro z_dotprod_y2x (y(),x())
local rows, cols, dot, i
rows = first(y)
cols =  second(y)
dot = 0
if rows = 1 then
for i = 0 to cols-1
dot = dot + x(0,i) * y(0,i) * y(0,i)
next i
endif
if cols = 1 then
for i = 0 to rows-1
dot = dot + x(i) * y(i) * y(i)
next i
endif
= dot
endmacro

macro z_extract (a,b)
local res
res = find(a,b)
=mid(a,res+2,100)+0
endmacro

macro z_factorial (a)
local factor
factor = 1
if a > 1 then 
factor =  a * z_factorial (a-1)
endif
= factor
endmacro

macro z_firstpart(x)
=x+0
endmacro

macro z_gamma (g)
local pi
pi = 3.14159265
if g = floor (g) then
= z_factorial (g-1)
else
if g - .5 = floor (g) then
if g >= 1.5 then 
= z_factorial (g-1) * .5 * pi^0.5
else
= pi^0.5
endif
else
= "unacceptable operand"
endif
endif
endmacro

macro z_interval(smallest,largest,classes)
local dump,power,diff,span,lowest,highest,interval,a,b
if smallest = largest  then ="0:0"
diff = abs(largest-smallest)
power=0
while(diff >= 10)
    power=power+1
    diff=diff/10
endwhile
while(diff < 1)
    power=power-1
    diff=diff*10
endwhile
if diff >= 8 then span=10
if diff < 8 & diff >= 7 then span = 8
if diff < 7 & diff >= 6 then span = 7
if diff < 6 & diff >= 5 then span = 6
if diff < 5 & diff >= 4 then span = 5
if diff < 4 & diff >= 3 then span = 4
if diff < 3 & diff >= 2 then span = 3
if diff < 2 & diff >= 1 then span = 2
if diff < 1 & diff >  0 then span = 1
span = span * 10^power
if power > 0 then
interval = int(span/classes)
else
interval = round((span/classes), abs(power)+1)
endif
lowest = floor(smallest/interval) * interval
highest = ceil(largest / interval) *interval
a= power
= "a=" cat a cat "b=" cat b cat "i=" cat interval cat "l=" cat lowest cat "h=" cat highest 
endmacro

macro z_lngamma (g)
rem  2 * g is an integer >0
local pi, sumx
pi = 3.14159265
sumx = 0
g = g-1
while g > 0
sumx = sumx + ln (g)
g = g-1
endwhile
if g < 0 then sumx = sumx + ln(pi ^ 0.5)
= sumx
endmacro

macro z_lngamma_2 (g)
local pi, cons1, cons2, cons3, cons4, cons5, cons6, cons7, gamm1, w, s
pi = 3.14159265
cons1=76.18009173
cons2=-86.50532033
cons3=24.01409822
cons4=-1.231739516
cons5=0.001208580
cons6=-0.000005364
cons7=2.506628275
gamm1=g-1.0
w=gamm1+5.5
w=(gamm1+0.5)*ln(w)-w
s=1.0+cons1/(gamm1+1)+cons2/(gamm1+2)+cons3/(gamm1+3)
s=s+cons4/(gamm1+4)+cons5/(gamm1+5)+cons6/(gamm1+6)
=w+ln(cons7*s)
endmacro

macro z_scale(smallest,largest,classes)
local dump,power,diff,span,lowest,highest,interval,a,b
if smallest = largest  then ="0:0"
diff = abs(largest-smallest)
power=0
while(diff >= 10)
    power=power+1
    diff=diff/10
endwhile
while(diff < 1)
    power=power-1
    diff=diff*10
endwhile
if diff >= 8 then span=10
if diff < 8 & diff >= 6 then span = 8
if diff < 7 & diff >= 5 then span = 7
if diff < 6 & diff >= 5 then span = 6
if diff < 5 & diff >= 4 then span = 5
if diff < 4 & diff >= 3 then span = 4
if diff < 3 & diff >= 2 then span = 3
if diff < 2 & diff >= 1 then span = 2
if diff < 1 & diff >  0 then span = 1
span = span * 10^power
if power>0 then
interval = int(span/classes)
else
interval = round((span/classes),abs(power)+1)
endif
lowest = floor(smallest/interval) * interval
highest = ceil(largest / interval) *interval
a=0
b=0
= "a=" cat a cat "b=" cat b cat "i=" cat interval cat "l=" cat lowest cat "h=" cat highest 
endmacro

macro z_secondpart(x)
local y
y=find(x,":")
=0+mid(x,y+1,100)
endmacro

macro z_Simpsons_int_f (a, b, df, dfden)
local h, n, i, integral, lastint, zsum1, zsum2, zsum4, error
error = .01
n = 2
h = 0.5 * (b-a)
zsum1 = ( f (a, df, dfden) + f (b, df, dfden) )
zsum2 = 0
zsum4 = f ( (a + h), df, dfden)
integral = h * (zsum1 + 4 * zsum4)
repeat
 lastint = integral
 n = n * 2
 h = 0.5 * h
 zsum2 = zsum2 + zsum4
 zsum4 = 0
 i = 1
  repeat
   zsum4 = zsum4 + f (a+i*h, df, dfden)
   i = i + 2
  until i > n
 integral = h * (zsum1 + 2 * zsum2 + 4 * zsum4)
until abs(integral - lastint) < error
= integral / 3
endmacro

macro z_t_single_sample (y(),m,mean_h0,data,obs)
local res,i,mean_y,variance_y,sample_size_y,mean_x,variance_x,sample_size_x
local midpty(m),fy(m)

if data = "raw" then
 mean_y = average(y)
 variance_y = estdev(y)^2
 sample_size_y = m
else
 if obs = "row" then
  for i = 0 to m-1
   midpty(i) = y(i,0)
   fy(i)   = y(i,1)
  next i
 else
  for i = 0 to m-1
   midpty(i) = y(0,i)
   fy(i)   = y(1,i)
  next i
 endif
 mean_y = freq_d_mean(midpty,fy)
 variance_y = freq_d_variance(midpty,fy)
 sample_size_y = sum(fy)
endif

res = stat_t_single_mean(mean_y,variance_y,sample_size_y,mean_h0)

= "finished"
endmacro

macro z_t_two_sample (y(),x(),m,n,data,obs)
local res,i,mean_y,variance_y,sample_size_y,mean_x,variance_x,sample_size_x
local midpty(m),fy(m),midptx(n),fx(n)

if data = "raw" then

  mean_y = average(y)
  mean_x = average(x)
  variance_y = estdev(y)^2
  variance_x = estdev(x)^2
  sample_size_y = m
  sample_size_x = n

else
 if obs = "row" then
  for i = 0 to m-1
   midpty(i) = y(i,0)
   fy(i)   = y(i,1)
  next i
 else
  for i = 0 to m-1
   midpty(i) = y(0,i)
   fy(i)   = y(1,i)
  next i
 endif
 mean_y = freq_d_mean(midpty,fy)
 variance_y = freq_d_variance(midpty,fy)
 sample_size_y = sum(fy)
 if obs = "row" then
  for i = 0 to n-1
   midptx(i) = x(i,0)
   fx(i)   = x(i,1)
  next i
 else
  for i = 0 to n-1
   midptx(i) = x(0,i)
   fx(i)   = x(1,i)
  next i
 endif
 mean_x = freq_d_mean(midptx,fx)
 variance_x = freq_d_variance(midptx,fx)
 sample_size_x = sum(fx)
endif

res = stat_t_two_means(mean_y,variance_y,sample_size_y,mean_x,variance_x,sample_size_x)

= "finished"
endmacro

macro z_t_two_sample_old (y(),x(),m,n,data,obs)
local res,i,mean_y,variance_y,sample_size_y,mean_x,variance_x,sample_size_x
local midpty(m),fy(m),midptx(n),fx(n)

if data = "raw" then
 if obs = "row" then
  mean_y = average(y)
  mean_x = average(x)
  variance_y = estdev(y)^2
  variance_x = estdev(x)^2
  sample_size_y = m
  sample_size_x = n
 else 

  for i = 0 to m-1
   midpty(i) = y(0,i)
   fy(i)   = 1
  next i
 
  mean_y = freq_d_mean(midpty,fy)
  variance_y = freq_d_variance(midpty,fy)
  sample_size_y = sum(fy)


  for i = 0 to n-1
   midptx(i) = x(0,i)
   fx(i)   = 1
  next i

  mean_x = freq_d_mean(midptx,fx)
  variance_x = freq_d_variance(midptx,fx)
  sample_size_x = sum(fx)

 endif 

else
 if obs = "row" then
  for i = 0 to m-1
   midpty(i) = y(i,0)
   fy(i)   = y(i,1)
  next i
 else
  for i = 0 to m-1
   midpty(i) = y(0,i)
   fy(i)   = y(1,i)
  next i
 endif
 mean_y = freq_d_mean(midpty,fy)
 variance_y = freq_d_variance(midpty,fy)
 sample_size_y = sum(fy)
 if obs = "row" then
  for i = 0 to n-1
   midptx(i) = x(i,0)
   fx(i)   = x(i,1)
  next i
 else
  for i = 0 to n-1
   midptx(i) = x(0,i)
   fx(i)   = x(1,i)
  next i
 endif
 mean_x = freq_d_mean(midptx,fx)
 variance_x = freq_d_variance(midptx,fx)
 sample_size_x = sum(fx)
endif

res = stat_t_two_means(mean_y,variance_y,sample_size_y,mean_x,variance_x,sample_size_x)

= "finished"
endmacro

