June 24, 2007

FP 做数学统计题

今天布置的家庭作业还真是烦人!做完《数学之友》上“统计案例”一章所有习题。我的妈呀,算什么卡方、线性回归,公式繁琐地要命,按计算器都能把人按死。还是让电脑帮我做吧。
不用什么 GNUplot、Mathematica,俺用 Haskell,自己动手,丰衣足食!

一. 独立性检验
卡方就是一个公式:
  
kaf a b c d =
  ((sqr ((a * d) - (b * c))) * (a + b + c + d)) /
    ((a + b) * (c + d) * (a + c) * (d + b))

-- sqr 函数是用来求一个数的平方的,再实现一下:
sqr n = n * n
  
这样就搞定一种题了!以后算卡方不要太简单,把要检测独立性的四个数据按顺序作为 kaf 函数的参数输入 Hugs 就搞定了!
Main> kaf 184 61 91 9
11.097807845216 :: Double

二. 线性回归分析
其实就是求一个 y = a + bx 中的 a,b。把 xn 值、yn 值各表示为一个列表,像这样来使用:
Main> lrb [1..8] [5.54, 7.52, 10.02, 11.73, 15.69, 16.12, 16.98, 21.06]
2.12142857142857 :: Double
Main> lra [1..8] [5.54, 7.52, 10.02, 11.73, 15.69, 16.12, 16.98, 21.06]
3.53607142857143 :: Double
那么,就可以根据公式把它们抄出来:
  
lrb xl yl =
  -- 这里要注意一下,数学中的 sigma 是一种记法而非函数
  (sigma
   (\(x, y) -> (x - xp) * (y - yp)) (zip xl yl)) /
    sigma (\x -> sqr (x - xp)) xl
  where
    xp = avg xl
    yp = avg yl

lra xl yl = avg yl - (lrb xl yl) * avg xl
  
其中,avg 函数用来求一个列表的平均值(这里也就是所有 x 取值的平均值啦),它需要知道列表长度:
  
avg ls = sum ls / len ls

len [] = 0
len (x:xl) = (len xl) + 1
  
最后一个 sigma 函数的实现也很简单,就是累加嘛,只是写地有些不好看罢了:
  
sigma f (x:xl) =
  if xl == []
  then f x
  else f x + sigma f xl
  
重申一遍,这里的 sigma 被表示为函数而非数学形式。那么,猜猜下面的 sigma 的数学形式是什么?
Main> sigma (\(x, y) -> x * y) [(1,4), (5,3), (3,8), (2,4)]
51 :: Integer

三. 线性相关系数
线性相关系数 r 的求法也只是个函数,再抄一遍:
  
lrr xl yl =
    (sigma
   (\(x, y) -> (x - xp) * (y - yp)) (zip xl yl)) /
    sqrt ((sigma (\x -> sqr (x - xp)) xl) *
     (sigma (\y -> sqr (y - yp)) yl))
  where
    xp = avg xl
    yp = avg yl
  
以后这样用就行了:
Main> lrr [1..8] [5.54, 7.52, 10.02, 11.73, 15.69, 16.12, 16.98, 21.06]
0.987345979074916 :: Double
哇!这组数据的线性相关性接近 1,好高啊!

四. 结语
现在发现 Haskell 还真是实用,这个周末数学作业不愁了,打魔兽去喽!
最后给楼下的看客留一个与 FP 没什么关系的小题目:
对于一组线性相关数据,往往要求它的所有产生数据—— a、b、r,还有线性回归函数,顺便还要用这个函数再求一下对未来的预期。我如果当真输入 3 遍数据然后手工计算预期岂不是很傻?!想想看怎么让我只输入一遍数据就得到所有的产生数据(设计一下数据结构而已)还有线性回归函数(其实只有这个才有那么一点点技术含量)。

No comments: