# regressione.py

# Coefficiente di correlazione, usando 20.31.
def corr (x,y): return cov(x,y)/(ds(x)*ds(y))

# Covarianza, usando 20.7.
def cov (x,y):
  n=len(x); s=0.0
  for i in xrange(n): s+=x[i]*y[i]
  return (s-n*media(x)*media(y))/(n-1)

# Deviazione standard di x.
def ds (x): return math.sqrt(cov(x,x))

def media (x):
  s=0.0; n=len(x)
  for xi in x: s+=xi
  return s/n

# Coefficienti della parabola di regressione.
def pareg (x,y):
  n=len(x); smx=smx2=smx3=smx4=smx2y=smxy=smy=0.0
  for i in xrange(n):
    xi=x[i]; yi=y[i]
    smx+=xi; smx2+=xi*xi; smx3+=xi*xi*xi; smx4+=xi*xi*xi*xi
    smx2y+=xi*xi*yi; smxy+=xi*yi; smy+=yi
  C1=[smx4,smx3,smx2]; C2=[smx3,smx2,smx]; C3=[smx2,smx,n]
  Cy=[smx2y,smxy,smy]; d=det3([C1,C2,C3])
  if not d: return None
  a=det3([Cy,C2,C3])/d; b=det3([C1,Cy,C3])/d; c=det3([C1,C2,Cy])/d
  return [a,b,c]

# Coefficienti della retta di regressione.
def rereg (x,y):
  n=len(x); smxy=smx=smy=smxx=0.0
  for i in xrange(n):
    xi=x[i]; yi=y[i]
    smxy+=xi*yi; smx+=xi; smy+=yi; smxx+=xi*xi
  a=(n*smxy-smx*smy)/(n*smxx-smx*smx)
  b=smy/n-a*smx/n; return [a,b]