From f92699af9aa47eda9bd2c8a983487a706da75089 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 2 Oct 2015 10:39:21 +0200 Subject: pinv --- examples/pinv.hs | 13 +- examples/pinv.ipynb | 722 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 728 insertions(+), 7 deletions(-) create mode 100644 examples/pinv.ipynb diff --git a/examples/pinv.hs b/examples/pinv.hs index 7de50b8..6f093b4 100644 --- a/examples/pinv.hs +++ b/examples/pinv.hs @@ -1,20 +1,19 @@ import Numeric.LinearAlgebra -import Graphics.Plot import Text.Printf(printf) -expand :: Int -> Vector Double -> Matrix Double +expand :: Int -> Vector R -> Matrix R expand n x = fromColumns $ map (x^) [0 .. n] -polynomialModel :: Vector Double -> Vector Double -> Int - -> (Vector Double -> Vector Double) +polynomialModel :: Vector R -> Vector R -> Int + -> (Vector R -> Vector R) polynomialModel x y n = f where - f z = expand n z <> ws + f z = expand n z #> ws ws = expand n x <\> y main = do - [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt" + [x,y] <- toColumns <$> loadMatrix "data.txt" let pol = polynomialModel x y let view = [x, y, pol 1 x, pol 2 x, pol 3 x] putStrLn $ " x y p 1 p 2 p 3" putStrLn $ format " " (printf "%.2f") $ fromColumns view - mplot view + diff --git a/examples/pinv.ipynb b/examples/pinv.ipynb new file mode 100644 index 0000000..532b8d0 --- /dev/null +++ b/examples/pinv.ipynb @@ -0,0 +1,722 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import Numeric.LinearAlgebra" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import IHaskell.Display" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + ":ext FlexibleInstances\n", + "\n", + "dec = 3\n", + "\n", + "instance IHaskellDisplay (Matrix C) where\n", + " display m = return $ Display [html (\"

$$\"++(latexFormat \"bmatrix\" . dispcf dec) m++\"$$

\")]\n", + "\n", + "instance IHaskellDisplay (Matrix R) where\n", + " display m = return $ Display [html (\"

$$\"++ (latexFormat \"bmatrix\" . dispf dec) m++\"$$

\")]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import Graphics.SVG\n", + "data RawSVG = RawSVG String\n", + "instance IHaskellDisplay RawSVG where\n", + " display (RawSVG s) = return $ Display [html $ \"
\"++ s ++ \"
\"]\n", + "\n", + "lplot = RawSVG . hPlot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# least squares polynomial model" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "expand :: Int -> Vector R -> Matrix R\n", + "expand n x = fromColumns $ map (x^) [0 .. n]\n", + "\n", + "polynomialModel :: Vector R -> Vector R -> Int -> (Vector R -> Vector R)\n", + "polynomialModel x y n = f\n", + " where\n", + " f z = expand n z #> ws\n", + " ws = expand n x <\\> y" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "[x,y] <- toColumns <$> loadMatrix \"data.txt\"" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[0.9,2.1,3.1,4.0,4.9,6.1,7.0,7.9,9.1,10.2]" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "[1.1,3.9,9.2,51.8,25.3,35.7,49.4,3.6,81.5,99.5]" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x\n", + "y" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "

$$\\begin{bmatrix}\n", + "1.000 & 0.900 & 0.810\n", + "\\\\\n", + "1.000 & 2.100 & 4.410\n", + "\\\\\n", + "1.000 & 3.100 & 9.610\n", + "\\\\\n", + "1.000 & 4.000 & 16.000\n", + "\\\\\n", + "1.000 & 4.900 & 24.010\n", + "\\\\\n", + "1.000 & 6.100 & 37.210\n", + "\\\\\n", + "1.000 & 7.000 & 49.000\n", + "\\\\\n", + "1.000 & 7.900 & 62.410\n", + "\\\\\n", + "1.000 & 9.100 & 82.810\n", + "\\\\\n", + "1.000 & 10.200 & 104.040\n", + "\\end{bmatrix}$$

" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "expand 2 x" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "pol = polynomialModel x y\n", + "view = [x, y, pol 1 x, pol 2 x, pol 3 x]" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + " x y p 1 p 2 p 3" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + " 0.90 1.10 -3.41 7.70 -6.99\n", + " 2.10 3.90 6.83 9.80 15.97\n", + " 3.10 9.20 15.36 13.39 25.09\n", + " 4.00 51.80 23.04 18.05 28.22\n", + " 4.90 25.30 30.72 24.05 28.86\n", + " 6.10 35.70 40.96 34.16 29.68\n", + " 7.00 49.40 48.64 43.31 33.17\n", + " 7.90 3.60 56.32 53.82 41.60\n", + " 9.10 81.50 66.57 69.92 64.39\n", + "10.20 99.50 75.95 86.80 101.01" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import Text.Printf\n", + "\n", + "putStrLn \" x y p 1 p 2 p 3\"\n", + "putStrLn $ format \" \" (printf \"%.2f\") $ fromColumns view" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "t = linspace 100 (minElement x, maxElement x)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + " polynomial models \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " \n", + "\n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " 2 \n", + " 4 \n", + " 6 \n", + " 8 \n", + " 10 \n", + "\n", + " \n", + "\n", + "\n", + "\n", + "\n", + " -20 \n", + " 0 \n", + " 20 \n", + " 40 \n", + " 60 \n", + " 80 \n", + " 100 \n", + " 120 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " \n", + "\n", + " \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " data \n", + " degree 1 \n", + " degree 2 \n", + " degree 3 \n", + " degree 9 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "lplot\n", + " [ plotMark x y \"none\" 1 circles \"red\" 3 \"data\"\n", + " , plot t (pol 1 t) \"blue\" 1 \"degree 1\"\n", + " , plot t (pol 2 t) \"green\" 1 \"degree 2\"\n", + " , plot t (pol 3 t) \"brown\" 1 \"degree 3\"\n", + " , plot t (pol 9 t) \"gray\" 1 \"degree 9\"\n", + " , MarginX 0.05, Title \"polynomial models\", LegendPos 0.05 0.95, MaxY 120\n", + " ] " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Haskell", + "language": "haskell", + "name": "haskell" + }, + "language_info": { + "codemirror_mode": "ihaskell", + "file_extension": ".hs", + "name": "haskell", + "version": "7.10.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} -- cgit v1.2.3