diff --git a/Algo.hs b/Algo.hs new file mode 100644 --- /dev/null +++ b/Algo.hs @@ -0,0 +1,20 @@ +module Algo where + +import Data.List +import qualified Data.Map as M + +-- Group data from a list of (key,value) pairs +-- to a Map of key->[value] lists. +groupFirst :: (Ord k) => [(k, v)] -> M.Map k [v] +groupFirst = M.map reverse . groupFirst' + +-- Group data from a list of (key,value) pairs +-- to a Map of key->[value] lists. +-- The values are present in reverse order. +groupFirst' :: (Ord k) => [(k, v)] -> M.Map k [v] +groupFirst' = foldl' add M.empty + where + add map_ (key, value) = M.insertWith prepend key [value] map_ + prepend (x:_) xs = x:xs +-- prepend _ _ = impossible + diff --git a/Args.hs b/Args.hs new file mode 100644 --- /dev/null +++ b/Args.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE DeriveDataTypeable #-} + +module Args + ( HIrt(..) + , getArgs + ) where + +import Types + +import System.Console.CmdArgs + +data HIrt = Convert + | Estimate + { responses :: FilePath + , iTaskParams :: Maybe FilePath + , iThetas :: Maybe FilePath + , algorithm :: Algorithm + , maxRounds :: Int + , precision :: Double + , intPrec :: Double + , thetaStats :: [StatisticType] + , taskStats :: [StatisticType] + , oResponses :: Maybe FilePath + , oRespFormat :: RespFormat + , oTaskParams :: Maybe FilePath + , oTheta :: Maybe FilePath + , oBayesPlot :: Maybe FilePath + , oAllTaskParams :: Maybe String + , oAllTheta :: Maybe String + } deriving (Show, Data, Typeable) + +convert :: HIrt +convert = Convert {} + +estimate :: HIrt +estimate = Estimate + { responses = def &= argPos 0 &= typ "RESPONSES" + , iTaskParams = def &= typFile &= help "initial task params" + , iThetas = def &= typFile &= help "initial contestant abilities" + , algorithm = JML &= help "estimation algorithm" + , maxRounds = 40 &= name "n" &= help "stop if more than rounds (defaults to 40)" + , precision = 1e-8 &= help "stop if change is less than (defaults to 1e-8)" + , intPrec = 1e-9 &= help "calculating precision of algorithms (default 1e-9)" + , thetaStats = def &= help "statistics for contestants" + , taskStats = def &= help "statistics for tasks" + , oResponses = def &= typFile &= help "output parsed responses" + , oRespFormat = List &= help "output response format" + , oTaskParams = def &= typFile &= help "output file for final task param estimate" + , oTheta = def &= typFile &= help "output file for final contestant ability estimate" + , oBayesPlot = def &= typFile &= help "output coordinate list for aposteriori bayes probability" + , oAllTaskParams = def &= typ "TEMPLATE" + &= help "template for output of all round task parameters (params -> params000)" + , oAllTheta = def &= typ "TEMPLATE" + &= help "template for output of all round thetas (thetas -> thetas000)" + } &= auto + +getArgs :: IO HIrt +getArgs = cmdArgs $ modes [estimate, convert] diff --git a/Driver.hs b/Driver.hs new file mode 100644 --- /dev/null +++ b/Driver.hs @@ -0,0 +1,111 @@ +{-# LANGUAGE DisambiguateRecordFields #-} +{-# LANGUAGE NamedFieldPuns #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE RecordWildCards #-} +module Driver where + +import Args +import Engine +import Input +import Output +import Parse +import Types +import Irt +import Statistics + +import Control.Monad +import Data.List +import qualified Data.Text as T +import Data.Text.Lazy.Builder (Builder) +import qualified Data.Vector.Generic as V +import Text.Printf +import System.IO + +setupEngine :: HIrt -> EngineParams +setupEngine s = EngineParams + { algorithm = Args.algorithm s + , maxRounds = Args.maxRounds s + , precision = Args.precision s + , intPrec = Args.intPrec s + } + +statistic :: (StatisticType -> State -> IO Statistic) -> StatisticType -> State -> IO [Column] +statistic stat t s = return . toCols =<< stat t s + where + toCols (SingleStatistic xs) = [column 11 (T.pack . show $ t) . map fixed' $ xs] + toCols (ListStatistic names xss) = + zipWith (column 11) names . transpose . map (map fixed') $ xss + +readParams :: FilePath -> State -> IO State +readParams name state = + let taskNames = V.toList . tasks . Engine.responses $ state + in return . updateParams state =<< fromFile (taskParamsParser $ taskNames) name + +readThetas :: FilePath -> State -> IO State +readThetas name state = return . updateThetas state =<< fromFile thetasParser name + +formatResp :: RespFormat -> Responses -> Builder +formatResp List = tableBody . responseColumns +formatResp HeadedList = buildTable . responseColumns +formatResp Dot = formatDottedResponses + +printStats :: [a] -> Responses -> IO () +printStats resp responses = + printf "Read %d responses of %d subjects on %d tasks, %d (%4.1f%%) unanswered %d ignored\n" + nread ncont ntask unans punan (nread - nresp) + :: IO () + where + nread = length resp + ncont = V.length . contestants $ responses + ntask = V.length . tasks $ responses + nresp = V.length . respAll $ responses + unans = ncont * ntask - nresp + punan = fromIntegral unans / fromIntegral (ncont * ntask) * 100.0 :: Double + +calcBayes :: State -> ([Column], [Column]) +calcBayes State {..} = (bounds, values) + where + merge c ys = map (\(x,y) -> (c,x,y)) ys + (bs,vs) = unzip . map (bayes 0.95) $ groupByContestant responses thetas params + (cs,xs,ys) = unzip3 . concat $ zipWith merge (V.toList $ contestants responses) vs + (lbs,ubs) = unzip bs + bounds = [ numberColumn 12 8 "BayesLB" lbs + , numberColumn 12 8 "BayesUB" ubs ] + values = [ textColumn "contestant" cs + , numberColumn 12 8 "x" xs + , numberColumn 12 8 "p" ys ] + + + +run :: HIrt -> IO () +run Convert {..} = return () +run args@Estimate {..} = do + putStrLn "Reading..." + resp <- readResponses responses + let responses = responsesFromList resp + istate = Engine.init responses + engine = setupEngine args in do + printStats resp responses + maybeWriteFile oResponses (formatResp oRespFormat responses) + putStrLn "Reading initial task parameters..." + istate <- maybe return readParams iTaskParams $ istate + putStrLn "Reading initial contestant abilities..." + istate <- maybe return readThetas iThetas $ istate + putStrLn "Fitting IRT model..." + results <- runEngine engine istate + thetas results `seq` params results `seq` return () + putStrLn "Calculating bayes probabilities..." + let (bayesBounds, bayesValues) = calcBayes results in do + putStrLn "Writing bayes probability values..." + maybeWriteFile oBayesPlot . buildTable $ bayesValues + putStrLn "Writing task parameters..." + let taskBase = tableTaskParams . getTaskParamsList $ results in do + taskStats <- zipWithM (statistic taskStatistic) taskStats (repeat results) + writeDefFile stdout oTaskParams . buildTable $ taskBase ++ concat taskStats + putStrLn "Writing contestant parameters..." + let thetaBase = tableThetas . getThetasList $ results in do + thetaStats <- zipWithM (statistic thetaStatistic) thetaStats (repeat results) + maybeWriteFile oTheta . buildTable $ thetaBase ++ concat thetaStats ++ bayesBounds + +main :: IO () +main = getArgs >>= run diff --git a/Engine.hs b/Engine.hs new file mode 100644 --- /dev/null +++ b/Engine.hs @@ -0,0 +1,141 @@ +{-# LANGUAGE NamedFieldPuns #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE RecordWildCards #-} +module Engine where + +import Irt +import Statistics +import Types +import qualified Listable as L + +import Data.List +import qualified Data.Map as M +import qualified Data.Vector.Generic as V +import Data.VectorSpace +import Text.Printf + +data State = State + { responses :: Responses + , params :: TaskParams + , thetas :: Thetas + , logLikelihood :: Double + , dparam :: TaskParam + , dtheta :: Double + , updatedTheta :: Bool + , updatedParams :: Bool + } deriving (Eq, Ord, Show) + +updateOne :: State -> State +updateOne s@State {..} = s { logLikelihood = lh } + where + lh = totalLogLikelihood responses params thetas + +updateTwo :: State -> State -> State +updateTwo old new = updateOne $ new { dparam = L.fromList dp, dtheta = dt} + where + [dt] = diff (thetas old) (thetas new) + dp = diff (params old) (params new) + diff xs ys = (max' . trans . sub xs) ys + where + max' :: [[Double]] -> [Double] + max' = map (maximum . map abs) + trans = transpose . map (L.toList) . V.toList + sub = V.zipWith (^-^) + +getTaskParamsList :: State -> [(Task, TaskParam)] +getTaskParamsList State {..} = zip (V.toList $ tasks responses) (V.toList params) + +getThetasList :: State -> [(Contestant, Theta)] +getThetasList State {..} = zip (V.toList $ contestants responses) (V.toList thetas) + +infinity :: Double +infinity = 1.0/0.0 + +init :: Responses -> State +init (responses@Responses {..}) = State + { responses = responses + , params = V.replicate (V.length tasks) defaultTask + , thetas = V.replicate (V.length contestants) 0 + , logLikelihood = -infinity + , dparam = L.fromList $ repeat infinity + , dtheta = infinity + , updatedTheta = False + , updatedParams = False + } + +data EngineParams = EngineParams + { algorithm :: Algorithm + , maxRounds :: Int + , precision :: Double + , intPrec :: Double + } + +updateThetas :: State -> [(Contestant, Theta)] -> State +updateThetas s@State {..} xs = s + { thetas = V.fromList . map (ts M.!) . V.toList $ contestants responses + , updatedTheta = True + } + where + ts = M.fromList xs + +updateParams :: State -> [(Task, TaskParam)] -> State +updateParams s@State {..} xs = s + { params = V.fromList . map (ps M.!) . V.toList $ tasks responses + , updatedParams = True + } + where + ps = M.fromList xs + +enhanceTaskParams :: EngineParams -> State -> State +enhanceTaskParams EngineParams {..} s@State {..} = s + { params = estimateAB $ groupByTask responses thetas params + , updatedTheta = False + , updatedParams = True + } + +enhanceThetas :: EngineParams -> State -> State +enhanceThetas EngineParams {..} s@State {..} = s + { thetas = estimateTheta $ groupByContestant responses thetas params + , updatedTheta = True + , updatedParams = False + } + +oneJML :: EngineParams -> State -> State +oneJML engine s@State { updatedTheta } + | updatedTheta = enhanceTaskParams engine s + | otherwise = enhanceThetas engine s + +oneBFGS :: EngineParams -> State -> State +oneBFGS EngineParams {..} s@State {..} = s + { thetas = ts + , params = ps + , updatedTheta = True + , updatedParams = True + } + where + (ts,ps) = estimateBfgs precision responses thetas params + +oneRound :: EngineParams -> State -> State +oneRound e@EngineParams { algorithm = JML } = oneJML e +oneRound e@EngineParams { algorithm = LBFGSB } = oneBFGS e + +dparamC :: TaskParam -> Double +dparamC a = paramC (a ^+^ a) - paramC a + +runEngine :: EngineParams -> State -> IO State +runEngine engine@EngineParams {..} s = loop 1 [] s + where + loop n states state@State {..} + | n > maxRounds = return state + | magnitude dparam + dtheta < precision = return state + | otherwise = do + printf "L:%+14.8f " logLikelihood :: IO () + printf "dA:%11.8f dB:%11.8f dC:%11.8f " (paramA dparam) (paramB dparam) (dparamC dparam) :: IO () + printf "dT:%11.8f\n" dtheta :: IO () + loop (n+1) (state:states) $ updateTwo state . oneRound engine $ state + +thetaStatistic :: StatisticType -> State -> IO Statistic +thetaStatistic t State {..} = statTheta t $ groupByContestant responses thetas params + +taskStatistic :: StatisticType -> State -> IO Statistic +taskStatistic t State {..} = statTask t $ groupByTask responses thetas params diff --git a/Input.hs b/Input.hs new file mode 100644 --- /dev/null +++ b/Input.hs @@ -0,0 +1,38 @@ +{-# LANGUAGE OverloadedStrings #-} + +module Input where + +import Types +import Parse + +import Data.Char +import Data.Attoparsec.Text (parseOnly) +import qualified Data.Text as T +import qualified Data.Text.IO as T +import Text.CSV + +fromCSV :: CSV -> [(Contestant, Task, Points)] +fromCSV xs = [(c,t,r > 0) | (c,t,r) <- concatMap row xs'] + where + row (name:ys) = zipWith (resp name) tasks ys + resp name task res = (T.pack name, task, readRes res) + + hasNames = (t1 == "") || (t1 == "name") + where t1 = (filter (not . isSpace) . head . head) xs + xs' | hasNames = tail xs + | otherwise = xs + + tasks :: [T.Text] + tasks | hasNames = (map T.pack . tail . head) xs + | otherwise = map (T.pack . reverse . take 2 . (++ repeat '0') . reverse . show) [1::Int ..] + readRes x = read x :: Int + + +readResponses :: FilePath -> IO [(Contestant, Task, Points)] +readResponses name = do + s <- T.readFile name + case parseOnly responsesParser s of + Right r -> return r + Left _ -> case parseCSV name . T.unpack $ s of + Right csv -> return $ fromCSV csv + Left _ -> error $ "could not parse " ++ name diff --git a/Irt.hs b/Irt.hs new file mode 100644 --- /dev/null +++ b/Irt.hs @@ -0,0 +1,162 @@ +{-# LANGUAGE CPP #-} +module Irt + ( groupByTask + , groupByContestant + , estimateTheta + , estimateAB + , estimateBfgs + , totalLogLikelihood + , paramVSumMap + , thetaVSumMap + , estTheta + ) where + +import Types +import Likelihood + +import Control.Arrow +import Control.Monad.ST.Safe +import Data.List +import Data.Monoid +--import Text.Printf +import Numeric.GSL.Minimization +import qualified Numeric.Lbfgsb as B +import qualified Data.Vector.Generic as V +import qualified Data.Vector.Generic.Mutable as M +import qualified Data.Vector.Unboxed as UV +import qualified Data.Vector.Unboxed.Mutable as UM + +delta :: Double +delta = 1e-8 + +maxOn :: (Ord b) => (a -> b) -> a -> a -> a +maxOn f x y | f x < f y = y + | otherwise = x + +maximize1 :: (Double -> (Double, Double)) -> Double -> Double +maximize1 fg x1 = r + where + bounds = [(Just *** Just) thetaBound] + [r] = B.minimize 5 1e5 delta [x1] bounds fg' + fg' [x] = let (fx, dx) = fg x in (-fx, [-dx]) + + +instance Monoid Double where + mempty = 0 + mappend = (+) + +sumMap :: (Monoid b, V.Vector v a) => (a -> b) -> v a -> b +sumMap f = foldl' mappend mempty . map f . V.toList + +thetaVSumMap :: (Monoid a) => (Points -> TaskParam -> Theta -> a) -> ContestantData -> a +thetaVSumMap f (t,v) = sumMap (\(d,r) -> f r d t) $ v + +paramVSumMap :: (Monoid a) => (Points -> TaskParam -> Theta -> a) -> TaskData -> a +paramVSumMap f (d,v) = sumMap (\(t,r) -> f r d t) $ v + +estTheta :: Theta -> UV.Vector (TaskParam, Points) -> Theta +estTheta it v = uncurry bound thetaBound $ maximize1 fg it + where + fg t = thetaVSumMap (((pick.).). logL) (t,v) + pick (lh, _, dt) = (lh, dt) + +estAB :: TaskParam -> UV.Vector (Theta, Points) -> TaskParam +#if(PL3) +estAB = error "not supported, try --algo lbfgsb" +#else +estAB (ia,ib) v = if f [ia', ib'] >= f [na, nb] + then (ia', ib') + else (na, nb) + where + (ia', ib') = clampTaskParam (ia, ib) + (na, nb) = clampTaskParam (ra, rb) + ([ra,rb],_x) = minimize NMSimplex2 delta 100 [5,10] (negate . f) [ia,ib] + f [a,b] = bound (-1e100) (1e100) $ paramVSumMap logLikelihood ((a,b),v) +#endif + +address :: (UV.Unbox a, UV.Unbox b) => UV.Vector b -> UV.Vector (Int, a) -> UV.Vector (b, a) +address vals = V.map (first (vals V.!)) + +groupByContestant :: Responses -> Thetas -> TaskParams -> ContestantsData +groupByContestant resp thetas diffs = zip (V.toList thetas) . map (address diffs) . V.toList . respCont $ resp + +groupByTask :: Responses -> Thetas -> TaskParams -> TasksData +groupByTask resp thetas diffs = zip (V.toList diffs) . map (address thetas) . V.toList . respTask $ resp + + +estimateTheta :: ContestantsData -> Thetas +estimateTheta = V.fromList . map (uncurry estTheta) + +estimateAB :: TasksData -> TaskParams +estimateAB = V.fromList . map (uncurry estAB) + +logLs :: Responses -> Thetas -> TaskParams -> (Double, Thetas, TaskParams) +logLs responses thetas params = runST $ do + z <- zero "fail" + r <- V.foldM plus z (respAll responses) + freeze r + where + nconts = V.length thetas + ntasks = V.length params + zero _ = do + a <- M.replicate nconts 0 + b <- M.replicate ntasks mempty + return (0, a, b) + plus (l, dthetas, dparams) (cont, task, r) = do + add dthetas cont dt + add dparams task dd + return (l+l', dthetas, dparams) + where + (l', dd, dt) = logL r (params V.! task) (thetas V.! cont) + add v i x = M.read v i >>= M.write v i . (`mappend` x) + freeze (l, v1, v2) = do + v1' <- V.freeze v1 + v2' <- V.freeze v2 + return (l, v1', v2') + +estimateBfgs :: Double -> Responses -> Thetas -> TaskParams -> (Thetas, TaskParams) +estimateBfgs prec responses ithetas iparams = unpack . V.fromList $ minimize' + where + nconts = V.length ithetas + ntasks = V.length iparams + haveParams = map (uncurry (/=)) [aBound, bBound, cBound] + unpack = second unpackParams . unpackTheta + unpackTheta = V.splitAt nconts + unpackParams xs = V.zipWith3 from3PL as bs cs + where + split haveParam xs | haveParam = V.splitAt ntasks xs + | otherwise = (zs, xs) + zs = V.replicate ntasks 0 + (as, ys) = split (haveParams !! 0) xs + (bs, us) = split (haveParams !! 1) ys + (cs, _ ) = split (haveParams !! 2) us + canonize x = (paramA x, paramB x, paramC x) + pack ts abs = V.concat $ ts : [ xs | (True, xs) <- zip haveParams [as,bs,cs]] + where + (as,bs,cs) = V.unzip3 . V.map canonize $ abs + bounds = map (Just *** Just) . concat $ + replicate nconts thetaBound : + [ replicate ntasks x | (True, x) <- zip haveParams [aBound, bBound, cBound] ] + + minimize' = B.minimize 25 1e5 prec (V.toList $ pack ithetas iparams) bounds fg' + + fg' :: [Double] -> (Double, [Double]) + fg' xs = let (a, b) = fg . V.fromList $ xs in (-a, map negate . V.toList $ b) + fg :: UV.Vector Double -> (Double, UV.Vector Double) + fg xs = let (l, ts, ps) = uncurry (logLs responses) . unpack $ xs in (l, pack ts ps) + +bound :: Double -> Double -> Double -> Double +bound min' max' x | x <= min' = min' + | x >= max' = max' + | otherwise = x + +#if(!PL3) +clampTaskParam :: TaskParam -> TaskParam +clampTaskParam = (uncurry bound aBound) *** (uncurry bound bBound) +#endif + +totalLogLikelihood :: Responses -> TaskParams -> Thetas -> Double +totalLogLikelihood responses difficulties thetas = V.sum . V.map (likelihood . val) . respAll $ responses + where + val (c, t, r) = (thetas V.! c, difficulties V.! t, r) + likelihood (t, d, r) = logLikelihood r d t diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c)2012, Ivan Labáth + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Ivan Labáth nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Likelihood.hs b/Likelihood.hs new file mode 100644 --- /dev/null +++ b/Likelihood.hs @@ -0,0 +1,55 @@ +{-# LANGUAGE CPP #-} +module Likelihood +where + +import Types +import Numeric.Extras (log1p) + +-- result x c fx df/dx df/dc +lp3 :: Bool -> Double -> Double -> (Double, Double, Double) +lp3 True x c = let e = exp x + t1 = (c-1)/(e+1) + t2 = e/(1+c*e) + in (-log1p ((1-c)*e/(1+c*e)), t1*t2, t2) +lp3 False x c = let e = exp (-x) + in (-log1p (e*(1+c/e)/(1-c)), e/(e+1), 1/(c-1)) + +logL3 :: Bool -> ABC -> Theta -> (Double, (Double, Double, Double), Double) +logL3 r (a,b,c) t = let (p, (pda, pdb), pdt) = pos (a,b) t + (lh, dlh, dc) = lp3 r p c + in (lh, (dlh*pda, dlh*pdb, dc), dlh*pdt) + +-- a,b t f df/da df/db df/dt +pos :: AB -> Theta -> (Double, (Double, Double), Double) +pos (a,b) t = (a*(b-t), (b-t, a), -a) + +-- result x fx df/dx +lp2 :: Bool -> Double -> (Double, Double) +lp2 False x = let (f, df) = lp2 True (-x) in (f,-df) +lp2 True x = let e = exp x in + (negate . log1p $ e, -e/(e+1)) + +logL2 :: Bool -> AB -> Theta -> (Double, (Double, Double), Double) +logL2 r (a,b) t = let (p, (pda, pdb), pdt) = pos (a,b) t + (lh, dlh) = lp2 r p in + (lh, (dlh*pda, dlh*pdb), dlh*pdt) + +logL :: Points -> TaskParam -> Theta -> (Double, TaskParam, Double) +#if(PL3) +logL r d t = logL3 r (a,b,c) t +#else +logL r d t = let (lh, (da,db,_), dt) = logL3 r (a,b,c) t in (lh, (da,db), dt) +#endif + where + a = paramA d + b = paramB d + c = paramC d + +logLikelihood :: Points -> TaskParam -> Theta -> Double +logLikelihood r d t = let (lh, _, _) = logL r d t in lh + +logLikelihood't :: Points -> TaskParam -> Theta -> Double +logLikelihood't r d t = let (_, _, dt) = logL r d t in dt + +logLikelihood'ab :: Points -> TaskParam -> Theta -> TaskParam +logLikelihood'ab r d t = let (_, ds, _) = logL r d t in ds diff --git a/Listable.hs b/Listable.hs new file mode 100644 --- /dev/null +++ b/Listable.hs @@ -0,0 +1,21 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} + +module Listable +where + +class Listable c a where + toList :: c -> [a] + fromList :: [a] -> c + +instance Listable a a where + toList a = [a] + fromList (a:_) = a + +instance Listable (a,a) a where + toList (a,b) = [a,b] + fromList (a:b:_) = (a,b) + +instance Listable (a,a,a) a where + toList (a,b,c) = [a,b,c] + fromList (a:b:c:_) = (a,b,c) diff --git a/Main.hs b/Main.hs new file mode 100644 --- /dev/null +++ b/Main.hs @@ -0,0 +1,1 @@ +import Driver (main) diff --git a/Output.hs b/Output.hs new file mode 100644 --- /dev/null +++ b/Output.hs @@ -0,0 +1,137 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE RecordWildCards #-} +module Output where + +import Types + +import Control.Arrow (second) +import Data.Char +import Data.List +import Data.Monoid +import Data.Text.Format +import Data.Text.Lazy.Builder +import System.IO +import qualified Data.Text as T +import qualified Data.Text.Lazy.IO as L +import qualified Data.Vector.Generic as V + +fixed' :: Double -> Builder +fixed' = left 11 ' ' . fixed 8 + +escape :: T.Text -> T.Text +escape = T.map nowhite + where + nowhite x | isSpace x = '_' + | otherwise = x + +showPts :: Points -> Char +showPts False = '0' +showPts True = '1' + +type Name = T.Text + +data Column = Column + { name :: Name + , width :: Int + , list :: [Builder] + } + +class Columnable a where + toColumns :: [a] -> [Column] + +column :: Int -> Name -> [Builder] -> Column +column w n xs = Column + { name = n + , width = max w (T.length n) + , list = xs + } + +mhead :: (Monoid a) => [a] -> a +mhead [] = mempty +mhead (x:_) = x + +mtail :: (Monoid a) => [a] -> [a] +mtail [] = [] +mtail (_:xs) = xs + +buildLine :: [Int] -> [Builder] -> Builder +buildLine widths = mconcat . (++ [singleton '\n']) . intersperse (singleton ' ') + . zipWith (flip left ' ') widths + +tableHead :: [Column] -> Builder +tableHead cs = buildLine (map width cs) . map (fromText . name) $ cs + +tableBody :: [Column] -> Builder +tableBody cs = build' . map list $ cs + where + widths = map width cs + build' xs + | all null xs = mempty + | otherwise = (buildLine widths . map mhead $ xs) `mappend` (build' . map mtail $ xs) + +buildTable :: [Column] -> Builder +buildTable [] = singleton '\n' +buildTable cs = tableHead cs `mappend` tableBody cs + +nameColumn :: [T.Text] -> Column +nameColumn = textColumn "" + +textColumn :: Name -> [T.Text] -> Column +textColumn name xs = column (maximum $ 0 : map T.length xs) name (map fromText xs) + +numberColumn :: Int -> Int -> T.Text -> [Double] -> Column +numberColumn w p n = column w n . map (fixed p) + + +thetaColumn :: [Theta] -> Column +thetaColumn = column 12 "theta" . map fixed' + +instance Columnable AB where + toColumns abs = let (as, bs) = unzip abs + in zipWith (\n -> column 11 n . map fixed') ["a", "b"] [as, bs] + +instance Columnable ABC where + toColumns abcs = let (as, bs, cs) = unzip3 abcs + in zipWith (\n -> column 11 n . map fixed') ["a", "b", "c"] [as, bs, cs] + +tableThetas :: [(Contestant, Theta)] -> [Column] +tableThetas xs = let (cs, ts) = unzip xs + in [nameColumn cs, thetaColumn ts] + +tableTaskParams :: [(Task, TaskParam)] -> [Column] +tableTaskParams xs = let (ts, ds) = unzip xs + in nameColumn ts : toColumns ds + +formatTextResponse :: (Contestant, Task, Points) -> Builder +formatTextResponse (c,t,r) = build "{} {} {}" (escape c, escape t, showPts r) + +responseColumns :: Responses -> [Column] +responseColumns xs = let (c,t,r) = (unzip3 . responsesToList) xs + in [ textColumn "contestant" c + , textColumn "task" t + , column 1 "result" (map (singleton . showPts) r) ] + +unlinesB :: [Builder] -> Builder +unlinesB = mconcat . concatMap (:[singleton '\n']) + +formatDottedResponses :: Responses -> Builder +formatDottedResponses Responses {..} = build' respCont + where + n :: Int + n = V.length tasks + line = V.update (V.replicate n '.') . V.map (second showPts) + buildLn xs = fromString $ V.toList xs ++ "\n" + build' = mconcat . V.toList . V.map (buildLn . line) + +writeFile :: FilePath -> Builder -> IO () +writeFile path = L.writeFile path . toLazyText + +maybeWriteFile :: Maybe FilePath -> Builder -> IO () +maybeWriteFile (Just path) = L.writeFile path . toLazyText +maybeWriteFile _ = const . return $ () + +writeDefFile :: Handle -> Maybe FilePath -> Builder -> IO () +writeDefFile _ (Just path) = L.writeFile path . toLazyText +writeDefFile h Nothing = L.hPutStr h . toLazyText diff --git a/Parse.hs b/Parse.hs new file mode 100644 --- /dev/null +++ b/Parse.hs @@ -0,0 +1,99 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TupleSections #-} +module Parse where + +import Types + +import Control.Applicative +import Control.Monad +import Data.Attoparsec.Text +import Data.Char +import Data.List +import Data.Text (Text) +import qualified Data.Text as T +import qualified Data.Text.IO as T + +skipLine :: Parser () +skipLine = skipWhile (not . isEndOfLine) >> endOfLine + +spaceyEndOfLine :: Parser () +spaceyEndOfLine = skipWhile (\x -> isSpace x && (not . isEndOfLine) x) >> endOfLine + +class Parseable a where + names :: a -> [T.Text] + parser :: Parser a + +instance Parseable AB where + names _ = ["a", "b"] + parser = (,) <$> double <* skipSpace <*> double + +instance Parseable ABC where + names _ = ["a", "b", "c"] + parser = (,,) <$> double <* skipSpace <*> double <* skipSpace <*> double + +headable :: Parser () -> Parser a -> Parser a +headable header list = (header >> list) <|> list + +line :: [T.Text] -> Parser () +line xs = (sequence . map name $ xs) >> skipLine + where name n = skipSpace >> string n + +responsesParser :: Parser [(Contestant, Task, Points)] +responsesParser = headable header list <* endOfInput + where + header = line ["contestant", "task", "result"] + list = many $ do + skipSpace + c <- takeTill isSpace <* skipSpace + t <- takeTill isSpace <* skipSpace + r <- double <* endOfLine + return (c, t, r > 0) + +thetasParser :: Parser [(Contestant, Theta)] +thetasParser = headable header list <* endOfInput + where + header = line ["theta"] + list = many $ do + skipSpace + c <- takeTill isSpace <* skipSpace + t <- double <* skipLine + return (c, t) + +taskParamParser :: Parser (Task, TaskParam) +taskParamParser = do + t <- takeTill isSpace <* skipSpace + p <- parser + return (t,p) + +taskParamsParser :: [Task] -> Parser [(Task, TaskParam)] +taskParamsParser tasks = headable header list <* endOfInput + where + header = line $ names (undefined :: TaskParam) + list = matchList paramsParser <|> (matchList $ many $ (taskParamParser <* skipLine)) + matchList p = do + xs <- p + when (length xs /= length tasks) $ + fail $ "expecting " ++ (show $ length tasks) + ++ " tasks while " ++ (show $ length xs) ++ " found." + when ((sort $ map fst xs) /= (sort tasks)) $ + fail "task names do not match" + return xs + + paramsParser = do + ps <- many $ parser <* endOfLine + return $ zip (tasks ++ repeat "") ps + +twologParser :: Parser [TaskParam] +twologParser = do + _ <- manyTill (skipWhile (not . isEndOfLine) >> endOfLine) endOfLine + many $ parser <* spaceyEndOfLine + +justParse :: Parser a -> Text -> a +justParse p s = case parseOnly p s of + Left e -> error e + Right r -> r + +fromFile :: Parser a -> FilePath -> IO a +fromFile p f = fmap (justParse p) (T.readFile f) diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/Statistics.hs b/Statistics.hs new file mode 100644 --- /dev/null +++ b/Statistics.hs @@ -0,0 +1,136 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE BangPatterns #-} +module Statistics + ( statTheta + , statTask + , bayes + ) where + +import Irt +import Likelihood +import Types + +import Control.Arrow +import Control.Monad +import Control.Monad.ST.Safe +import Data.List +import qualified Data.Vector.Generic as V +import qualified Data.Vector.Unboxed as UV +import Statistics.Sample +import System.Random.MWC + +statTheta :: StatisticType -> ContestantsData -> IO Statistic +statTheta TaskCount xs = return . SingleStatistic . map (fromIntegral . V.length . snd) $ xs +statTheta SolvedProp xs = return . SingleStatistic . map (prop . snd) $ xs + where + prop = mean . V.map (ok . snd) + ok True = 1 + ok False = 0 +statTheta LogLikelihood xs = return . SingleStatistic . map (thetaVSumMap logLikelihood) $ xs +statTheta DLogLikelihood xs = return . ListStatistic ["ELogLikelihood", "DLogLikelihood"] . + map (thetaVSumMap edlogl) $ xs + where + edlogl _ d t = [elog, dlog] + where + lt = logLikelihood False d t + lf = logLikelihood True d t + elog = lt * exp lt + lf * exp lf + dlog = (lt-lf)^(2::Int) * (exp $ lt + lf) / (1+paramC d)^(3::Int) +statTheta FisherSEM xs = return . SingleStatistic . map fish $ xs + where + fish = (1 /) . sqrt . thetaVSumMap inf2 + inf2 _r d t = (-paramA d) * (paramA d) * logLikelihood't True d t + * logLikelihood't False d t +statTheta Bootstrap xss = withSystemRandom $ asGenST calc + where + calc g = do + seed <- save g + return $ ListStatistic ["BootstrapSEM", "BootstrapLb", "BootstrapUb"] + . map (bootstrap seed 100000) $ xss + + stat :: [Double] -> (Double, Double, Double) + stat xs = (sem, lb, ub) + where + n = length xs + sem = sqrt . varianceUnbiased . UV.fromList $ xs + confidence = 0.95 :: Rational + dist = (n - round (fromIntegral n * confidence)) `div` 2 + select = (!! dist) + (lb, ub) = (select &&& (select . reverse)) . sort $ xs + + + bootstrap seed n xs = [err,lb,ub] + where + (err,lb,ub) = runST $ do + g <- restore seed + v <- bootstrap' g n xs + return $ stat v + + bootstrap' g n (t,xs) = replicateM n $ do + xs' <- resample g xs + return $! estTheta t xs' + + resample g xs = return . V.map (xs V.!) =<< V.replicateM l (uniformR (0,l-1) g) + where l = V.length xs + + +statTask :: StatisticType -> TasksData -> IO Statistic +statTask TaskCount xs = return . SingleStatistic . map (fromIntegral . V.length . snd) $ xs +statTask SolvedProp xs = return . SingleStatistic . map (prop . snd) $ xs + where + prop = mean . V.map (ok . snd) + ok True = 1 + ok False = 0 +statTask LogLikelihood xs = return . SingleStatistic . map (paramVSumMap logLikelihood) $ xs +statTask DLogLikelihood xs = return . ListStatistic ["ELogLikelihood", "DLogLikelihood"] . + map (paramVSumMap edlogl) $ xs + where + edlogl _ d t = [elog, dlog] + where + lt = logLikelihood False d t + lf = logLikelihood True d t + elog = lt * exp lt + lf * exp lf + dlog = (lt-lf)^(2::Int) * (exp $ lt + lf) / (1+paramC d)^(3::Int) + +bayes :: Double -> ContestantData -> ((Double, Double), [(Double, Double)]) +bayes confidence (theta, tasks) = ((lbound, ubound), distribution) + where + pick (lh, _, dt) = (p, p * dt) + where + p = exp lh + fg t = pick $ thetaVSumMap logL (t,tasks) + maxp = fst . fg $ theta + + aintegral = integrate $ evaluate 0.3 (0.1*maxp) (fst thetaBound) (snd thetaBound) + points = evaluate 0.3 (0.002*aintegral) (fst thetaBound) (snd thetaBound) + integral = integrate points + distribution :: [(Double, Double)] + distribution = map (second (/integral)) points + percentile = (1 - confidence) / 2 + lbound = bound percentile distribution + ubound = bound percentile (reverse distribution) + + evaluate maxx maxv x0 x1 = evaluate' l1 r1 ++ [second fst r1] + where + l1 = (x0, fg x0) + r1 = (x1, fg x1) + split l r x = evaluate' l m ++ evaluate' m r + where m = (x, fg x) + evaluate' l@(xl, (vl, dl)) r@(xr, (vr, dr)) + | xr - xl > maxx * 1.01 = split l r (xl + maxx) + | (xr - xl) * (abs (dv - dl) + abs (dv - dr)) > maxv = split l r ((xl + xr) / 2) + | otherwise = [(xl, vl)] + where + dv = (vr - vl) / (xr - xl) + + integrate :: [(Double, Double)] -> Double + integrate = integrate' 0 + where + integrate' !c ((xl,vl):next@((xr,vr):_)) = integrate' (c + (xr - xl) * (vl + vr) / 2) next + integrate' !c _ = c + + bound !c ((xl, vl):next@((xr,vr):_)) + | c > v = bound (c-v) next + | otherwise = xl + (xr - xl) * c / v + where + v = abs $ (xr - xl) * (vl + vr) / 2 diff --git a/Types.hs b/Types.hs new file mode 100644 --- /dev/null +++ b/Types.hs @@ -0,0 +1,117 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Types where + +import Algo + +import Control.Arrow +import Data.Data +import Data.List +import Data.Text (Text) +import qualified Data.Map as M +import qualified Data.Vector.Generic as V +import qualified Data.Vector as BV +import qualified Data.Vector.Unboxed as UV + +type Task = Text +type Contestant = Text + +type Theta = Double +type AB = (Double, Double) +type ABC = (Double, Double, Double) +#if(PL3) +type TaskParam = ABC +#else +type TaskParam = AB +#endif + +type ContestantData = (Theta, UV.Vector (TaskParam, Points)) +type ContestantsData = [ContestantData] +type TaskData = (TaskParam, UV.Vector (Theta, Points)) +type TasksData = [TaskData] + + +thetaBound :: (Double, Double) +thetaBound = (-10, 10) + +defaultTask :: TaskParam +paramA :: TaskParam -> Double +paramB :: TaskParam -> Double +paramC :: TaskParam -> Double +from3PL :: Double -> Double -> Double -> TaskParam +aBound :: (Double, Double) +bBound :: (Double, Double) +cBound :: (Double, Double) + +aBound = (-3, 5) +bBound = (-10, 10) +#if(PL3) +defaultTask = (1,0,0) +paramA (a,_,_) = a +paramB (_,b,_) = b +paramC (_,_,c) = c +from3PL a b c = (a,b,c) + +cBound = (0, 0.99) +#else +defaultTask = (1,0) +paramA (a,_) = a +paramB (_,b) = b +paramC _ = 0 +from3PL a b 0 = (a,b) +from3PL _ _ _ = error "compiled for 2PL model" +cBound = (0, 0) +#endif + +type Points = Bool + +type TaskParams = UV.Vector TaskParam +type Thetas = UV.Vector Theta + +data StatisticType = TaskCount | SolvedProp | LogLikelihood + | DLogLikelihood | FisherSEM | Bootstrap + deriving (Eq, Show, Data, Typeable) + +data Statistic = SingleStatistic [Double] + | ListStatistic [Text] [[Double]] + +data Algorithm = JML | LBFGSB + deriving (Eq, Show, Data, Typeable) + +data RespFormat = List | HeadedList | Dot + deriving (Eq, Show, Data, Typeable) + +data Responses = Responses + { tasks :: BV.Vector Task + , contestants :: BV.Vector Contestant + , respAll :: UV.Vector (Int, Int, Points) -- all responses (contestant, task, result) + , respCont :: BV.Vector (UV.Vector (Int, Points)) -- responses arranged by contestant + , respTask :: BV.Vector (UV.Vector (Int, Points)) -- responses arranged by task + } deriving (Eq, Ord, Show) + +responsesFromList :: [(Contestant, Task, Points)] -> Responses +responsesFromList list = Responses tasks contestants resp respCont respTask + where + sorted = sort list + byCont = groupFirst' . map (\(c, t, r) -> (c, (t, r))) $ sorted + byTask = groupFirst' . map (\(c, t, r) -> (t, (c, r))) $ sorted + + contestants = V.fromList . M.keys $ byCont + tasks = V.fromList . M.keys $ byTask + + indexes = M.fromList . map (\(a,b) -> (b,a)) . V.toList . V.indexed + contIndex = indexes contestants + taskIndex = indexes tasks + + index2 (c, t, r) = (contIndex M.! c, taskIndex M.! t, r) + resp = V.fromList . map index2 $ sorted + + index inds = V.fromList . M.elems . M.map (V.fromList . map (first (inds M.!))) + respCont = index taskIndex byCont + respTask = index contIndex byTask + +responsesToList :: Responses -> [(Contestant, Task, Points)] +responsesToList Responses {..} = map f . V.toList $ respAll + where + f (c, t, r) = (contestants V.! c, tasks V.! t, r) diff --git a/hirt.cabal b/hirt.cabal new file mode 100644 --- /dev/null +++ b/hirt.cabal @@ -0,0 +1,85 @@ +Name: hirt +Version: 0.0.1.0 +Synopsis: Calculates IRT 2PL and 3PL models +Description: + Program for fitting Item Response Theory (IRT) two (2PL) and + three (3PL) parameter logistic models. + . + Implements Joint Maximum Likelihood (JML) algorithm + (currently only supported in 2PL model) and via + generic function optimization using L-BFGS-B (both 2PL and 3PL). + . + Calculates item parameter and subject ability estimates and + log likelihood statistics. For contestant abilities supports + error estimates via Fisher information, and via two algorithms + of the author, namely bootstrap and Bayes a posteriori probability. + . + Supports outputting coordinate list for a plot of Bayes + a posteriori probability of individual subject abilities. + . + Part of a masters thesis of the author "http://people.ksp.sk/~ivan/irt/ebook.pdf" . + . + As a side note, it is currently a lacking proper documentation and user friendliness. + +Homepage: https://people.ksp.sk/~ivan/hirt +License: BSD3 +License-file: LICENSE +Author: Ivan Labáth +Maintainer: ivan@hirt.ksp.sk +Category: Math +Build-type: Simple +Cabal-version: >=1.6 + +Flag PL3 + Description: Compile for 3PL model, doesn't support JML yet. + Model needs to be selected at compile time. + Default: False + +-- Extra-source-files: + +Executable hirt + Main-is: Main.hs + + if flag(PL3) + cpp-options: -DPL3 + + Build-depends: base >= 4 && < 5, + vector >= 0.9 && < 0.10, + containers >= 0.4 && < 0.5, + text >= 0.11.1.13 && < 0.12, + attoparsec >= 0.10.1 && < 0.11, + text-format >= 0.3.0.7 && < 0.4, + csv >= 0.1.2 && < 0.2, + hmatrix >= 0.13.1.0 && < 0.14, + numeric-extras >= 0.0.2.2 && < 0.1, + cmdargs >= 0.9.3 && < 0.10, + random >= 1.0.1.1 && < 1.1, + statistics >= 0.10 && < 0.11, + mwc-random >= 0.12 && < 0.13, + vector-space >= 0.8 && < 0.9, + hlbfgsb >= 0.0.1.0 && < 0.1.0.0 + + Other-modules: + Algo + Args + Driver + Engine + Input + Irt + Likelihood + Listable + Main + Output + Parse + Setup + Statistics + Types + +source-repository head + type: darcs + location: http://people.ksp.sk/~ivan/hirt + +source-repository this + type: darcs + location: http://people.ksp.sk/~ivan/hirt + tag: 0.0.1.0