@@ -14,9 +14,9 @@ h2. Overview and motivation _Incanter is a <a href="http://clojure.org">Clojure</a>-based, <a href="http://www.r-project.org/">R</a>-like statistical computing and graphics environment for the JVM._ At the core of Incanter are the <a href="http://sites.google.com/site/piotrwendykier/software/parallelcolt">Parallel Colt</a> numerics library, a multithreaded version of <a href="http://acs.lbl.gov/~hoschek/colt/">Colt</a>, and the <a href="http://www.jfree.org/jfreechart/">JFreeChart</a> charting library, as well as several other Java and Clojure libraries. -The motivation for creating Incanter is to provide a JVM-based statistical computing and graphics platform with R-like semantics and interactive-programming environment. Running on the JVM provides access to the large number of existing Java libraries for data access, data processing, and presentation. Clojure’s seamless integration with Java makes leveraging these libraries much simpler than is possible in R, and Incanter’s R-like semantics makes statistical programming much simpler than is possible in pure Java. +The motivation for creating Incanter is to provide a JVM-based statistical computing and graphics platform with R-like semantics and interactive-programming environment. Running on the JVM provides access to the large number of existing Java libraries for data access, data processing, and presentation. Clojure’s seamless integration with Java makes leveraging these libraries much simpler than is possible in R, and Incanter’s R-like semantics makes statistical programming much simpler than is possible in pure Java. -Motivation for a Lisp-based R-like statistical environment can be found in the paper <i><a href="http://books.google.com/books?id=8Cf16JkKz30C&pg=PA21">Back to the Future: Lisp as a Base for a Statistical Computing System by Ihaka and Lang (2008)</a></i>. Incanter is also inspired by the now dormant Lisp-Stat (see the special volume in the Journal of Statistical Software on <a href="http://www.jstatsoft.org/v13">Lisp-Stat: Past, Present, and Future</a> from 2005). +Motivation for a Lisp-based R-like statistical environment can be found in the paper <i><a href="http://books.google.com/books?id=8Cf16JkKz30C&pg=PA21">Back to the Future: Lisp as a Base for a Statistical Computing System by Ihaka and Lang (2008)</a></i>. Incanter is also inspired by the now dormant Lisp-Stat (see the special volume in the Journal of Statistical Software on <a href="http://www.jstatsoft.org/v13">Lisp-Stat: Past, Present, and Future</a> from 2005). Motivation for a JVM-based Lisp can be found at the <i><a href="http://clojure.org">Clojure website</a></u></i>, and screencasts of several excellent Clojure talks by the language’s creator, Rich Hickey, can be found at <i><a href="http://clojure.blip.tv/">clojure.blip.tv</a></i>. @@ -47,7 +47,7 @@ From the Clojure shell, load the Incanter libraries: <pre><code>(use '(incanter Try an example: sample 1,000 values from a standard-normal distribution and view a histogram: <pre><code>(view (histogram (sample-normal 1000)))</code></pre> <img src="http://incanter.org/images/examples/normal_hist.png" /> -The online documentation for most Incanter functions contain usage examples. The documentation can be viewed using Clojure's @doc@ function. For example, to view the documentation and usage examples for the @linear-model@ function, call @(doc linear-model)@ from the Clojure shell. Use @(find-doc "search term")@ to search the online documentation from the Clojure shell. The API documentation can also be found in the @docs/api/@ directory and online at <a href="http://incanter.org/docs/api">http://incanter.org/docs/api/</a>. +The online documentation for most Incanter functions contain usage examples. The documentation can be viewed using Clojure's @doc@ function. For example, to view the documentation and usage examples for the @linear-model@ function, call @(doc linear-model)@ from the Clojure shell. Use @(find-doc "search term")@ to search the online documentation from the Clojure shell. The API documentation can also be found in the @docs/api/@ directory and online at <a href="http://incanter.org/docs/api">http://incanter.org/docs/api/</a>. *More Incanter examples* @@ -87,7 +87,7 @@ h2. Building Incanter The Incanter distribution includes a pre-built jar file, but if you would like to build it, test it, and generate the api documentation yourself, run the following commands from the Incanter directory: * Build incanter.jar: @ant@ -* Run unit tests: @tests/runtests.clj@ +* Run unit tests: @tests/runtests.clj@ * Generate API docs: @bin/generate_docs.clj@ (Note: the @bin/clj@ and @*.clj@ scripts won't run on Windows systems as is. Although, the Clojure code contained in the @*.clj@ scripts can be run from the Clojure shell.) @@ -109,6 +109,3 @@ See @lib/license.txt@ file for license terms of included libraries. * "JGraph":http://www.jgraph.com/ (for future development) * "Clojure-Contrib":http://code.google.com/p/clojure-contrib/ (test-is, gen-html-docs, repl-utils) * "JLine":http://jline.sourceforge.net/ - - - README.textile @@ -13,7 +13,7 @@ ;;; Commentary: ;; The `incanter' function starts or switches to a SLIME session with -;; Incanter loaded. +;; Incanter loaded. ;; You'll need clojure-mode and SLIME installed. Install clojure-mode ;; manually or through ELPA (http://tromey.com/elpa) and then hit M-x @@ -45,7 +45,7 @@ (slime)))) (defun incanter-init () - (save-excursion + (save-excursion (switch-to-buffer "*slime-repl clojure*") (rename-buffer "*incanter*") (slime-eval-async `(swank:eval-and-grab-output bin/incanter.el @@ -1,10 +1,10 @@ /* * Matrix.java -- Matrix class for Clojure that extends DenseColDoubleMatrix2D * in the CERN Colt Library -* +* * by David Edgar Liebke http://incanter.org * March 11, 2009 -* +* * Copyright (c) David Edgar Liebke, 2009. All rights reserved. The use * and distribution terms for this software are covered by the Eclipse * Public License 1.0 (http://opensource.org/licenses/eclipse-1.0.php) @@ -105,11 +105,11 @@ public class Matrix extends DenseColDoubleMatrix2D implements ISeq, Counted { return(mat); } - + public Matrix like(int rows, int columns) { - return new Matrix(rows, columns); + return new Matrix(rows, columns); } - + /************************************** * ISeq METHODS **************************************/ @@ -131,7 +131,7 @@ public class Matrix extends DenseColDoubleMatrix2D implements ISeq, Counted { public ISeq next() { - if(this.rows == 0 || this.columns == 0) + if(this.rows == 0 || this.columns == 0) return(null); else if(!this.oneDimensional && this.rows == 1) { return(null); @@ -203,33 +203,33 @@ public class Matrix extends DenseColDoubleMatrix2D implements ISeq, Counted { return(null); } - public int count() { - if(this.oneDimensional) + public int count() { + if(this.oneDimensional) return this.elements.length; else - return this.rows; + return this.rows; } - public Matrix seq() { + public Matrix seq() { if(this.elements.length > 0) return this; else - return null; + return null; } public IPersistentCollection empty() { - return(new Matrix(0, 0)); + return(new Matrix(0, 0)); } public boolean equiv(Object o) { return equals(o); - + } - + public String toString() { StringBuffer buf = new StringBuffer(); for(int i = 0; i < this.rows; i++) { - for(int j = 0; j < this.columns; j++) { + for(int j = 0; j < this.columns; j++) { buf.append(this.get(i, j)); buf.append(" "); } @@ -238,6 +238,3 @@ public class Matrix extends DenseColDoubleMatrix2D implements ISeq, Counted { return(buf.toString()); } } - - - src/incanter/Matrix.java @@ -16,11 +16,11 @@ -(ns incanter.bayes - (:use [incanter.core :only (matrix mmult mult div minus trans ncol nrow +(ns incanter.bayes + (:use [incanter.core :only (matrix mmult mult div minus trans ncol nrow plus to-list decomp-cholesky solve half-vectorize vectorize symmetric-matrix identity-matrix kronecker - bind-columns)] + bind-columns)] [incanter.stats :only (sample-normal sample-gamma sample-dirichlet sample-inv-wishart sample-mvn mean)])) @@ -28,9 +28,9 @@ -(defn sample-model-params +(defn sample-model-params " Returns a sample of the given size of the the parameters (coefficients and - error variance) of the given linear-model. The sample is generated using + error variance) of the given linear-model. The sample is generated using Gibbs sampling. See also: @@ -44,13 +44,13 @@ (def y (sel ols-data (range 0 2313) 10)) (def lm (linear-model y x :intercept false)) (def param-samp (sample-model-params 5000 lm)) - + ;; view trace plots - (view (trace-plot (:var param-samp ))) + (view (trace-plot (:var param-samp ))) (view (trace-plot (sel (:coefs param-samp) :cols 0))) ;; view histograms - (view (histogram (:var param-samp))) + (view (histogram (:var param-samp))) (view (histogram (sel (:coefs param-samp) :cols 0))) ;; calculate statistics @@ -71,15 +71,15 @@ shape (/ (- (nrow x) (ncol x)) 2) rate (mult 1/2 (mmult (trans resid) resid)) s-sq (div 1 (sample-gamma size :shape shape :rate rate))] - {:coefs - (matrix + {:coefs + (matrix ;(pmap ;; run a parallel map over the values of s-sq (map - (fn [s2] + (fn [s2] (to-list (plus (trans pars) - (mmult (trans (sample-normal (ncol x))) + (mmult (trans (sample-normal (ncol x))) (decomp-cholesky (mult s2 xtxi)))))) - (to-list (trans s-sq)))) + (to-list (trans s-sq)))) :var s-sq}))) @@ -87,7 +87,7 @@ (defn sample-proportions " sample-proportions has been renamed sample-multinomial-params" - ([size counts] + ([size counts] (throw (Exception. "sample-proportions has been renamed sample-multinomial-params")))) @@ -97,7 +97,7 @@ The counts are assumed to have a multinomial distribution. A uniform prior distribution is assigned to the multinomial vector theta, then the posterior distribution of theta is - proportional to a dirichlet distribution with parameters + proportional to a dirichlet distribution with parameters (plus counts 1). @@ -116,13 +116,13 @@ (mean (sel samp-props :cols 2)) (quantile (sel samp-props :cols 2) :probs [0.0275 0.975]) (view (histogram (sel samp-props :cols 2))) - + ;; view a histogram of the difference in proportions between the first ;; two candidates (view (histogram (minus (sel samp-props :cols 0) (sel samp-props :cols 1)))) - - + + " ([size counts] (sample-dirichlet size (plus counts 1)))) @@ -131,7 +131,7 @@ (defn- sample-mvn-params -" Returns samples of means (sampled from an mvn distribution) and vectorized covariance +" Returns samples of means (sampled from an mvn distribution) and vectorized covariance matrices (sampled from an inverse-wishart distribution) for the given mvn data. Arguments: @@ -166,14 +166,14 @@ (use '(incanter core stats bayes charts)) - (def y (sample-mvn 500 :sigma (symmetric-matrix [10 5 10]) :mean [5 2])) + (def y (sample-mvn 500 :sigma (symmetric-matrix [10 5 10]) :mean [5 2])) (def samp (sample-mvn-params 1000 y)) (symmetric-matrix (map mean (trans (:sigmas samp))) :lower false) (map mean (trans (:means samp))) (use '(incanter core stats bayes charts datasets)) - (def cars (to-matrix (get-dataset :cars))) + (def cars (to-matrix (get-dataset :cars))) (def cars-std (trans (map #(sweep (sweep %) :fun div :stat sd) (trans cars)))) (def samp (sample-mvn-params 1000 cars-std)) (symmetric-matrix (map mean (trans (:sigmas samp))) :lower false) @@ -188,19 +188,19 @@ (let [opts (when options (apply assoc {} options)) means (map mean (trans y)) n (count y) - S (reduce plus - (map #(mmult (minus (to-list %) means) - (trans (minus (to-list %) means))) + S (reduce plus + (map #(mmult (minus (to-list %) means) + (trans (minus (to-list %) means))) y)) - sigma-samp (matrix (for [_ (range size)] + sigma-samp (matrix (for [_ (range size)] (half-vectorize (sample-inv-wishart :df (dec n) :scale (solve S))))) mu-samp (matrix (for [sigma sigma-samp] - (sample-mvn 1 - :mean means + (sample-mvn 1 + :mean means :sigma (div (symmetric-matrix sigma :lower false) n)))) ] {:means mu-samp :sigmas sigma-samp}))) - + @@ -250,8 +250,3 @@ v (mmult (trans e) e) s-new (sample-inv-wishart :df df :scale v)] (recur (inc i) (conj coefs b) (conj sigmas (vectorize s-new))))))))) - - - - - src/incanter/bayes.clj @@ -1,27 +1,25 @@ - - -(ns incanter.censored - (:use [incanter.core :only (mult pow)] +(ns incanter.censored + (:use [incanter.core :only (mult pow)] [incanter.stats :only (pdf-normal cdf-normal sd)])) - -(defn- lambda-two-sided - ([a-std b-std] - (/ (- (pdf-normal a-std) (pdf-normal b-std)) + +(defn- lambda-two-sided + ([a-std b-std] + (/ (- (pdf-normal a-std) (pdf-normal b-std)) (- (cdf-normal b-std) (cdf-normal a-std))))) -(defn- lambda-lower +(defn- lambda-lower ([a-std] (/ (pdf-normal a-std) (- 1 (cdf-normal a-std))))) -(defn- lambda-upper - ([b-std] +(defn- lambda-upper + ([b-std] (/ (- (pdf-normal b-std)) (cdf-normal b-std)))) (defn- psi [y mu sigma] (/ (- y mu) sigma)) -(defn censored-mean-two-sided -" Returns the mean of a normal distribution (with mean mu and standard +(defn censored-mean-two-sided +" Returns the mean of a normal distribution (with mean mu and standard deviation sigma) with the lower tail censored at 'a' and the upper tail censored at 'b' " @@ -36,7 +34,7 @@ (defn censored-variance-two-sided -" Returns the variance of a normal distribution (with mean mu and standard +" Returns the variance of a normal distribution (with mean mu and standard deviation sigma) with the lower tail censored at 'a' and the upper tail censored at 'b' " @@ -48,7 +46,7 @@ cdf-a (cdf-normal a-std) cdf-b (cdf-normal b-std)] (+ (* cdf-a (* a a)) - (* (- cdf-b cdf-a) + (* (- cdf-b cdf-a) (+ sigma-sq (* mu mu) (* 2 mu sigma (lambda-two-sided a-std b-std)))) (* sigma-sq (- (* a-std (pdf-normal a-std)) (* b-std (pdf-normal b-std)))) (- (* (- 1 cdf-b) (* b b)) (* Ey Ey)))))) @@ -56,19 +54,19 @@ (defn censored-mean-lower -" Returns the mean of a normal distribution (with mean mu and standard - deviation sigma) with the lower tail censored at 'a' +" Returns the mean of a normal distribution (with mean mu and standard + deviation sigma) with the lower tail censored at 'a' " ([a mu sigma] (let [a-std (psi a mu sigma)] (+ (* (cdf-normal a-std) a) - (* (- 1 (cdf-normal a-std)) + (* (- 1 (cdf-normal a-std)) (+ mu (* sigma (lambda-lower a-std)))))))) (defn censored-variance-lower -" Returns the variance of a normal distribution (with mean mu and standard - deviation sigma) with the lower tail censored at 'a' +" Returns the variance of a normal distribution (with mean mu and standard + deviation sigma) with the lower tail censored at 'a' " ([a mu sigma] (let [a-std (psi a mu sigma) @@ -84,7 +82,7 @@ (defn censored-mean-upper -" Returns the mean of a normal distribution (with mean mu and standard +" Returns the mean of a normal distribution (with mean mu and standard deviation sigma) with the upper tail censored at 'b' " ([b mu sigma] @@ -94,7 +92,7 @@ (defn censored-variance-upper -" Returns the variance of a normal distribution (with mean mu and standard +" Returns the variance of a normal distribution (with mean mu and standard deviation sigma) with the upper tail censored at 'b' " ([b mu sigma] @@ -110,7 +108,7 @@ -(defn truncated-variance +(defn truncated-variance " Returns the variance of a normal distribution truncated at a and b. Options: @@ -120,9 +118,9 @@ :b (default +infinity) upper truncation point Examples: - + (use '(incanter core stats)) - (truncated-variance :a -1.96 :b 1.96) + (truncated-variance :a -1.96 :b 1.96) (truncated-variance :a 0) (truncated-variance :b 0) @@ -138,9 +136,9 @@ References: - DeMaris, A. (2004) Regression with social data: modeling continuous and limited response variables. + DeMaris, A. (2004) Regression with social data: modeling continuous and limited response variables. Wiley-IEEE. - + http://en.wikipedia.org/wiki/Truncated_normal_distribution " ([& options] @@ -152,11 +150,11 @@ b (or (:b opts) Double/POSITIVE_INFINITY) lambda (fn [alpha] (/ (pdf-normal alpha) (- 1 (cdf-normal alpha)))) delta (fn [alpha] (* (lambda alpha) (- (lambda alpha) alpha))) - ;one-tail-var (fn [alpha s-sq] - ; (* s-sq - ; (- 1 (cdf-normal alpha)) - ; (+ (- 1 (delta alpha)) - ; (* (pow (- alpha (lambda alpha)) 2) + ;one-tail-var (fn [alpha s-sq] + ; (* s-sq + ; (- 1 (cdf-normal alpha)) + ; (+ (- 1 (delta alpha)) + ; (* (pow (- alpha (lambda alpha)) 2) ; (cdf-normal alpha))))) a-std (if (= a Double/NEGATIVE_INFINITY) Double/NEGATIVE_INFINITY (/ (- a mu) sd)) b-std (if (= b Double/POSITIVE_INFINITY) Double/POSITIVE_INFINITY (/ (- b mu) sd)) @@ -164,7 +162,7 @@ pdf-b (if (= b Double/POSITIVE_INFINITY) 0 (pdf-normal b-std)) cdf-a (if (= a Double/NEGATIVE_INFINITY) 0 (cdf-normal a-std)) cdf-b (if (= b Double/POSITIVE_INFINITY) 1 (cdf-normal b-std))] - (cond + (cond (and (= b Double/POSITIVE_INFINITY) (= a Double/NEGATIVE_INFINITY)) sigma-sq (and (= b Double/POSITIVE_INFINITY) (> a Double/NEGATIVE_INFINITY)) @@ -174,9 +172,7 @@ (* sigma-sq (- 1 (delta (- 1 b-std)))) ;(- sigma-sq (one-tail-var b-std sigma-sq)) :else - (* sigma-sq + (* sigma-sq (+ 1 (/ (- (* a-std pdf-a) (* b-std pdf-b)) (- cdf-b cdf-a)) (- (pow (/ (- pdf-a pdf-b) (- cdf-b cdf-a)) 2)))))))) - - src/incanter/censored.clj @@ -16,9 +16,9 @@ -(ns incanter.charts +(ns incanter.charts ;(:use (incanter core stats)) - (:use [incanter.core :only (matrix? to-list plus minus div group-by + (:use [incanter.core :only (matrix? to-list plus minus div group-by bind-columns view save)] [incanter.stats :only (quantile quantile-normal cumulative-mean sd)]) (:import (java.io File) @@ -36,14 +36,14 @@ (org.jfree.data.category DefaultCategoryDataset) (org.jfree.chart.renderer.xy XYLineAndShapeRenderer) (org.jfree.ui TextAnchor) - (org.jfree.chart.annotations XYPointerAnnotation - XYTextAnnotation + (org.jfree.chart.annotations XYPointerAnnotation + XYTextAnnotation XYPolygonAnnotation) )) -(defmacro histogram +(defmacro histogram " Returns a JFreeChart object representing the histogram of the given data. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -75,7 +75,7 @@ # plot some gamma data (def gam-hist (histogram (sample-gamma 1000) :density true :nbins 30)) - (view gam-hist) + (view gam-hist) (def x (range 0 8 0.01)) (add-lines gam-hist x (pdf-gamma x)) @@ -84,7 +84,7 @@ References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([x & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -93,7 +93,7 @@ density?# (true? (:density opts#)) main-title# (or (:title opts#) "Histogram") x-lab# (or (:x-label opts#) (str '~x)) - y-lab# (or (:y-label opts#) + y-lab# (or (:y-label opts#) (if density?# "Density" "Frequency")) series-lab# (or (:series-label opts#) (str '~x)) legend?# (true? (:legend opts#)) @@ -102,13 +102,13 @@ (do (.addSeries dataset# series-lab# (double-array data#) nbins#) (when density?# (.setType dataset# org.jfree.data.statistics.HistogramType/SCALE_AREA_TO_1)) - (org.jfree.chart.ChartFactory/createHistogram + (org.jfree.chart.ChartFactory/createHistogram main-title# x-lab# y-lab# - dataset# - org.jfree.chart.plot.PlotOrientation/VERTICAL - legend?# ; no legend + dataset# + org.jfree.chart.plot.PlotOrientation/VERTICAL + legend?# ; no legend true ; tooltips false))))) @@ -116,7 +116,7 @@ (declare add-points) -(defmacro scatter-plot +(defmacro scatter-plot " Returns a JFreeChart object representing a scatter-plot of the given data. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -155,13 +155,13 @@ ;; plot the first two columns grouped by the fifth column (view (scatter-plot (sel iris :cols 0) (sel iris :cols 1) :group-by (sel iris :cols 4))) - + ;; see INCANTER_HOME/examples/probability_plots.clj for more examples of plots References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([x y & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -182,31 +182,31 @@ chart# (do (doseq [i# (range (count x#))] (.add data-series# (nth x# i#) (nth y# i#))) (.addSeries dataset# data-series#) - (org.jfree.chart.ChartFactory/createScatterPlot + (org.jfree.chart.ChartFactory/createScatterPlot main-title# x-lab# y-lab# - dataset# - org.jfree.chart.plot.PlotOrientation/VERTICAL - legend?# + dataset# + org.jfree.chart.plot.PlotOrientation/VERTICAL + legend?# true ; tooltips false)) _# (when x-groups# - (doseq [i# (range 1 (count x-groups#))] - (add-points chart# (nth x-groups# i#) + (doseq [i# (range 1 (count x-groups#))] + (add-points chart# (nth x-groups# i#) (nth y-groups# i#) :series-label (format "%s, %s (%s)" '~x '~y i#))))] chart#))) -(defn line-plot +(defn line-plot "WARNING: line-plot has been renamed xy-plot." ([x y & options] (throw (Exception. "line-plot has been renamed xy-plot")))) (declare add-lines) -(defmacro xy-plot +(defmacro xy-plot " Returns a JFreeChart object representing a xy-plot of the given data. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -227,13 +227,13 @@ (use '(incanter core stats charts)) ;; plot the cosine function - (def x (range -1 5 0.01)) + (def x (range -1 5 0.01)) (def y (cos (mult 2 Math/PI x))) (view (xy-plot x y)) ;; plot gamma pdf with different parameters (def x2 (range 0 20 0.1)) - (def gamma-plot (xy-plot x2 (pdf-gamma x2 :shape 1 :rate 2) + (def gamma-plot (xy-plot x2 (pdf-gamma x2 :shape 1 :rate 2) :legend true :title \"Gamma PDF\" :y-label \"Density\")) @@ -245,17 +245,17 @@ ;; use :group-by option (use '(incanter core charts datasets)) - (def data (to-matrix (get-dataset :chick-weight))) + (def data (to-matrix (get-dataset :chick-weight))) (let [[weight time chick] (trans data)] (view (xy-plot time weight :group-by chick))) - + ;; see INCANTER_HOME/examples/probability_plots.clj for more examples of plots - + References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([x y & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -276,18 +276,18 @@ chart# (do (doseq [i# (range (count x#))] (.add data-series# (nth x# i#) (nth y# i#))) (.addSeries dataset# data-series#) - (org.jfree.chart.ChartFactory/createXYLineChart + (org.jfree.chart.ChartFactory/createXYLineChart main-title# x-lab# y-lab# - dataset# - org.jfree.chart.plot.PlotOrientation/VERTICAL - legend?# + dataset# + org.jfree.chart.plot.PlotOrientation/VERTICAL + legend?# true ; tooltips false)) _# (when x-groups# - (doseq [i# (range 1 (count x-groups#))] - (add-lines chart# (nth x-groups# i#) + (doseq [i# (range 1 (count x-groups#))] + (add-lines chart# (nth x-groups# i#) (nth y-groups# i#) :series-label (format "%s, %s (%s)" '~x '~y i#))))] chart#))) @@ -296,8 +296,8 @@ (defmacro function-plot -" Returns a xy-plot object of the given function over the range indicated - by the min-range and max-range arguments. Use the 'view' function to +" Returns a xy-plot object of the given function over the range indicated + by the min-range and max-range arguments. Use the 'view' function to display the chart, or the 'save' function to write it to a file. Options: @@ -321,7 +321,7 @@ (defn cubic [x] (+ (* x x x) (* 2 x x) (* 2 x) 3)) (view (function-plot cubic -10 10)) - + " ([function min-range max-range & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -332,11 +332,11 @@ y-lab# (or (:y-label opts#) (str '~function)) series-lab# (or (:series-label opts#) (format "%s" '~function)) legend?# (true? (:legend opts#))] - (xy-plot x# (map ~function x#) - :x-label x-lab# - :y-label y-lab# - :title main-title# - :series-label series-lab# + (xy-plot x# (map ~function x#) + :x-label x-lab# + :y-label y-lab# + :title main-title# + :series-label series-lab# :legend legend?#)))) @@ -345,7 +345,7 @@ (declare add-box-plot) -(defmacro box-plot +(defmacro box-plot " Returns a JFreeChart object representing a box-plot of the given data. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -364,24 +364,24 @@ Examples: (use '(incanter core stats charts)) - (def gamma-box-plot (box-plot (sample-gamma 1000 :shape 1 :rate 2) + (def gamma-box-plot (box-plot (sample-gamma 1000 :shape 1 :rate 2) :title \"Gamma Boxplot\" - :legend true)) + :legend true)) (view gamma-box-plot) (add-box-plot gamma-box-plot (sample-gamma 1000 :shape 2 :rate 2)) (add-box-plot gamma-box-plot (sample-gamma 1000 :shape 3 :rate 2)) ;; use the group-by options (use '(incanter core stats datasets charts)) - (def iris (to-matrix (get-dataset :iris))) + (def iris (to-matrix (get-dataset :iris))) (view (box-plot (sel iris :cols 0) :group-by (sel iris :cols 4) :legend true)) ;; see INCANTER_HOME/examples/probability_plots.clj for more examples of plots - + References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([x & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -396,11 +396,11 @@ category-label# (or (:category-label opts#) 0) legend?# (true? (:legend opts#)) dataset# (DefaultBoxAndWhiskerCategoryDataset.) - chart# (org.jfree.chart.ChartFactory/createBoxAndWhiskerChart + chart# (org.jfree.chart.ChartFactory/createBoxAndWhiskerChart main-title# x-label# y-label# - dataset# + dataset# legend?#)] (do (.. chart# getCategoryPlot getRenderer (setMaximumBarWidth 0.25)) @@ -412,7 +412,7 @@ -(defmacro bar-chart +(defmacro bar-chart " Returns a JFreeChart object representing a bar-chart of the given data. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -420,14 +420,14 @@ Arguments: categories -- a sequence of categories values -- a sequence of numeric values - + Options: :title (default 'Histogram') main title :x-label (default 'Categories') :y-label (default 'Value') :legend (default false) prints legend :vertical (default true) the orientation of the plot - :group-by (default nil) -- a vector of values used to group the values into + :group-by (default nil) -- a vector of values used to group the values into series within each category. @@ -443,10 +443,10 @@ (def grass-type (sel data :cols 1)) (def treatment-type (sel data :cols 2)) (def uptake (sel data :cols 4)) - (view (bar-chart grass-type uptake + (view (bar-chart grass-type uptake :title \"CO2 Uptake\" - :group-by treatment-type - :x-label \"Grass Types\" :y-label \"Uptake\" + :group-by treatment-type + :x-label \"Grass Types\" :y-label \"Uptake\" :legend true)) @@ -460,7 +460,7 @@ (def data (get-dataset :austres)) (view data) - (def plot (bar-chart (sel data :cols 0) (sel data :cols 1) + (def plot (bar-chart (sel data :cols 0) (sel data :cols 1) :group-by (sel data :cols 2) :legend true)) (view plot) (save plot \"/tmp/austres_plot.png\" :width 1000) @@ -474,18 +474,18 @@ (view (bar-chart [\"a\" \"b\" \"c\"] [10 20 30])) (view (bar-chart [\"a\" \"a\" \"b\" \"b\" \"c\" \"c\" ] [10 20 30 10 40 20] - :legend true + :legend true :group-by [\"I\" \"II\" \"I\" \"II\" \"I\" \"II\"])) - (view (bar-chart (sample \"abcdefghij\" :size 10 :replacement true) + (view (bar-chart (sample \"abcdefghij\" :size 10 :replacement true) (sample-uniform 10 :max 50) :legend true)) - + References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([categories values & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -501,16 +501,16 @@ main-title# x-label# y-label# - dataset# - (if vertical?# + dataset# + (if vertical?# org.jfree.chart.plot.PlotOrientation/VERTICAL org.jfree.chart.plot.PlotOrientation/HORIZONTAL) legend?# true false)] (do - (doseq [i# (range 0 (count values#))] (.addValue dataset# - (nth values# i#) + (doseq [i# (range 0 (count values#))] (.addValue dataset# + (nth values# i#) (if group-by# (nth group-by# i#) (str '~values)) @@ -520,7 +520,7 @@ -(defmacro line-chart +(defmacro line-chart " Returns a JFreeChart object representing a line-chart of the given values and categories. Use the 'view' function to display the chart, or the 'save' function to write it to a file. @@ -528,13 +528,13 @@ Arguments: categories -- a sequence of categories values -- a sequence of numeric values - + Options: :title (default 'Histogram') main title :x-label (default 'Categories') :y-label (default 'Value') :legend (default false) prints legend - :group-by (default nil) -- a vector of values used to group the values into + :group-by (default nil) -- a vector of values used to group the values into series within each category. @@ -544,7 +544,7 @@ Examples: (use '(incanter core stats charts datasets)) - + (def data (get-dataset :airline-passengers)) (def years (sel data :cols 0)) (def months (sel data :cols 2)) @@ -560,15 +560,15 @@ (view (line-chart [\"a\" \"b\" \"c\" \"d\" \"e\" \"f\"] [10 20 30 10 40 20])) - (view (line-chart (sample \"abcdefghij\" :size 10 :replacement true) + (view (line-chart (sample \"abcdefghij\" :size 10 :replacement true) (sample-uniform 10 :max 50) :legend true)) - + References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([categories values & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -584,16 +584,16 @@ main-title# x-label# y-label# - dataset# - (if vertical?# + dataset# + (if vertical?# org.jfree.chart.plot.PlotOrientation/VERTICAL org.jfree.chart.plot.PlotOrientation/HORIZONTAL) legend?# true false)] (do - (doseq [i# (range 0 (count values#))] (.addValue dataset# - (nth values# i#) + (doseq [i# (range 0 (count values#))] (.addValue dataset# + (nth values# i#) (if group-by# (nth group-by# i#) (str '~values)) @@ -602,7 +602,7 @@ -(defmacro add-histogram +(defmacro add-histogram " Adds a histogram to an existing histogram plot, returns the modified chart object. @@ -614,15 +614,15 @@ Examples: (use '(incanter core charts stats)) - (doto (histogram (sample-normal 1000) + (doto (histogram (sample-normal 1000) :legend true) - view + view (add-histogram (sample-normal 1000 :sd 0.5))) References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([chart x & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -640,7 +640,7 @@ chart#)))) -(defmacro add-box-plot +(defmacro add-box-plot " Adds an additional box to an existing box-plot, returns the modified chart object. @@ -658,17 +658,17 @@ References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([chart x & options] `(let [opts# (when '~options (assoc {} ~@options)) _x# (if (matrix? ~x) (to-list ~x) ~x) chart# ~chart data-plot# (.getCategoryPlot chart#) - n-col# (.getColumnCount (.getDataset data-plot#)) + n-col# (.getColumnCount (.getDataset data-plot#)) n-row# (.getRowCount (.getDataset data-plot#)) series-label# (or (:series-label opts#) (str '~x)) - category-label# (or (:category-label opts#) + category-label# (or (:category-label opts#) (str n-col#)) ] (do @@ -702,7 +702,7 @@ References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([chart categories values & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -710,11 +710,11 @@ chart# ~chart group-by# (:group-by opts#) data-plot# (.getCategoryPlot chart#) - n-col# (.getColumnCount (.getDataset data-plot#)) + n-col# (.getColumnCount (.getDataset data-plot#)) n-row# (.getRowCount (.getDataset data-plot#))] (do (doseq [i# (range 0 (count values#))] (.addValue (.getDataset data-plot#) - (nth values# i#) + (nth values# i#) (if group-by# (nth group-by# i#) (str '~values)) @@ -726,90 +726,90 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn set-alpha -" Sets the alpha level (transparancy) of the plot's foreground, +(defn set-alpha +" Sets the alpha level (transparancy) of the plot's foreground, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart alpha] + ([chart alpha] (.setForegroundAlpha (.getPlot chart) alpha) chart)) -(defn set-background-alpha -" Sets the alpha level (transparancy) of the plot's background, +(defn set-background-alpha +" Sets the alpha level (transparancy) of the plot's background, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart alpha] + ([chart alpha] (.setBackgroundAlpha (.getPlot chart) alpha) chart)) -(defn clear-background -" Sets the alpha level (transparancy) of the plot's background to zero, +(defn clear-background +" Sets the alpha level (transparancy) of the plot's background to zero, removing the default grid, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart] + ([chart] (.setBackgroundAlpha (.getPlot chart) 0.0) chart)) -(defn set-title +(defn set-title " Sets the main title of the plot, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart title] + ([chart title] (.setTitle chart title) chart)) -(defn set-x-label +(defn set-x-label " Sets the label of the x-axis, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart label] + ([chart label] (.setLabel (.getDomainAxis (.getPlot chart)) label) chart)) -(defn set-y-label +(defn set-y-label " Sets the label of the y-axis, returns the modified chart object. References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " - ([chart label] + ([chart label] (.setLabel (.getRangeAxis (.getPlot chart)) label) chart)) - + (defmacro add-lines " Plots lines on the given scatter or line plot of the (x,y) points. Equivalent to R's lines function, returns the modified chart object. @@ -825,7 +825,7 @@ (def x (sel cars :cols 1)) (def plot1 (scatter-plot x y :legend true)) (view plot1) - + ;; add regression line to scatter plot (def lm1 (linear-model y x)) (add-lines plot1 x (:fitted lm1)) @@ -845,7 +845,7 @@ References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([chart x y & options] @@ -866,7 +866,7 @@ (doseq [i# (range (count _x#))] (.add data-series# (nth _x# i#) (nth _y# i#))) ;(.addSeries (.getDataset data-plot) data-series) (.setSeriesRenderingOrder (.getPlot chart#) org.jfree.chart.plot.SeriesRenderingOrder/FORWARD) - (.setDatasetRenderingOrder data-plot# org.jfree.chart.plot.DatasetRenderingOrder/FORWARD) + (.setDatasetRenderingOrder data-plot# org.jfree.chart.plot.DatasetRenderingOrder/FORWARD) (.addSeries data-set# data-series#) (.setDataset data-plot# (.getDatasetCount data-plot#) data-set#) ;(.setSeriesPaint line-renderer 1 java.awt.Color/blue) @@ -916,7 +916,7 @@ (doto (function-plot cubic -10 10) (add-function deriv-cubic -10 10) view) - + " ([chart function min-range max-range & options] `(let [opts# (when '~options (assoc {} ~@options)) @@ -924,7 +924,7 @@ x# (range ~min-range ~max-range step-size#) series-lab# (or (:series-label opts#) (format "%s" '~function)) ] - (add-lines ~chart x# (map ~function x#) + (add-lines ~chart x# (map ~function x#) :series-label series-lab#)))) @@ -942,7 +942,7 @@ (def cars (to-matrix (get-dataset :cars))) (def y (sel cars :cols 0)) (def x (sel cars :cols 1)) - + ;; add regression line to scatter plot (def lm1 (linear-model y x)) ;; model the data without an intercept @@ -957,7 +957,7 @@ References: http://www.jfree.org/jfreechart/api/javadoc/ http://www.jfree.org/jfreechart/api/javadoc/org/jfree/chart/JFreeChart.html - + " ([chart x y & options] @@ -1002,7 +1002,7 @@ Options: :text -- (default \"\") text to include at the end of the arrow - :angle -- (default :nw) either a number indicating the angle of the arrow, + :angle -- (default :nw) either a number indicating the angle of the arrow, or a keyword indicating a direction (:north :nw :west :sw :south :se :east :ne) @@ -1045,7 +1045,7 @@ angle (if (nil? (:angle opts)) nw (if (keyword? (:angle opts)) - (cond + (cond (= (:angle opts) :north) n (= (:angle opts) :south) s (= (:angle opts) :west) w @@ -1104,17 +1104,17 @@ ;; project the first four dimension of the iris data onto the first ;; two principal components - (def x1 (mmult (sel iris :cols (range 4)) pc1)) - (def x2 (mmult (sel iris :cols (range 4)) pc2)) - + (def x1 (mmult (sel iris :cols (range 4)) pc1)) + (def x2 (mmult (sel iris :cols (range 4)) pc2)) + ;; now plot the transformed data, coloring each species a different color - (def plot (scatter-plot x1 x2 + (def plot (scatter-plot x1 x2 :group-by (sel iris :cols 4) :x-label \"PC1\" :y-label \"PC2\" :title \"Iris PCA\")) (view plot) ;; add some text annotations (add-text plot -2.5 -6.5 \"Setosa\") - (add-text plot -5 -5.5 \"Versicolor\") + (add-text plot -5 -5.5 \"Versicolor\") (add-text plot -8 -5.5 \"Virginica\") @@ -1127,7 +1127,7 @@ (defn add-polygon -" Adds a polygon outline defined by a given coordinates. The last coordinate will +" Adds a polygon outline defined by a given coordinates. The last coordinate will close with the first. If only two points are given, it will plot a line. Arguments: @@ -1140,7 +1140,7 @@ (def x (range -3 3 0.01)) (def plot (xy-plot x (pdf-normal x))) (view plot) - + ;; add polygon to the chart (add-polygon plot [[-1.96 0] [1.96 0] [1.96 0.4] [-1.96 0.4]]) ;; the coordinates can also be passed in a matrix @@ -1161,11 +1161,11 @@ ;; project the first four dimension of the iris data onto the first ;; two principal components - (def x1 (mmult (sel iris :cols (range 4)) pc1)) - (def x2 (mmult (sel iris :cols (range 4)) pc2)) - + (def x1 (mmult (sel iris :cols (range 4)) pc1)) + (def x2 (mmult (sel iris :cols (range 4)) pc2)) + ;; now plot the transformed data, coloring each species a different color - (def plot (scatter-plot x1 x2 + (def plot (scatter-plot x1 x2 :group-by (sel iris :cols 4) :x-label \"PC1\" :y-label \"PC2\" :title \"Iris PCA\")) @@ -1174,7 +1174,7 @@ (add-polygon plot [[-3.2 -6.3] [-2 -6.3] [-2 -3.78] [-3.2 -3.78]]) ;; add some text annotations (add-text plot -2.5 -6.5 \"Setosa\") - (add-text plot -5 -5.5 \"Versicolor\") + (add-text plot -5 -5.5 \"Versicolor\") (add-text plot -8 -5.5 \"Virginica\") @@ -1203,7 +1203,7 @@ (def x (sel ols-data (range 0 2313) (range 1 10))) (def y (sel ols-data (range 0 2313) 10)) (def sample-params (sample-model-params 5000 (linear-model y x :intercept false))) - (view (trace-plot (:var sample-params))) + (view (trace-plot (:var sample-params))) (view (trace-plot (sel (:coefs sample-params) :cols 0))) @@ -1216,10 +1216,10 @@ series-lab (or (:series-label opts) "Value") ;legend? (or (:series-label opts) true) n (count x) - chart (xy-plot (range n) + chart (xy-plot (range n) x ;)] - :title title - :x-label x-label + :title title + :x-label x-label :y-label y-label :series-label series-lab)] (do @@ -1239,7 +1239,7 @@ http://en.wikipedia.org/wiki/QQ_plot Examples: - + (use '(incanter core stats charts)) (view (qq-plot (sample-normal 100))) (view (qq-plot (sample-exp 100))) @@ -1252,7 +1252,7 @@ quants (for [k (range 1 n)] (/ k (inc n))) norm-quants (quantile-normal quants) y (quantile x :probs quants)] - (scatter-plot norm-quants y + (scatter-plot norm-quants y :title "QQ-Plot" :x-label "Normal theoretical quantiles" :y-label "Data quantiles" @@ -1279,7 +1279,7 @@ " ([x1 x2] - (let [plot (scatter-plot (div (plus x1 x2) 2) (minus x1 x2) + (let [plot (scatter-plot (div (plus x1 x2) 2) (minus x1 x2) :title "Bland Altman Plot" :legend false) x-axis (div (plus x1 x2) 2) @@ -1300,7 +1300,7 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defmethod view org.jfree.chart.JFreeChart +(defmethod view org.jfree.chart.JFreeChart ([chart & options] (let [opts (when options (apply assoc {} options)) window-title (or (:window-title opts) "Incanter Plot") @@ -1320,6 +1320,3 @@ height (or (:height opts) 400)] (ChartUtilities/saveChartAsPNG (File. filename) chart width height) nil))) - - - src/incanter/charts.clj @@ -16,18 +16,18 @@ -(ns incanter.core +(ns incanter.core (:use (incanter internal)) (:import (incanter Matrix) - (cern.colt.matrix.tdouble DoubleMatrix2D - DoubleFactory2D + (cern.colt.matrix.tdouble DoubleMatrix2D + DoubleFactory2D DoubleFactory1D) (cern.colt.matrix.tdouble.algo DoubleAlgebra DoubleFormatter) (cern.colt.matrix.tdouble.algo.decomposition DoubleCholeskyDecomposition DoubleSingularValueDecompositionDC DoubleEigenvalueDecomposition - DoubleLUDecomposition + DoubleLUDecomposition DoubleQRDecomposition) (cern.jet.math.tdouble DoubleFunctions DoubleArithmetic) (cern.colt.function.tdouble DoubleDoubleFunction DoubleFunction) @@ -37,7 +37,7 @@ (java.util Vector))) -(defn matrix +(defn matrix " Returns an instance of an incanter.Matrix, which is an extension of cern.colt.matrix.tdouble.impl.DenseColDoubleMatrix2D that implements the Clojure @@ -57,17 +57,17 @@ (rest (first A)) ; produces a row matrix [2 3] ; since (plus row1 row2) adds the two rows element-by-element, - (reduce plus A) ; produces the sums of the columns, + (reduce plus A) ; produces the sums of the columns, - ; and since (sum row1) sums the elements of the row, - (map sum A) ; produces the sums of the rows, + ; and since (sum row1) sums the elements of the row, + (map sum A) ; produces the sums of the rows, ; you can filter the rows using Clojure's filter function (filter #(> (nth % 1) 4) A) ; returns the rows where the second column is greater than 4. References: http://incanter.org/docs/parallelcolt/api/cern/colt/matrix/tdouble/DoubleMatrix2D.html - + " ([data] (make-matrix data)) @@ -79,23 +79,23 @@ (make-matrix init-val rows cols))) -(defn matrix? +(defn matrix? " Test if obj is 'derived' incanter.Matrix." ([obj] (is-matrix obj))) -(defn nrow +(defn nrow " Returns the number of rows in the given matrix. Equivalent to R's nrow function." ([mat] - (cond + (cond (matrix? mat) (.rows #^Matrix mat) (coll? mat) (count mat)))) -(defn ncol +(defn ncol " Returns the number of columns in the given matrix. Equivalent to R's ncol function." ([mat] - (cond + (cond (matrix? mat) (.columns #^Matrix mat) (coll? mat) 1 ))) @@ -108,9 +108,9 @@ -(defn identity-matrix +(defn identity-matrix " Returns an n-by-n identity matrix. - + Examples: (identity-matrix 4) @@ -118,11 +118,11 @@ ([n] (Matrix. (.identity DoubleFactory2D/dense n)))) -(defn diag +(defn diag " If given a matrix, diag returns a sequence of its diagonal elements. - If given a sequence, it returns a matrix with the sequence's elements + If given a sequence, it returns a matrix with the sequence's elements on its diagonal. Equivalent to R's diag function. - + Examples: (diag [1 2 3 4]) @@ -134,7 +134,7 @@ " ([m] - (cond + (cond (matrix? m) (seq (.toArray (.diagonal DoubleFactory2D/dense m))) (coll? m) @@ -143,19 +143,19 @@ m))) -(defn #^Matrix trans +(defn #^Matrix trans " Returns the transpose of the given matrix. Equivalent to R's t function Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) (trans A) - + " ([mat] - (cond + (cond (matrix? mat) (.viewDice #^Matrix mat) (coll? mat) @@ -163,7 +163,7 @@ -(defn- except-for +(defn- except-for " Returns a lazy list of numbers ranging from 0 to n, except for the given exceptions. Examples: @@ -172,27 +172,27 @@ (except-for 10 [5 7]) " - ([n exceptions] + ([n exceptions] (let [except (if (coll? exceptions) exceptions [exceptions])] (for [i (range n) :when (reduce #(and %1 %2) (map #(not= i %) except))] i)))) -(defmulti sel +(defmulti sel " - Returns an element or subset of the given matrix, or dataset. + Returns an element or subset of the given matrix, or dataset. Argument: a matrix object or dataset. Options: - :rows (default true) + :rows (default true) returns all rows by default, can pass a row index or sequence of row indices - :cols (default true) + :cols (default true) returns all columns by default, can pass a column index or sequence of column indices :except-rows (default nil) can pass a row index or sequence of row indices to exclude :except-cols (default nil) can pass a column index or sequence of column indices to exclude - :filter (default nil) + :filter (default nil) a function can be provided to filter the rows of the matrix Examples: @@ -249,16 +249,16 @@ (let [opts (when options (apply assoc {} options)) except-rows (:except-rows opts) except-columns (:except-cols opts) - rows (cond - (:rows opts) - (:rows opts) + rows (cond + (:rows opts) + (:rows opts) except-rows (except-for (.rows mat) except-rows) :else true) - cols (cond - (:cols opts) - (:cols opts) + cols (cond + (:cols opts) + (:cols opts) except-columns (except-for (.columns mat) except-columns) :else @@ -285,15 +285,15 @@ (and (true? rows) (true? cols)) mat)))) - -(defn bind-rows -" Returns the matrix resulting from concatenating the given matrices + +(defn bind-rows +" Returns the matrix resulting from concatenating the given matrices and/or sequences by their rows. Equivalent to R's rbind. Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) (def B (matrix [[10 11 12] @@ -306,8 +306,8 @@ " ([& args] (reduce - (fn [A B] - (cond + (fn [A B] + (cond (nil? (seq A)) B (nil? (seq B)) @@ -321,17 +321,17 @@ (and (coll? A) (coll? B)) (conj (matrix A (count A)) (matrix B (count B))) :else - (throw (Exception. "Incompatible types")))) + (throw (Exception. "Incompatible types")))) args))) -(defn bind-columns -" Returns the matrix resulting from concatenating the given matrices +(defn bind-columns +" Returns the matrix resulting from concatenating the given matrices and/or sequences by their columns. Equivalent to R's cbind. Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) (def B (matrix [10 11 12])) @@ -343,12 +343,12 @@ " ([& args] - (reduce + (reduce (fn [A B] (.viewDice (bind-rows (trans A) (trans B)))) args))) -;(defn inner-product [& args] (apply + (apply map * args))) +;(defn inner-product [& args] (apply + (apply map * args))) ;(inner-product [1 2 3] [4 5 6]) ; = 32 @@ -359,16 +359,16 @@ -(defn plus -" Performs element-by-element addition on multiple matrices, sequences, +(defn plus +" Performs element-by-element addition on multiple matrices, sequences, and/or numbers. Equivalent to R's + operator. - + Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) - (plus A A A) + (plus A A A) (plus A 2) (plus 2 A) (plus [1 2 3] [1 2 3]) @@ -380,19 +380,19 @@ ([& args] (reduce (fn [A B] (combine-with A B clojure.core/+ plus)) args))) -(defn minus -" Performs element-by-element subtraction on multiple matrices, sequences, - and/or numbers. If only a single argument is provided, returns the +(defn minus +" Performs element-by-element subtraction on multiple matrices, sequences, + and/or numbers. If only a single argument is provided, returns the negative of the given matrix, sequence, or number. Equivalent to R's - operator. Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) (minus A) - (minus A A A) + (minus A A A) (minus A 2) (minus 2 A) (minus [1 2 3] [1 2 3]) @@ -408,37 +408,37 @@ -(defn mult -" Performs element-by-element multiplication on multiple matrices, sequences, +(defn mult +" Performs element-by-element multiplication on multiple matrices, sequences, and/or numbers. Equivalent to R's * operator. - + Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) - (mult A A A) + (mult A A A) (mult A 2) (mult 2 A) (mult [1 2 3] [1 2 3]) (mult [1 2 3] 2) (mult 2 [1 2 3]) - + " ([& args] (reduce (fn [A B] (combine-with A B clojure.core/* mult)) args))) -(defn div -" Performs element-by-element division on multiple matrices, sequences, +(defn div +" Performs element-by-element division on multiple matrices, sequences, and/or numbers. Equivalent to R's / operator. - + Examples: - (def A (matrix [[1 2 3] - [4 5 6] + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) - (div A A A) + (div A A A) (div A 2) (div 2 A) (div [1 2 3] [1 2 3]) @@ -453,98 +453,98 @@ (reduce (fn [A B] (combine-with A B clojure.core// div)) args)))) -(defn pow +(defn pow " This is an element-by-element exponent function, raising the first argument, by the exponents in the remaining arguments. Equivalent to R's ^ operator." ([& args] (reduce (fn [A B] (combine-with A B #(Math/pow %1 %2) pow)) args))) -(defn atan2 +(defn atan2 "Returns the atan2 of the elements in the given matrices, sequences or numbers. Equivalent to R's atan2 function." - ([& args] (reduce (fn [A B] (combine-with A B #(Math/atan2 %1 %2) + ([& args] (reduce (fn [A B] (combine-with A B #(Math/atan2 %1 %2) cern.jet.math.tdouble.DoubleFunctions/atan2)) args))) -(defn sqrt +(defn sqrt "Returns the square-root of the elements in the given matrix, sequence or number. Equivalent to R's sqrt function." ([A] (pow A 1/2))) -(defn sq +(defn sq "Returns the square of the elements in the given matrix, sequence or number. Equivalent to R's sq function." ([A] (mult A A))) -(defn log +(defn log "Returns the natural log of the elements in the given matrix, sequence or number. Equvalent to R's log function." ([A] (transform-with A #(Math/log %) log))) -(defn log2 +(defn log2 "Returns the log base 2 of the elements in the given matrix, sequence or number. Equivalent to R's log2 function." ([A] (transform-with A #(/ (Math/log %) (Math/log 2)) log2))) -(defn log10 +(defn log10 "Returns the log base 10 of the elements in the given matrix, sequence or number. Equivalent to R's log10 function." ([A] (transform-with A #(Math/log10 %) (lg 10.0)))) -(defn exp +(defn exp "Returns the exponential of the elements in the given matrix, sequence or number. Equivalent to R's exp function." ([A] (transform-with A #(Math/exp %) exp))) -(defn abs +(defn abs "Returns the absolute value of the elements in the given matrix, sequence or number. Equivalent to R's abs function." ([A] (transform-with A #(Math/abs (float %)) abs))) -(defn sin +(defn sin "Returns the sine of the elements in the given matrix, sequence or number. Equivalent to R's sin function." ([A] (transform-with A #(Math/sin %) sin))) -(defn asin +(defn asin "Returns the arc sine of the elements in the given matrix, sequence or number. Equivalent to R's asin function." ([A] (transform-with A #(Math/asin %) asin))) -(defn cos +(defn cos "Returns the cosine of the elements in the given matrix, sequence or number. Equivalent to R's cos function." ([A] (transform-with A #(Math/cos %) cos))) -(defn acos +(defn acos "Returns the arc cosine of the elements in the given matrix, sequence or number. Equivalent to R's acos function." ([A] (transform-with A #(Math/acos %) acos))) -(defn tan +(defn tan "Returns the tangent of the elements in the given matrix, sequence or number. Equivalent to R's tan function." ([A] (transform-with A #(Math/tan %) tan))) -(defn atan +(defn atan "Returns the arc tangent of the elements in the given matrix, sequence or number. Equivalent to R's atan function." ([A] (transform-with A #(Math/atan %) atan))) -(defn factorial +(defn factorial " Returns the factorial of k (k must be a positive integer). Equivalent to R's factorial function. @@ -561,10 +561,10 @@ -(defn choose +(defn choose " - Returns number of k-combinations (each of size k) from a set S with - n elements (size n), which is the binomial coefficient (also known + Returns number of k-combinations (each of size k) from a set S with + n elements (size n), which is the binomial coefficient (also known as the 'choose function') [wikipedia] choose = n!/(k!(n - k)!) @@ -586,13 +586,13 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn to-list +(defn to-list " Returns a list-of-lists if the given matrix is two-dimensional, and a flat list if the matrix is one-dimensional." - ([#^Matrix mat] + ([#^Matrix mat] (cond (and (coll? mat) (not (matrix? mat))) - mat + mat (= (.columns mat) 1) (first (map #(seq %) (seq (.toArray (.viewDice mat))))) (= (.rows mat) 1) @@ -601,19 +601,19 @@ (map #(seq %) (seq (.toArray mat)))))) -(defn #^Matrix copy +(defn #^Matrix copy "Returns a copy of the given matrix." ([#^Matrix mat] (.copy mat))) -(defn mmult +(defn mmult " Returns the matrix resulting from the matrix multiplication of the the given arguments. Equivalent to R's %*% operator. Examples: - - (def A (matrix [[1 2 3] - [4 5 6] + + (def A (matrix [[1 2 3] + [4 5 6] [7 8 9]])) (mmult A (trans A)) (mmult A (trans A) A) @@ -630,7 +630,7 @@ result (Matrix. (.zMult #^Matrix a #^Matrix b nil))] (if (and (= (.rows result) 1) (= (.columns result) 1)) (.getQuick result 0 0) - result))) + result))) args))) @@ -649,12 +649,12 @@ " ([& args] (reduce (fn [A B] - (let [a (cond - (matrix? A) A + (let [a (cond + (matrix? A) A (number? A) (matrix [A]) :else (matrix A)) - b (cond - (matrix? B) B + b (cond + (matrix? B) B (number? B) (matrix [B]) :else (matrix B)) rows (* (nrow a) (nrow b)) @@ -664,19 +664,19 @@ (mult (sel a i j) b))))))) args))) - -(defn solve + +(defn solve " Returns a matrix solution if A is square, least squares solution otherwise. Equivalent to R's solve function. - + Examples: (solve (matrix [[2 0 0] [0 2 0] [0 0 2]])) References: http://en.wikipedia.org/wiki/Matrix_inverse - + " ([#^Matrix A & B] @@ -686,7 +686,7 @@ -(defn det +(defn det " Returns the determinant of the given matrix using LU decomposition. Equivalent to R's det function. @@ -698,7 +698,7 @@ ([mat] (.det DoubleAlgebra/DEFAULT mat))) -(defn trace +(defn trace " Returns the trace of the given matrix. References: @@ -709,17 +709,17 @@ -(defn vectorize +(defn vectorize " Returns the vectorization (i.e. vec) of the given matrix. - The vectorization of an m-by-n matrix A, denoted by vec(A), - is the m*n-by-1 column vector obtain by stacking the columns + The vectorization of an m-by-n matrix A, denoted by vec(A), + is the m*n-by-1 column vector obtain by stacking the columns of the matrix A on top of one another. For instance: (= (vectorize (matrix [[a b] [c d]])) (matrix [a c b d])) Examples: - (def A (matrix [[1 2] [3 4]])) + (def A (matrix [[1 2] [3 4]])) (vectorize A) References: @@ -730,17 +730,17 @@ ;(reduce #(concat %1 (to-list %2)) '() (trans mat)))) -(defn half-vectorize +(defn half-vectorize " Returns the half-vectorization (i.e. vech) of the given matrix. - The half-vectorization, vech(A), of a symmetric nxn matrix A - is the n(n+1)/2 x 1 column vector obtained by vectorizing only + The half-vectorization, vech(A), of a symmetric nxn matrix A + is the n(n+1)/2 x 1 column vector obtained by vectorizing only the upper triangular part of A. For instance: (= (half-vectorize (matrix [[a b] [b d]])) (matrix [a b d])) Examples: - (def A (matrix [[1 2] [2 4]])) + (def A (matrix [[1 2] [2 4]])) (half-vectorize A) References: @@ -751,21 +751,21 @@ -(defn sum-of-squares +(defn sum-of-squares "Returns the sum-of-squares of the given sequence." ([x] (let [xx (if (or (nil? x) (empty? x)) [0] (to-list x))] (DoubleDescriptive/sumOfSquares (DoubleArrayList. (double-array xx)))))) -(defn sum +(defn sum "Returns the sum of the given sequence." ([x] (let [xx (if (or (nil? x) (empty? x)) [0] (to-list x))] (DoubleDescriptive/sum (DoubleArrayList. (double-array xx)))))) -(defn prod +(defn prod "Returns the product of the given sequence." ([x] (let [xx (if (or (nil? x) (empty? x)) [0] (to-list x))] @@ -773,8 +773,8 @@ -(defn cumulative-sum - " Returns a sequence of cumulative sum for the given collection. For instance +(defn cumulative-sum + " Returns a sequence of cumulative sum for the given collection. For instance The first value equals the first value of the argument, the second value is the sum of the first two arguments, the third is the sum of the first three arguments, etc. @@ -783,7 +783,7 @@ (use 'incanter.core) (cumulative-sum (range 100)) " - ([coll] + ([coll] (let [n (count coll)] (loop [in-coll (rest coll) cumu-sum [(first coll)]] @@ -806,10 +806,10 @@ Returns: a matrix of the triangular factor (note: the result from - cern.colt.matrix.linalg.CholeskyDecomposition is transposed so + cern.colt.matrix.linalg.CholeskyDecomposition is transposed so that it matches the result return from R's chol function. - + Examples: @@ -827,7 +827,7 @@ " ([#^Matrix mat] (.viewDice (.getL (DoubleCholeskyDecomposition. mat))))) - ;(Matrix. (.viewDice (.getL (CholeskyDecomposition. mat))))) + ;(Matrix. (.viewDice (.getL (CholeskyDecomposition. mat))))) @@ -847,7 +847,7 @@ (use 'incanter.core) (def foo (matrix (range 9) 3)) (decomp-foo foo) - + References: http://en.wikipedia.org/wiki/Singular_value_decomposition @@ -949,7 +949,7 @@ (use 'incanter.core) (def foo (matrix (range 9) 3)) (condition foo) - + References: http://en.wikipedia.org/wiki/Condition_number @@ -968,7 +968,7 @@ (use 'incanter.core) (def foo (matrix (range 9) 3)) (rank foo) - + References: @@ -985,7 +985,7 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn to-vect +(defn to-vect " Returns a vector-of-vectors if the given matrix is two-dimensional, and a flat vector if the matrix is one-dimensional. This is a bit slower than the to-list function. " @@ -999,15 +999,15 @@ (map #(into [] (seq %)) (seq (.toArray mat))))))) -(defn length - " A version of count that works on collections, matrices, and numbers. +(defn length + " A version of count that works on collections, matrices, and numbers. The length of a number is one, the length of a collection is its count, and the length of a matrix is the number of elements it contains (nrow*ncol). Equivalent to R's length function. " ([coll] (cond - (number? coll) + (number? coll) 1 (matrix? coll) (* (.rows #^Matrix coll) (.columns #^Matrix coll)) @@ -1015,11 +1015,11 @@ (count coll) :else (throw (Exception. "Argument must be a collection or matrix!"))))) - + (defn group-by -" Groups the given matrix by the values in the columns indicated by the +" Groups the given matrix by the values in the columns indicated by the 'on-cols' argument, returning a sequence of matrices. The returned matrices are sorted by the value of the group column ONLY when there is only a single (non-vector) on-col argument. @@ -1037,7 +1037,7 @@ (def plant-growth-dummies (to-matrix (get-dataset :plant-growth) :dummies true)) (group-by plant-growth-dummies [1 2]) ;; return only the first column - (group-by plant-growth-dummies [1 2] :cols 0) + (group-by plant-growth-dummies [1 2] :cols 0) ;; don't return the last two columns (group-by plant-growth-dummies [1 2] :except-cols [1 2]) @@ -1059,17 +1059,17 @@ groups (if (coll? on-cols) (into #{} (to-list (sel mat :cols on-cols))) (sort (into #{} (to-list (sel mat :cols on-cols))))) - filter-fn (fn [group] - (cond + filter-fn (fn [group] + (cond (and (coll? on-cols) (> (count on-cols) 1)) - (fn [row] - (reduce #(and %1 %2) + (fn [row] + (reduce #(and %1 %2) (map (fn [i g] (= (nth row i) g)) on-cols group))) (and (coll? on-cols) (= (count on-cols) 1)) - (fn [row] + (fn [row] (= (nth row (first on-cols)) group)) :else - (fn [row] + (fn [row] (= (nth row on-cols) group)))) ] (cond @@ -1087,12 +1087,12 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn dataset +(defn dataset " Returns a map of type ::dataset constructed from the given column-names and data. The data is a sequence of sequences. " - ([column-names & data] - (with-meta + ([column-names & data] + (with-meta {:column-names column-names :rows (map #(apply assoc {} (interleave column-names %)) (first data))} {:type ::dataset}))) @@ -1118,22 +1118,22 @@ ([dataset & options] (let [opts (when options (apply assoc {} options)) rows (or (:rows opts) true) - cols (if (:cols opts) + cols (if (:cols opts) (if (coll? (:cols opts)) - (:cols opts) + (:cols opts) [(:cols opts)]) (:column-names dataset)) row-filter (:filter opts) - selected-rows (cond - (true? rows) (:rows dataset) + selected-rows (cond + (true? rows) (:rows dataset) (number? rows) (list (nth (:rows dataset) rows)) - (coll? rows) (map #(nth (:rows dataset) %) rows)) + (coll? rows) (map #(nth (:rows dataset) %) rows)) data (map (fn [row] (map #(row (get-column-id dataset %)) cols)) selected-rows) result (if (nil? row-filter) data (filter row-filter data))] - (if (= (count cols) 1) + (if (= (count cols) 1) (mapcat identity result) - (with-meta (hash-map :column-names cols - :rows (map #(apply assoc {} (interleave cols %)) result)) + (with-meta (hash-map :column-names cols + :rows (map #(apply assoc {} (interleave cols %)) result)) {:type ::dataset}))))) @@ -1162,7 +1162,7 @@ (let [opts (when args (apply assoc {} args)) data (:data opts) ordered? (if (false? (:ordered? opts)) true false) - labels (or (:labels opts) + labels (or (:labels opts) (if (nil? data) (:levels opts) (sort (into #{} data)))) @@ -1174,7 +1174,7 @@ :to-levels (apply assoc {} (interleave labels levels))}))) -(defn to-levels +(defn to-levels " " ([coll & options] @@ -1184,10 +1184,10 @@ (for [label coll] (to-levels label))))) -(defn to-labels +(defn to-labels " " - ([coll cat-var] + ([coll cat-var] (let [to-labels (:to-labels cat-var)] (for [level coll] (to-labels level))))) @@ -1214,27 +1214,27 @@ (defn- string-to-categorical [dataset column-key dummies?] (let [col (first (get-columns dataset [column-key]))] - (if (some string? col) + (if (some string? col) (if dummies? (matrix (to-dummies col)) (matrix (to-levels col))) (matrix col)))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn to-matrix +(defn to-matrix " Converts a dataset into a matrix. Equivalent to R's as.matrix function for datasets. Options: - :dummies (default false) -- if true converts non-numeric variables into sets + :dummies (default false) -- if true converts non-numeric variables into sets of binary dummy variables, otherwise converts them into numeric codes. " ([dataset & options] (let [opts (when options (apply assoc {} options)) dummies? (if (true? (:dummies opts)) true false)] - (reduce bind-columns - (map #(string-to-categorical dataset % dummies?) + (reduce bind-columns + (map #(string-to-categorical dataset % dummies?) (range (count (keys (:column-names dataset))))))))) @@ -1248,7 +1248,7 @@ ;; GAMMA BASED FUNCTIONS FUNCTIONS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn gamma +(defn gamma " Equivalent to R's gamma function. @@ -1258,7 +1258,7 @@ ([x] (Gamma/gamma x))) -(defn beta +(defn beta " Equivalent to R's beta function. @@ -1268,7 +1268,7 @@ ([a b] (Gamma/beta a b))) -(defn incomplete-beta +(defn incomplete-beta " Returns the non-regularized incomplete beta value. @@ -1280,7 +1280,7 @@ -(defn regularized-beta +(defn regularized-beta " Returns the regularized incomplete beta value. Equivalent to R's pbeta function. @@ -1289,17 +1289,17 @@ http://en.wikipedia.org/wiki/Regularized_incomplete_beta_function http://mathworld.wolfram.com/RegularizedBetaFunction.html " - ([x a b] + ([x a b] (Gamma/incompleteBeta a b x))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; SYMMETRIC MATRIX +;; SYMMETRIC MATRIX ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn solve-quadratic -" +" Returns a vector with the solution to x from the quadratic equation, a*x^2 + b*x + c. @@ -1316,7 +1316,7 @@ http://en.wikipedia.org/wiki/Quadratic_formula " - ([a b c] + ([a b c] (let [t1 (- 0 b) t2 (sqrt (- (* b b) (* 4 a c))) t3 (* 2 a)] @@ -1335,15 +1335,15 @@ :lower (default true) -- lower-triangular. Set :lower to false to reverse the half-vectorize function. Examples: - + (use 'incanter.core) (symmetric-matrix [1 2 3 4 5 6 7 8 9 10]) - - (half-vectorize + + (half-vectorize (symmetric-matrix [1 2 3 4 5 6 @@ -1364,18 +1364,18 @@ (let [[i j] (nth indices idx)] (.set mat i j (nth data idx)) (.set mat j i (nth data idx)))) - mat))) + mat))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; VIEW METHODS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defmulti view - " This is a general 'view' function. When given an Incanter matrix/dataset - or a Clojure numeric collection, it will display it in a Java Swing - JTable. When given an Incanter chart object, it will display it in a new - window. When given a URL string, it will open the location with the +(defmulti view + " This is a general 'view' function. When given an Incanter matrix/dataset + or a Clojure numeric collection, it will display it in a Java Swing + JTable. When given an Incanter chart object, it will display it in a new + window. When given a URL string, it will open the location with the platform's default web browser. Examples: @@ -1411,8 +1411,8 @@ " - (fn [obj & options] (cond - (and (not (matrix? obj)) + (fn [obj & options] (cond + (and (not (matrix? obj)) (not (dataset? obj)) (coll? obj)) ::coll @@ -1426,16 +1426,16 @@ (defmethod view incanter.Matrix - ([obj & options] + ([obj & options] (let [opts (when options (apply assoc {} options)) column-names (or (:column-names opts) (range (ncol obj))) m (ncol obj) n (nrow obj)] (doto (JFrame. "Incanter Matrix") - (.add (JScrollPane. - (JTable. + (.add (JScrollPane. + (JTable. (cond - (and (> m 1) (> n 1)) + (and (> m 1) (> n 1)) (Vector. (map #(Vector. %) (to-list obj))) (or (and (> m 1) (= n 1)) (and (= m 1) (= n 1))) (Vector. (map #(Vector. %) [(to-list obj) []])) @@ -1447,11 +1447,11 @@ (defmethod view :incanter.core/dataset - ([obj & options] + ([obj & options] (let [column-names (:column-names obj) column-vals (map (fn [row] (map #(row %) column-names)) (:rows obj))] (doto (JFrame. "Incanter Dataset") - (.add (JScrollPane. (JTable. (Vector. (map #(Vector. %) column-vals)) + (.add (JScrollPane. (JTable. (Vector. (map #(Vector. %) column-vals)) (Vector. column-names)))) (.setSize 400 600) (.setVisible true))))) @@ -1459,26 +1459,26 @@ -;; URL view method code lifted from clojure.contrib.javadoc.browse/open-url-in-browser +;; URL view method code lifted from clojure.contrib.javadoc.browse/open-url-in-browser (defmethod view java.lang.String ([url] - (try + (try (when (clojure.lang.Reflector/invokeStaticMethod "java.awt.Desktop" "isDesktopSupported" (to-array nil)) (-> (clojure.lang.Reflector/invokeStaticMethod "java.awt.Desktop" "getDesktop" (to-array nil)) (.browse (java.net.URI. url))) url) - (catch ClassNotFoundException e nil)))) + (catch ClassNotFoundException e nil)))) (defn quit -" Exits the Clojure shell." +" Exits the Clojure shell." ([] (System/exit 0))) -(defmulti save -" Save is a multi-function that is used to write matrices, datasets and +(defmulti save +" Save is a multi-function that is used to write matrices, datasets and charts (in png format) to a file. Arguments: @@ -1486,11 +1486,11 @@ filename -- the filename to create. Matrix and dataset options: - :delim (default \\,) column delimiter + :delim (default \\,) column delimiter :header (default nil) an sequence of strings to be used as header line, for matrices the default value is nil, for datasets, the default is the dataset's column-names array. - :append (default false) determines whether this given file should be + :append (default false) determines whether this given file should be appended to. If true, a header will not be written to the file again. Chart options: @@ -1504,9 +1504,9 @@ (def A (matrix (range 12) 3)) ; creates a 3x4 matrix (save A \"A.dat\") ; writes A to the file A.dat, with no header and comma delimited (save A \"A.dat\" :delim \\tab) ; writes A to the file A.dat, with no header and tab delimited - + ;; writes A to the file A.dat, with a header and tab delimited - (save A \"A.dat\" :delim \\, :header [\"col1\" \"col2\" \"col3\"]) + (save A \"A.dat\" :delim \\, :header [\"col1\" \"col2\" \"col3\"]) Dataset Example: @@ -1521,13 +1521,10 @@ (use '(incanter core io stats charts)) (save (histogram (sample-normal 1000)) \"hist.png\") - + " - (fn [obj filename & options] + (fn [obj filename & options] (if (.isInstance processing.core.PApplet obj) :sketch (type obj)))) - - - src/incanter/core.clj @@ -16,13 +16,13 @@ -(ns incanter.datasets +(ns incanter.datasets (:use [incanter.io :only (read-dataset)])) -(def **datasets** - { +(def **datasets** + { :iris {:filename "data/iris.dat" :delim \space :header true} @@ -100,15 +100,15 @@ Datasets: :iris -- the Fisher's or Anderson's Iris data set gives the - measurements in centimeters of the variables sepal - length and width and petal length and width, - respectively, for 50 flowers from each of 3 species + measurements in centimeters of the variables sepal + length and width and petal length and width, + respectively, for 50 flowers from each of 3 species of iris. - :cars -- The data give the speed of cars and the distances taken + :cars -- The data give the speed of cars and the distances taken to stop. Note that the data were recorded in the 1920s. - :survey -- survey data used in Scott Lynch's 'Introduction to Applied Bayesian Statistics + :survey -- survey data used in Scott Lynch's 'Introduction to Applied Bayesian Statistics and Estimation for Social Scientists' :us-arrests -- This data set contains statistics, in arrests per 100,000 @@ -118,34 +118,34 @@ :flow-meter -- flow meter data used in Bland Altman Lancet paper. - :co2 -- has 84 rows and 5 columns of data from an experiment on the cold tolerance + :co2 -- has 84 rows and 5 columns of data from an experiment on the cold tolerance of the grass species _Echinochloa crus-galli_. - :chick-weight -- has 578 rows and 4 columns from an experiment on the effect of diet + :chick-weight -- has 578 rows and 4 columns from an experiment on the effect of diet on early growth of chicks. :plant-growth -- Results from an experiment to compare yields (as measured by dried weight of plants) obtained under a control and two different treatment conditions. - :pontius -- These data are from a NIST study involving calibration of load cells. - The response variable (y) is the deflection and the predictor variable + :pontius -- These data are from a NIST study involving calibration of load cells. + The response variable (y) is the deflection and the predictor variable (x) is load. See http://www.itl.nist.gov/div898/strd/lls/data/Pontius.shtml - :filip -- NIST data set for linear regression certification, + :filip -- NIST data set for linear regression certification, see http://www.itl.nist.gov/div898/strd/lls/data/Filip.shtml - :longely -- This classic dataset of labor statistics was one of the first used to - test the accuracy of least squares computations. The response variable - (y) is the Total Derived Employment and the predictor variables are GNP - Implicit Price Deflator with Year 1954 = 100 (x1), Gross National Product - (x2), Unemployment (x3), Size of Armed Forces (x4), Non-Institutional + :longely -- This classic dataset of labor statistics was one of the first used to + test the accuracy of least squares computations. The response variable + (y) is the Total Derived Employment and the predictor variables are GNP + Implicit Price Deflator with Year 1954 = 100 (x1), Gross National Product + (x2), Unemployment (x3), Size of Armed Forces (x4), Non-Institutional Population Age 14 & Over (x5), and Year (x6). See http://www.itl.nist.gov/div898/strd/lls/data/Longley.shtml - :Chwirut -- These data are the result of a NIST study involving ultrasonic calibration. - The response variable is ultrasonic response, and the predictor variable is + :Chwirut -- These data are the result of a NIST study involving ultrasonic calibration. + The response variable is ultrasonic response, and the predictor variable is metal distance. See http://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Chwirut1.dat @@ -157,7 +157,7 @@ :airline-passengers -- Monthly Airline Passenger Numbers 1949-1960 - :math-prog -- Pass/fail results for a high school mathematics assessment test + :math-prog -- Pass/fail results for a high school mathematics assessment test and a freshmen college programming course. :iran-election -- Vote counts for 30 provinces from the 2009 Iranian election. @@ -171,6 +171,3 @@ delim (ds :delim) header (ds :header)] (read-dataset filename :delim delim :header header)))) - - - src/incanter/datasets.clj @@ -50,13 +50,13 @@ (.. vv getRenderContext (setVertexLabelTransformer (new ToStringLabeller))) (.. vv getRenderContext (setEdgeLabelTransformer (new ToStringLabeller))) -(.. vv getRenderer - getVertexLabelRenderer +(.. vv getRenderer + getVertexLabelRenderer (setPosition (edu.uci.ics.jung.visualization.renderers.Renderer$VertexLabel$Position/CNTR))) (def frame (new JFrame "Simple Graph View")) ;(.setDefaultCloseOperation frame JFrame/EXIT_ON_CLOSE) ; this option kills the REPL -(.. frame getContentPane (add vv)) +(.. frame getContentPane (add vv)) (.pack frame) (.setVisible frame true) @@ -87,10 +87,10 @@ (.. vv getRenderContext (setEdgeStrokeTransformer edge-stroke-transformer)) (.. vv getRenderContext (setVertexLabelTransformer (new ToStringLabeller))) (.. vv getRenderContext (setEdgeLabelTransformer (new ToStringLabeller))) - ;(.. vv getRenderer getVertexLabelRenderer + ;(.. vv getRenderer getVertexLabelRenderer ;(setPosition (edu.uci.ics.jung.visualization.renderers.Renderer$VertexLabel$Position/CNTR))) ;(.setDefaultCloseOperation frame JFrame/EXIT_ON_CLOSE) ; this option kills the REPL - (.. frame getContentPane (add vv)) + (.. frame getContentPane (add vv)) (.pack frame) (.setVisible frame true)))) @@ -185,4 +185,3 @@ (. org.jgraph.graph.GraphConstants (setDashPattern (.getAttributes edge) (float-array [10 10]))) (.cellsChanged model (to-array [edge])) ;(. model (cellsChanged (to-array [edge]))) - src/incanter/graphs.clj @@ -16,7 +16,7 @@ -(ns incanter.internal +(ns incanter.internal (:import (incanter Matrix) (cern.colt.matrix.tdouble.algo DoubleFormatter) (cern.jet.math.tdouble DoubleFunctions DoubleArithmetic) @@ -26,17 +26,17 @@ ;;(derive DoubleMatrix2D ::matrix) ; commented out to track down non-ISeq matrices (derive Matrix ::matrix) -(defn is-matrix +(defn is-matrix " Test if obj is 'derived' from ::matrix (e.g. class incanter.Matrix)." ([obj] (isa? (class obj) ::matrix))) -(defn make-matrix +(defn make-matrix ([data] - (cond - (coll? (first data)) + (cond + (coll? (first data)) (Matrix. (into-array (map double-array data))) - (number? (first data)) + (number? (first data)) (Matrix. (double-array data)))) ([data ncol] (cond @@ -51,60 +51,60 @@ (defmacro #^Matrix transform-with [A op fun] - `(cond + `(cond (is-matrix ~A) (.assign (.copy ~A) (. DoubleFunctions ~fun)) - (and (coll? ~A) (coll? (first ~A))) + (and (coll? ~A) (coll? (first ~A))) (.assign #^Matrix (make-matrix ~A) (. DoubleFunctions ~fun)) (coll? ~A) (map ~op ~A) (number? ~A) (~op ~A))) - + (defmacro combine-with [A B op fun] - `(cond + `(cond (and (number? ~A) (number? ~B)) (~op ~A ~B) (and (is-matrix ~A) (is-matrix ~B)) - (.assign #^Matrix (.copy #^Matrix ~A) - #^Matrix ~B + (.assign #^Matrix (.copy #^Matrix ~A) + #^Matrix ~B #^DoubleDoubleFunction (. DoubleFunctions ~fun)) (and (is-matrix ~A) (number? ~B)) - (.assign #^Matrix (.copy #^Matrix ~A) + (.assign #^Matrix (.copy #^Matrix ~A) (make-matrix ~B (.rows ~A) (.columns ~A)) #^DoubleDoubleFunction (. DoubleFunctions ~fun)) ;;#^DoubleDoubleFunction (. DoubleFunctions (~fun ~B))) (and (number? ~A) (is-matrix ~B)) - (.assign #^Matrix (make-matrix ~A (.rows ~B) (.columns ~B)) - #^Matrix ~B + (.assign #^Matrix (make-matrix ~A (.rows ~B) (.columns ~B)) + #^Matrix ~B #^DoubleDoubleFunction (. DoubleFunctions ~fun)) (and (coll? ~A) (is-matrix ~B)) - (.assign #^Matrix (make-matrix ~A (.columns ~B)) - #^Matrix (make-matrix ~B) + (.assign #^Matrix (make-matrix ~A (.columns ~B)) + #^Matrix (make-matrix ~B) #^DoubleDoubleFunction (. DoubleFunctions ~fun)) (and (is-matrix ~A) (coll? ~B)) (.assign #^Matrix (.copy ~A) - #^Matrix (make-matrix ~B) + #^Matrix (make-matrix ~B) #^DoubleDoubleFunction (. DoubleFunctions ~fun)) - (and (coll? ~A) (coll? ~B) (coll? (first ~A))) - (.assign (make-matrix ~A) + (and (coll? ~A) (coll? ~B) (coll? (first ~A))) + (.assign (make-matrix ~A) (make-matrix ~B) (. DoubleFunctions ~fun)) - (and (coll? ~A) (number? ~B) (coll? (first ~A))) - (.assign (make-matrix ~A) - (make-matrix ~B) + (and (coll? ~A) (number? ~B) (coll? (first ~A))) + (.assign (make-matrix ~A) + (make-matrix ~B) (. DoubleFunctions ~fun)) ;;(. DoubleFunctions (~fun ~B))) - (and (number? ~A) (coll? ~B) (coll? (first ~B))) - (.assign (make-matrix ~A (.rows ~B) (.columns ~B)) + (and (number? ~A) (coll? ~B) (coll? (first ~B))) + (.assign (make-matrix ~A (.rows ~B) (.columns ~B)) (make-matrix ~B) (. DoubleFunctions ~fun)) (and (coll? ~A) (coll? ~B)) - (map ~op ~A ~B) - (and (number? ~A) (coll? ~B)) + (map ~op ~A ~B) + (and (number? ~A) (coll? ~B)) (map ~op (replicate (count ~B) ~A) ~B) - (and (coll? ~A) (number? ~B)) + (and (coll? ~A) (number? ~B)) (map ~op ~A (replicate (count ~A) ~B)) )) @@ -112,7 +112,7 @@ ;; PRINT METHOD FOR COLT MATRICES (defmethod print-method Matrix [o, #^java.io.Writer w] (let [formatter (DoubleFormatter. "%1.4f")] - (do + (do (.setPrintShape formatter false) (.write w "[") (.write w (.toString formatter o)) @@ -121,12 +121,9 @@ ;; PRINT METHOD FOR INCANTER DATASETS (defmethod print-method :incanter.core/dataset [o, #^java.io.Writer w] - (do + (do (.write w (str (:column-names o))) (.write w "\n") (doseq [row (:rows o)] (.write w (str (apply vector (map #(get row %) (:column-names o))))) (.write w "\n")))) - - - src/incanter/internal.clj @@ -19,21 +19,21 @@ ;; DATA IO FUNCTIONS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(ns incanter.io +(ns incanter.io (:import (java.io FileReader) (au.com.bytecode.opencsv CSVReader)) (:use [incanter.core :only (dataset save)]) ) -(defn- parse-string [value] - (try (Integer/parseInt value) - (catch NumberFormatException _ +(defn- parse-string [value] + (try (Integer/parseInt value) + (catch NumberFormatException _ (try (Double/parseDouble value) (catch NumberFormatException _ value))))) -(defn read-dataset +(defn read-dataset " Returns a dataset read from a file. @@ -43,14 +43,14 @@ :skip (default 0) the number of lines to skip at the top of the file. :header (default false) indicates the file has a header line " - ([filename & options] + ([filename & options] (let [opts (when options (apply assoc {} options)) delim (or (:delim opts) \,) ; space delim default quote-char (or (:quote opts) \") skip (or (:skip opts) 0) header? (or (:header opts) false) - reader (au.com.bytecode.opencsv.CSVReader. - (java.io.FileReader. filename) + reader (au.com.bytecode.opencsv.CSVReader. + (java.io.FileReader. filename) delim quote-char skip) @@ -59,21 +59,21 @@ parsed-data (into [] (map (fn [row] (into [] (map #(parse-string %) row))) raw-data)) ] ;;(if header? (dataset (first parsed-data) (rest parsed-data) (dataset parsed-data)))))) - (if header? + (if header? ; have header row - (dataset (first parsed-data) (rest parsed-data)) + (dataset (first parsed-data) (rest parsed-data)) ; no header row so build a default one (let [col-count (count (first parsed-data)) col-names (apply vector (map str (repeat col-count "col") (iterate inc 0)))] (dataset col-names parsed-data)))))) - + (defmethod save incanter.Matrix [mat filename & options] (let [opts (when options (apply assoc {} options)) - delim (or (:delim opts) \,) - header (or (:header opts) nil) + delim (or (:delim opts) \,) + header (or (:header opts) nil) append? (if (true? (:append opts)) true false) file-writer (java.io.FileWriter. filename append?)] (do @@ -92,13 +92,13 @@ (.write file-writer (str \newline))))) (.flush file-writer) (.close file-writer)))) - + (defmethod save :incanter.core/dataset [dataset filename & options] (let [opts (when options (apply assoc {} options)) - delim (or (:delim opts) \,) + delim (or (:delim opts) \,) header (or (:header opts) (:column-names dataset)) append? (if (true? (:append opts)) true false) file-writer (java.io.FileWriter. filename append?) @@ -118,4 +118,3 @@ (.write file-writer (str \newline)))) (.flush file-writer) (.close file-writer)))) - src/incanter/io.clj @@ -16,13 +16,13 @@ -(ns incanter.optimize +(ns incanter.optimize (:use [incanter.core :only (plus minus div mult mmult symmetric-matrix ncol solve abs sel trans bind-columns to-list identity-matrix)])) -(defn integrate +(defn integrate " Integrate a function f from a to b @@ -39,7 +39,7 @@ (def std 1) (def mu 0) (defn normal [x] - (/ 1 + (/ 1 (* (* std (Math/sqrt (* 2 Math/PI))) (Math/exp (/ (Math/pow (- (- x mu)) 2) (* 2 (Math/pow std 2))))))) @@ -64,7 +64,7 @@ -(defn derivative +(defn derivative " Returns a function that approximates the derivative of the given function. @@ -88,7 +88,7 @@ ;; get the second derivative function (def cube-deriv2 (derivative cube-deriv)) (add-lines plot x (map cube-deriv2 x)) - + ;; plot the normal pdf and its derivatives (def plot (xy-plot x (pdf-normal x))) (view plot) @@ -115,9 +115,9 @@ Examples: (defn quad-fx [[x y]] (+ (* x x) (* x y) (* y y))) - (def quad-dfx0 (partial-derivative quad-fx 0)) + (def quad-dfx0 (partial-derivative quad-fx 0)) (def quad-dfx1 (partial-derivative quad-fx 1)) - (quad-dfx0 [1 1]) + (quad-dfx0 [1 1]) (use '(incanter core optimize charts)) (def x (range -3 3 0.1)) @@ -127,12 +127,12 @@ References: http://en.wikipedia.org/wiki/Partial_derivative - + " ([fx i & options] (let [opts (when options (apply assoc {} options)) dx (or (:dx opts) 0.0001)] - (fn [theta] + (fn [theta] (let [theta-next (assoc theta i (+ (theta i) dx))] (/ (- (fx theta-next) (fx theta)) dx)))))) @@ -161,7 +161,7 @@ (defn- second-partial-derivative " Examples: - + (use '(incanter core optimize charts)) (defn quad-fx [[x y]] (+ (* x x) (* x y) (* y y))) (def quad-dfx00 (second-partial-derivative quad-fx 0 0)) @@ -172,7 +172,7 @@ (quad-dfx10 [1 1]) (quad-dfx01 [1 1]) (quad-dfx11 [1 1]) - + " ([fx i j] (partial-derivative (partial-derivative fx i) j))) @@ -194,7 +194,7 @@ " ([fx n & options] - (let [funs (for [i (range n) j (range n) :when (<= j i)] + (let [funs (for [i (range n) j (range n) :when (<= j i)] (second-partial-derivative fx i j))] (fn [theta] (symmetric-matrix (map #(% theta) funs)))))) @@ -211,8 +211,8 @@ Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -240,7 +240,7 @@ (reduce bind-columns (for [i (range p)] (div (map - (map (partial f (plus theta (mult (nth e i) (nth dx i)))) x) - (map (partial f theta) x)) + (map (partial f theta) x)) (nth dx i)))))))) @@ -248,19 +248,19 @@ (defn gradient -" - Returns a function that calculates a 5-point approximation to +" + Returns a function that calculates a 5-point approximation to the gradient of the given function. The vector of start values are used to determine the number of parameters required by the function, and - to scale the step-size. The generated function accepts a vector of + to scale the step-size. The generated function accepts a vector of parameter values and a vector of x data points and returns a matrix, where each row is the gradient evaluated at the corresponding x value. Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -286,21 +286,21 @@ (reduce bind-columns (for [i (range p)] (let [h (mult (nth e i) dx)] - (div - (map + (map (partial f (minus theta (mult 2 h))) x) - (map - (mult 8 (map (partial f (minus theta h)) x))) - (mult 8 (map (partial f (plus theta h)) x)) - (map - (map (partial f (plus theta (mult 2 h))) x))) + (div + (map + (map (partial f (minus theta (mult 2 h))) x) + (map - (mult 8 (map (partial f (minus theta h)) x))) + (mult 8 (map (partial f (plus theta h)) x)) + (map - (map (partial f (plus theta (mult 2 h))) x))) (* 12 (nth dx i)))))))))) (defn hessian -" Returns a function that calculates an approximation to the Hessian matrix - of the given function. The vector of start values are used to determine - the number of parameters required by the function, and to scale the - step-size. The generated function accepts a vector of +" Returns a function that calculates an approximation to the Hessian matrix + of the given function. The vector of start values are used to determine + the number of parameters required by the function, and to scale the + step-size. The generated function accepts a vector of parameter values and a vector of x data points and returns a matrix, where each row with p*(p+1)/2 columns, one for each unique entry in the Hessian evaluated at the corresponding x value. @@ -308,8 +308,8 @@ Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -337,11 +337,11 @@ (let [hi (mult (nth e i) dx) hj (mult (nth e j) dx) hij (mult (plus (nth e i) (nth e j)) dx)] - (div + (div (map + - (map - (map (partial f theta) x) - (map (partial f (plus theta hi)) x) - (map (partial f (plus theta hj)) x)) + (map - (map (partial f theta) x) + (map (partial f (plus theta hi)) x) + (map (partial f (plus theta hj)) x)) (map (partial f (plus theta hij)) x)) (* (nth dx i) (nth dx j)))))))))) @@ -352,7 +352,7 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ROUTINES FOR NON-LINEAR LEAST SQUARES -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn- nls-rss @@ -363,8 +363,8 @@ Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -396,8 +396,8 @@ Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -408,8 +408,8 @@ ;; view the data (view (scatter-plot x y)) - (nls-gradient f - (gradient f theta-init) + (nls-gradient f + (gradient f theta-init) theta-init x y) @@ -434,8 +434,8 @@ Examples: (use '(incanter core optimize datasets charts)) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -446,9 +446,9 @@ ;; view the data (view (scatter-plot x y)) - (time (doall (nls-hessian f - (gradient f theta-init) - (hessian f theta-init) + (time (doall (nls-hessian f + (gradient f theta-init) + (hessian f theta-init) theta-init x y))) @@ -479,8 +479,8 @@ (use '(incanter core optimize datasets charts)) ;; define the Michaelis-Menton model function ;; y = a + (b - a)*x/(c + x) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) @@ -497,12 +497,12 @@ (def d2f (hessian f start)) (def result (nls-newton-raphson f df d2f start x y)) - ;(add-lines plot x (f (:theta result) x)) - (add-lines plot x (map (partial f (:theta result)) x)) + ;(add-lines plot x (f (:theta result) x)) + (add-lines plot x (map (partial f (:theta result)) x)) -" +" ([f df d2f start x y & options] (let [opts (when options (apply assoc {} options)) max-iter (or (:max-iter opts) 200) @@ -513,9 +513,9 @@ (let [H (solve (nls-hessian f df d2f th x y)) g (nls-gradient f df th x y)] (if (or (< (reduce max (abs g)) tol) (= i max-iter)) - {:theta th - :iterations i - :gradient g + {:theta th + :iterations i + :gradient g :hessian H :rss (nls-rss f th x y)} (recur (inc i) (map - th (mmult H g))))))))) @@ -534,12 +534,12 @@ (use '(incanter core optimize datasets charts)) ;; define the Michaelis-Menton model function ;; y = a + (b - a)*x/(c + x) - (defn f [theta x] - (+ (nth theta 0) + (defn f [theta x] + (+ (nth theta 0) (div (* x (- (nth theta 1) (nth theta 0))) (+ (nth theta 2) x)))) -; (defn f [theta x] +; (defn f [theta x] ; (let [[a b c] theta] ; (plus a (div (mult x (minus b a)) (plus c x))))) @@ -554,12 +554,12 @@ ;(def df (gradient f start)) (def result (nls-gauss-newton f df start x y)) - ;(add-lines plot x (f (:theta result) x)) - (add-lines plot x (map (partial f (:theta result)) x)) + ;(add-lines plot x (f (:theta result) x)) + (add-lines plot x (map (partial f (:theta result)) x)) -" +" ([f df start x y & options] (let [opts (when options (apply assoc {} options)) max-iter (or (:max-iter opts) 200) @@ -574,8 +574,8 @@ update (mmult (solve (mmult (trans g) g)) (trans g) resid) ] (if (or (< (reduce max (abs update)) tol) (= i max-iter)) - {:theta th - :iterations i + {:theta th + :iterations i :gradient g :rss (nls-rss f th x y)} (recur (inc i) (map + th update)))))))) @@ -585,9 +585,9 @@ (defn non-linear-model -" +" Determine the nonlinear least-squares estimates of the - parameters of a nonlinear model. + parameters of a nonlinear model. Based on R's nls (non-linear least squares) function. Arguments: @@ -599,7 +599,7 @@ Options: :method (default :gauss-newton) other option :newton-raphson - :tol (default 1E-5) + :tol (default 1E-5) :max-iter (default 200) Returns: a hash-map containing the following fields: @@ -620,7 +620,7 @@ (use '(incanter core optimize datasets charts)) ;; define the Michaelis-Menton model function ;; y = a + (b - a)*x/(c + x) - (defn f [theta x] + (defn f [theta x] (let [[a b c] theta] (plus a (div (mult x (minus b a)) (plus c x))))) @@ -651,9 +651,9 @@ (def plot (scatter-plot x y :legend true)) (view plot) - - ;; the newton-raphson algorithm fails to converge to the correct solution - ;; using first set of start values from NIST, but the default gauss-newton + + ;; the newton-raphson algorithm fails to converge to the correct solution + ;; using first set of start values from NIST, but the default gauss-newton ;; algorith converges to the correct solution. (def start1 [0.1 0.01 0.02]) @@ -687,6 +687,3 @@ :fitted fitted :x x :y y}))) - - - src/incanter/optimize.clj @@ -7,7 +7,7 @@ ;; which can be found in the file CPL.TXT at the root of this ;; distribution. By using this software in any fashion, you are ;; agreeing to be bound by the terms of this license. You must not -;; remove this notice, or any other, from this software. +;; remove this notice, or any other, from this software. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Modified by David Edgar Liebke 30 June 2009 @@ -37,7 +37,7 @@ ;;(defn ambient -;; ([sketch rgb] +;; ([sketch rgb] ;; (if (string? rgb) ;; (.ambient sketch (color sketch rgb)) ;; (.ambient sketch (float rgb)))) @@ -55,7 +55,7 @@ ;; ([sketch rgb] ;; (.ambientLight sketch rgb)) ;; ([sketch rgb x y z] - (.ambientLight sketch rgb (float x) (float y) (float z)))) + (.ambientLight sketch rgb (float x) (float y) (float z)))) (defn ambient-light ([sketch red green blue] @@ -67,16 +67,16 @@ (defn apply-matrix ([sketch n00 n01 n02 n10 n11 n12] - (.applyMatrix sketch (float n00) (float n01) (float n02) - (float n10) (float n11) (float n12))) + (.applyMatrix sketch (float n00) (float n01) (float n02) + (float n10) (float n11) (float n12))) ([sketch n00 n01 n02 n03 - n10 n11 n12 n13 - n20 n21 n22 n23 - n30 n31 n32 n33] - (.applyMatrix sketch (float n00) (float n01) (float n02) (float 03) - (float n10) (float n11) (float n12) (float 13) - (float n20) (float n21) (float n22) (float 23) - (float n30) (float n31) (float n32) (float 33)))) + n10 n11 n12 n13 + n20 n21 n22 n23 + n30 n31 n32 n33] + (.applyMatrix sketch (float n00) (float n01) (float n02) (float 03) + (float n10) (float n11) (float n12) (float 13) + (float n20) (float n21) (float n22) (float 23) + (float n30) (float n31) (float n32) (float 33)))) (defn arc [sketch a b c d start stop] @@ -143,8 +143,8 @@ " ;; ([#^PApplet sketch rgb] (.background sketch rgb))) - ([#^PApplet sketch gray] - (cond + ([#^PApplet sketch gray] + (cond (string? gray) (.background sketch (color gray)) (integer? gray) @@ -153,22 +153,22 @@ (.background sketch (float gray)) (= java.awt.Color (type gray)) (.background sketch #^java.awt.Color gray))) - ([#^PApplet sketch gray alpha] - (cond + ([#^PApplet sketch gray alpha] + (cond (string? gray) (.background sketch (color gray) (int alpha)) (integer? gray) (.background sketch (int gray) (int alpha)) (float? gray) (.background sketch (float gray) (float alpha)))) - ([#^PApplet sketch r g b] - (cond + ([#^PApplet sketch r g b] + (cond (or (integer? r) (integer? g) (integer? b)) (.background sketch (int r) (int g) (int b)) (or (float? r) (float? g) (float? b)) (.background sketch (float r) (float g) (float b)))) - ([#^PApplet sketch r g b a] - (cond + ([#^PApplet sketch r g b a] + (cond (or (integer? r) (integer? g) (integer? b)) (.background sketch (int r) (int g) (int b) (int a)) (or (float? r) (float? g) (float? b)) @@ -177,13 +177,13 @@ ;; DEL ;;(defn background -;; ([#^PApplet sketch color] -;; (if (float? color) -;; (.background sketch (float color)) +;; ([#^PApplet sketch color] +;; (if (float? color) +;; (.background sketch (float color)) ;; (.background sketch (int color)))) -;; ([#^PApplet sketch color alpha] -;; (if (float? color) -;; (.background sketch (float color) (float alpha)) +;; ([#^PApplet sketch color alpha] +;; (if (float? color) +;; (.background sketch (float color) (float alpha)) ;; (.background sketch (int color) (float alpha)))) ;; ([#^PApplet sketch r g b] (.background sketch (float r) (float g) (float b))) ;; ([#^PApplet sketch r g b alpha] (.background sketch (float r) (float g) (float b) (float alpha)))) @@ -196,8 +196,8 @@ (defn begin-raw ([#^PApplet sketch #^PGraphics rawGfx] (.beginRaw sketch rawGfx)) - ([#^PApplet sketch #^java.lang.String renderer #^java.lang.String filename] - (.beginRaw sketch renderer filename))) + ([#^PApplet sketch #^java.lang.String renderer #^java.lang.String filename] + (.beginRaw sketch renderer filename))) ;; $$beginRecord @@ -207,17 +207,17 @@ (defn bezier ([#^PApplet sketch x1 y1 x2 y2 x3 y3 x4 y4] - (.bezier sketch - (float x1) (float y1) - (float x2) (float y2) - (float x3) (float y3) - (float x4) (float y4))) + (.bezier sketch + (float x1) (float y1) + (float x2) (float y2) + (float x3) (float y3) + (float x4) (float y4))) ([#^PApplet sketch x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4] - (.bezier sketch - (float x1) (float y1) (float z1) - (float x2) (float y2) (float z2) - (float x3) (float y3) (float z3) - (float x4) (float y4) (float z4)))) + (.bezier sketch + (float x1) (float y1) (float z1) + (float x2) (float y2) (float z2) + (float x3) (float y3) (float z3) + (float x4) (float y4) (float z4)))) (defn bezier-detail [#^PApplet sketch detail] (.bezierDetail sketch (int detail))) @@ -230,26 +230,26 @@ (defn bezier-vertex ([#^PApplet sketch x2 y2 x3 y3 x4 y4] - (.bezierVertex sketch - (float x2) (float y2) - (float x3) (float y3) - (float x4) (float y4))) + (.bezierVertex sketch + (float x2) (float y2) + (float x3) (float y3) + (float x4) (float y4))) ([#^PApplet sketch x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4] - (.bezierVertex sketch - (float x2) (float y2) (float z2) - (float x3) (float y3) (float z3) - (float x4) (float y4) (float z4)))) + (.bezierVertex sketch + (float x2) (float y2) (float z2) + (float x3) (float y3) (float z3) + (float x4) (float y4) (float z4)))) ;; $$binary (defn blend ([#^PApplet sketch sx1 sy1 sx2 sy2 dx1 dy1 dx2 dy2 mode] - (.blend sketch (int sx1) (int sy1) (int sx2) (int sy2) - (int dx1) (int dy1) (int dx2) (int dy2) + (.blend sketch (int sx1) (int sy1) (int sx2) (int sy2) + (int dx1) (int dy1) (int dx2) (int dy2) (int mode))) ([#^PApplet sketch #^PImage src sx1 sy1 sx2 sy2 dx1 dy1 dx2 dy2 mode] - (.blend sketch src (int sx1) (int sy1) (int sx2) (int sy2) - (int dx1) (int dy1) (int dx2) (int dy2) + (.blend sketch src (int sx1) (int sy1) (int sx2) (int sy2) + (int dx1) (int dy1) (int dx2) (int dy2) (int mode)))) (defn blend-color @@ -266,7 +266,7 @@ (defn camera ([#^PApplet sketch ] (.camera sketch)) ([#^PApplet sketch eyeX eyeY eyeZ centerX centerY centerZ upX upY upZ] - (.camera sketch (float eyeX) (float eyeY) (float eyeZ) (float centerX) (float centerY) (float centerZ) (float upX) (float upY) (float upZ)))) + (.camera sketch (float eyeX) (float eyeY) (float eyeZ) (float centerX) (float centerY) (float centerZ) (float upX) (float upY) (float upZ)))) (defn can-draw? [#^PApplet sketch ] (.canDraw sketch)) @@ -285,7 +285,7 @@ ;; ([#^PApplet sketch r g b a] (.background sketch (int r) (int g) (int b) (int a)))) -;; modified by DEL +;; modified by DEL (defn color " @@ -304,26 +304,26 @@ (color 1.0 0.0 1.0 1.0) ;; with alpha value " -([rgb] +([rgb] (.getRGB (java.awt.Color. (int rgb)))) -;; (cond +;; (cond ;; (string? rgb) ;; (.getRGB (java.awt.Color/decode rgb)) ;; (or (integer? rgb) (float? rgb)) ;; (.getRGB (java.awt.Color. (int rgb))))) -([rgb alpha?] +([rgb alpha?] ;;(.getRGB (java.awt.Color. (color rgb) alpha?))) (.getRGB (java.awt.Color. (int rgb) alpha?))) -([x y z] +([x y z] (if (or (float? x) (float? y) (float? z)) (.getRGB (java.awt.Color. (float x) (float y) (float z))) (.getRGB (java.awt.Color. (int x) (int y) (int z))))) -([x y z alpha] +([x y z alpha] (if (or (float? x) (float? y) (float? z)) (.getRGB (java.awt.Color. (float x) (float y) (float z) (float alpha))) (.getRGB (java.awt.Color. (int x) (int y) (int z) (int alpha)))))) -;;([#^PApplet sketch rgb] +;;([#^PApplet sketch rgb] ;; (if (string? rgb) ;; (.getRGB (java.awt.Color/decode rgb)) ;; (.color sketch (int rgb)))) @@ -333,54 +333,54 @@ -;; ([#^PApplet sketch gray] -;; (cond -;; (string? gray) +;; ([#^PApplet sketch gray] +;; (cond +;; (string? gray) ;; (.color sketch (Integer/parseInt gray 16)) -;; (integer? gray) +;; (integer? gray) ;; (.color sketch (int gray)) -;; (float? gray) +;; (float? gray) ;; (.color sketch (float gray)))) -;; ([#^PApplet sketch gray alpha] -;; (cond -;; (string? gray) +;; ([#^PApplet sketch gray alpha] +;; (cond +;; (string? gray) ;; (.color sketch (Integer/parseInt gray 16) (Integer/parseInt alpha 16)) -;; (integer? gray) +;; (integer? gray) ;; (.color sketch (int gray) (int alpha)) -;; (float? gray) +;; (float? gray) ;; (.color sketch (float gray) (int alpha)))) -;; ([#^PApplet sketch r g b] +;; ([#^PApplet sketch r g b] ;; (if (or (integer? r) (integer? g) (integer? b)) ;; (.color sketch (int r) (int g) (int b)) ;; (.color sketch (float r) (float g) (float b)))) -;; ([#^PApplet sketch r g b alpha] +;; ([#^PApplet sketch r g b alpha] ;; (if (or (integer? r) (integer? g) (integer? b)) ;; (.color sketch (int r) (int g) (int b) (int alpha)) ;; (.color sketch (float r) (float g) (float b) (float alpha))))) ;; -;; ([rgb] +;; ([rgb] ;; (if (string? rgb) ;; (.getRGB (java.awt.Color/decode rgb)) ;; (.getRGB (java.awt.Color. (int rgb))))) -; (let [rgb (cond -; (> rgb 255) 255 +; (let [rgb (cond +; (> rgb 255) 255 ; (< rgb 0) 0 ; :else rgb)] -; (reduce bit-or [(int 0xff000000) -; (bit-shift-left (int rgb) 16) -; (bit-shift-left (int rgb) 8) +; (reduce bit-or [(int 0xff000000) +; (bit-shift-left (int rgb) 16) +; (bit-shift-left (int rgb) 8) ; (int rgb)])))) - + (defn color-mode - ([#^PApplet sketch mode] + ([#^PApplet sketch mode] (.colorMode sketch (int mode))) - ([#^PApplet sketch mode max] + ([#^PApplet sketch mode max] (.colorMode sketch (int mode) (float max))) - ([#^PApplet sketch mode max-x max-y max-z] + ([#^PApplet sketch mode max-x max-y max-z] (.colorMode sketch (int mode) (float max-x) (float max-y) (float max-z))) - ([#^PApplet sketch mode max-x max-y max-z max-a] + ([#^PApplet sketch mode max-x max-y max-z max-a] (.colorMode sketch (int mode) (float max-x) (float max-y) (float max-z) (float max-a)))) ;; $$concat @@ -403,11 +403,11 @@ (defn copy-pixels " Processing copy function. " ([#^PApplet sketch [sx1 sy1 sx2 sy2] [dx1 dy1 dx2 dy2]] - (.copy sketch (int sx1) (int sy1) (int sx2) (int sy2) - (int dx1) (int dy1) (int dx2) (int dy2))) + (.copy sketch (int sx1) (int sy1) (int sx2) (int sy2) + (int dx1) (int dy1) (int dx2) (int dy2))) ([#^PApplet sketch #^PImage img [sx1 sy1 sx2 sy2] [dx1 dy1 dx2 dy2]] - (.copy sketch img (int sx1) (int sy1) (int sx2) (int sy2) - (int dx1) (int dy1) (int dx2) (int dy2)))) + (.copy sketch img (int sx1) (int sy1) (int sx2) (int sy2) + (int dx1) (int dy1) (int dx2) (int dy2)))) ;;(defn cos [angle] (PApplet/cos (float angle))) @@ -415,7 +415,7 @@ ([#^PApplet sketch name size] (.createFont sketch name (float size))) ([#^PApplet sketch name size smooth] (.createFont sketch name (float size) smooth)) ([#^PApplet sketch name size smooth #^chars charset] - (.createFont sketch name (float size) smooth charset))) + (.createFont sketch name (float size) smooth charset))) ;; added by DEL @@ -439,8 +439,8 @@ (export-font \"Ariel\" 48 \"/tmp/ariel_48.vlw\") (view - (sketch - (setup [] + (sketch + (setup [] (let [font (load-font this \"/tmp/ariel_48.vlw\")] (doto this (text-font font) @@ -467,22 +467,22 @@ (defn create-graphics ([#^PApplet sketch w h renderer] - (.createGraphics sketch (int w) (int h) renderer)) + (.createGraphics sketch (int w) (int h) renderer)) ([#^PApplet sketch w h renderer path] - (.createGraphics sketch (int w) (int h) renderer path))) + (.createGraphics sketch (int w) (int h) renderer path))) (defn create-image [#^PApplet sketch w h format] (.createImage sketch (int w) (int h) (int format))) (defn create-input [filename] - (PApplet/createInput (java.io.File. filename))) + (PApplet/createInput (java.io.File. filename))) (defn create-input-raw - "Call openStream() without automatic gzip decompression." - [#^PApplet sketch filename] - (.createInputRaw sketch filename)) + "Call openStream() without automatic gzip decompression." + [#^PApplet sketch filename] + (.createInputRaw sketch filename)) (defn create-output [filename] - (PApplet/createOutput (java.io.File. filename))) + (PApplet/createOutput (java.io.File. filename))) (defn create-path [filename] (PApplet/createPath filename)) @@ -500,17 +500,17 @@ (defn curve ([#^PApplet sketch x1 y1 x2 y2 x3 y3 x4 y4] - (.curve sketch - (float x1) (float y1) - (float x2) (float y2) - (float x3) (float y3) - (float x4) (float y4))) + (.curve sketch + (float x1) (float y1) + (float x2) (float y2) + (float x3) (float y3) + (float x4) (float y4))) ([#^PApplet sketch x1 y1 z1 x2 y2 z2 x3 y3 z3 x4 y4 z4] - (.curve sketch - (float x1) (float y1) (float z1) - (float x2) (float y2) (float z2) - (float x3) (float y3) (float z3) - (float x4) (float y4) (float z4)))) + (.curve sketch + (float x1) (float y1) (float z1) + (float x2) (float y2) (float z2) + (float x3) (float y3) (float z3) + (float x4) (float y4) (float z4)))) (defn curve-detail [#^PApplet sketch detail] (.curveDetail sketch (int detail))) @@ -562,11 +562,11 @@ ;; modified by DEL (defn emissive - ([#^PApplet sketch gray] + ([#^PApplet sketch gray] (if (integer? gray) (.emissive sketch (int gray)) (.emissive sketch (float gray)))) - ([#^PApplet sketch x y z] + ([#^PApplet sketch x y z] (if (integer? gray) (.emissive sketch (int x) (int y) (int z)) (.emissive sketch (float x) (float y) (float z))))) @@ -592,21 +592,21 @@ ;; modified by DEL (defn fill " - Sets the color used to fill shapes. For example, if you run (fill 204 102 0), all - subsequent shapes will be filled with orange. This color is either specified in - terms of the RGB or HSB color depending on the current colorMode() (the default - color space is RGB, with each value in the range from 0 to 255). + Sets the color used to fill shapes. For example, if you run (fill 204 102 0), all + subsequent shapes will be filled with orange. This color is either specified in + terms of the RGB or HSB color depending on the current colorMode() (the default + color space is RGB, with each value in the range from 0 to 255). - When using hexadecimal notation to specify a color, use '16r' before - the values (e.g. 16rCCFFAA, 16rFFCCFFAA). Use six digits to specify - a color (the way colors are specified in HTML and CSS). When eight - digits are used, the first two characters define the alpha component and the - remainder the red, green, and blue components. + When using hexadecimal notation to specify a color, use '16r' before + the values (e.g. 16rCCFFAA, 16rFFCCFFAA). Use six digits to specify + a color (the way colors are specified in HTML and CSS). When eight + digits are used, the first two characters define the alpha component and the + remainder the red, green, and blue components. - The value for the parameter \"gray\" must be less than or equal to the current - maximum value as specified by (colorMode). The default maximum value is 255. + The value for the parameter \"gray\" must be less than or equal to the current + maximum value as specified by (colorMode). The default maximum value is 255. - To change the color of an image (or a texture), use (tint). + To change the color of an image (or a texture), use (tint). @@ -633,7 +633,7 @@ Returns: None - + References: http://processing.org/reference/fill_.html @@ -646,27 +646,27 @@ (fill sktch 255 0 255) " - ([#^PApplet sketch gray] - (cond + ([#^PApplet sketch gray] + (cond (string? gray) (.fill sketch (color sketch gray)) (integer? gray) (.fill sketch (int gray)) (float? gray) (.fill sketch (float gray)))) - ([#^PApplet sketch gray alpha] - (cond + ([#^PApplet sketch gray alpha] + (cond (string? gray) (.fill sketch (color sketch gray) (int alpha)) (integer? gray) (.fill sketch (int gray) (int alpha)) (float? gray) (.fill sketch (float gray) (float alpha)))) - ([#^PApplet sketch x y z] + ([#^PApplet sketch x y z] (if (or (float? x) (float? y) (float? z)) (.fill sketch (float x) (float y) (float z)) (.fill sketch (int x) (int y) (int z)))) - ([#^PApplet sketch x y z alpha] + ([#^PApplet sketch x y z alpha] (if (or (float? x) (float? y) (float? z)) (.fill sketch (float x) (float y) (float z) (float alpha)) (.fill sketch (int x) (int y) (int z) (int alpha))))) @@ -685,7 +685,7 @@ (defn framerate [#^PApplet sketch new-rate] (.frameRate sketch (float new-rate))) -;; added by DEL +;; added by DEL (defn frame-count [#^PApplet sketch] (.frameCount sketch)) @@ -728,11 +728,11 @@ ;; added by DEL -(defn lerp-color +(defn lerp-color " - Calculates a color or colors between two color at a specific increment. - The amt parameter is the amount to interpolate between the two values - where 0.0 equal to the first point, 0.1 is very near the first point, + Calculates a color or colors between two color at a specific increment. + The amt parameter is the amount to interpolate between the two values + where 0.0 equal to the first point, 0.1 is very near the first point, 0.5 is half-way in between, etc. Parameters: @@ -742,9 +742,9 @@ Returns: float - + " - ([c1 c2 amt] + ([c1 c2 amt] (PApplet/lerpColor c1 c2 (float amt) 1))) ;; use RGB mode @@ -772,14 +772,14 @@ (defn load-pixels [#^PApplet sketch ] (.loadPixels sketch)) (defn load-shape - "Load a geometry from a file as a PShape." - [#^PApplet sketch filename] - (.loadShape sketch filename)) + "Load a geometry from a file as a PShape." + [#^PApplet sketch filename] + (.loadShape sketch filename)) -(defn load-strings - "Load data from a file and shove it into a String array." - [#^PApplet sketch filename] - (.loadStrings sketch filename)) +(defn load-strings + "Load data from a file and shove it into a String array." + [#^PApplet sketch filename] + (.loadStrings sketch filename)) ;; $$log @@ -848,9 +848,9 @@ (defn no-loop [#^PApplet sketch ] (.noLoop sketch)) (defn norm - "Normalize a value to exist between 0 and 1 (inclusive)." - [val start stop] - (PApplet/norm (float val) (float start) (float stop))) + "Normalize a value to exist between 0 and 1 (inclusive)." + [val start stop] + (PApplet/norm (float val) (float start) (float stop))) (defn normal [#^PApplet sketch nx ny nz] (.normal sketch (float nx) (float ny) (float nz))) @@ -864,7 +864,7 @@ ;; $$open -- overload -(defn ortho +(defn ortho ([#^PApplet sketch ] (.ortho sketch)) ([#^PApplet sketch l r b t near far] (.ortho sketch (float l) (float r) (float b) (float t) (float near) (float far)))) @@ -879,7 +879,7 @@ (defn perspective ([#^PApplet sketch ] (.perspective sketch)) ([#^PApplet sketch fovy aspect z-near z-far] - (.perspective sketch (float fovy) (float aspect) (float z-near) (float z-far)))) + (.perspective sketch (float fovy) (float aspect) (float z-near) (float z-far)))) ;; pmouse-x doesn't work, always returns 0 ;;(defn pmouse-x [#^PApplet sketch ] (.pmouseX sketch)) @@ -943,8 +943,8 @@ ;; $$registerSize (defn request-image - ([#^PApplet sketch filename] (.requestImage sketch filename)) - ([#^PApplet sketch filename extension] (.requestImage sketch filename extension))) + ([#^PApplet sketch filename] (.requestImage sketch filename)) + ([#^PApplet sketch filename extension] (.requestImage sketch filename extension))) (defn reset-matrix [#^PApplet sketch ] (.resetMatrix sketch)) @@ -967,9 +967,9 @@ (defn saturation [#^PApplet sketch what] (.saturation sketch (int what))) ;; modified by DEL -(defmethod save :sketch ([#^PApplet sketch filename] - (doto sketch - ;.redraw +(defmethod save :sketch ([#^PApplet sketch filename] + (doto sketch + ;.redraw (.save filename)))) ;; $$saveBytes @@ -1051,9 +1051,9 @@ (defn spotlight ([#^PApplet sketch r g b x y z nx ny nz angle concentration] - (.spotLight sketch r g b x y z nx ny nz angle concentration)) + (.spotLight sketch r g b x y z nx ny nz angle concentration)) ([#^PApplet sketch [r g b] [x y z] [nx ny nz] angle concentration] - (.spotLight sketch r g b x y z nx ny nz angle concentration))) + (.spotLight sketch r g b x y z nx ny nz angle concentration))) ;;(defn sq [a] (PApplet/sq (float a))) @@ -1066,25 +1066,25 @@ ;; modified by DEL (defn stroke - ([#^PApplet sketch gray] - (cond + ([#^PApplet sketch gray] + (cond (string? gray) (.stroke sketch (color sketch gray)) (integer? gray) (.stroke sketch (int gray)) (float? gray) (.stroke sketch (float gray)))) - ([#^PApplet sketch gray alpha] - (cond + ([#^PApplet sketch gray alpha] + (cond (string? gray) (.stroke sketch (color sketch gray) (int alpha)) (integer? gray) (.stroke sketch (int gray) (int alpha)) (float? gray) (.stroke sketch (float gray) (float alpha)))) - ([#^PApplet sketch x y z] + ([#^PApplet sketch x y z] (.stroke sketch (float x) (float y) (float z))) - ([#^PApplet sketch x y z alpha] + ([#^PApplet sketch x y z alpha] (.stroke sketch (float x) (float y) (float z) (float alpha)))) @@ -1100,16 +1100,16 @@ ;; added by DEL (defn text - ([#^PApplet sketch s] + ([#^PApplet sketch s] (.text sketch s)) - ([#^PApplet sketch s x y] + ([#^PApplet sketch s x y] (.text sketch s (float x) (float y))) - ([#^PApplet sketch s x y z] + ([#^PApplet sketch s x y z] (.text sketch s (float x) (float y) (float z))) ([#^PApplet sketch s x1 y1 x2 y2] - (.text sketch s (float x1) (float y1) (float x2) (float y2))) + (.text sketch s (float x1) (float y1) (float x2) (float y2))) ([#^PApplet sketch s x1 y1 x2 y2 z] - (.text sketch s (float x1) (float y1) (float x2) (float y2) (float z)))) + (.text sketch s (float x1) (float y1) (float x2) (float y2) (float z)))) (defn char->text @@ -1128,9 +1128,9 @@ (defn string->text-in ([#^PApplet sketch #^java.lang.String s x1 y1 x2 y2] - (.text sketch s (float x1) (float y1) (float x2) (float y2))) + (.text sketch s (float x1) (float y1) (float x2) (float y2))) ([#^PApplet sketch #^java.lang.String s x1 y1 x2 y2 z] - (.text sketch s (float x1) (float y1) (float x2) (float y2) (float z)))) + (.text sketch s (float x1) (float y1) (float x2) (float y2) (float z)))) (defn text-align ([#^PApplet sketch align] (.textAlign sketch (int align))) @@ -1162,36 +1162,36 @@ (defmethod text-width false [#^PApplet sketch #^java.lang.String s] (.textWidth sketch s)) -;; modified by DEL +;; modified by DEL (defn tint - ([#^PApplet sketch gray] - (cond + ([#^PApplet sketch gray] + (cond (string? gray) (.tint sketch (color sketch gray)) (integer? gray) (.tint sketch (int gray)) (float? gray) (.tint sketch (float gray)))) - ([#^PApplet sketch gray alpha] - (cond + ([#^PApplet sketch gray alpha] + (cond (string? gray) (.tint sketch (color sketch gray) (int alpha)) (integer? gray) (.tint sketch (int gray) (int alpha)) (float? gray) (.tint sketch (float gray) (float alpha)))) - ([#^PApplet sketch x y z] + ([#^PApplet sketch x y z] (if (or (integer? x) (integer? y) (integer? z)) (.tint sketch (int x)(int y) (int z)) (.tint sketch (float x)(float y) (float z)))) - ([#^PApplet sketch x y z a] + ([#^PApplet sketch x y z a] (if (or (integer? x) (integer? y) (integer? z)) (.tint sketch (int x)(int y) (int z) (int alpha)) (.tint sketch (float x)(float y) (float z) (float alpha))))) (defn translate - ([v] (apply translate v)) + ([v] (apply translate v)) ([#^PApplet sketch tx ty] (.translate sketch (float tx) (float ty))) ([#^PApplet sketch tx ty tz] (.translate sketch (float tx) (float ty) (float tz)))) @@ -1219,44 +1219,44 @@ ([#^PApplet sketch x y z] (.vertex sketch (float x) (float y) (float z))) ([#^PApplet sketch x y u v] (.vertex sketch (float x) (float y) (float u) (float v))) ([#^PApplet sketch x y z u v] - (.vertex sketch (float x) (float y) (float z) (float u) (float v)))) + (.vertex sketch (float x) (float y) (float z) (float u) (float v)))) ;; added by DEL (defn width [sketch] (.width sketch)) (defn height [sketch] (.height sketch)) - - + + (defn year [] (PApplet/year)) ;; utility functions. clj-processing specific (defmacro with-translation - "Berforms body with translation, restores current transformation on exit." - [translation-vector & body] - `(let [tr# ~translation-vector] - (push-matrix) - (translate tr#) - ~@body - (pop-matrix))) - -(defmacro with-rotation - "Berforms body with rotation, restores current transformation on exit. + "Berforms body with translation, restores current transformation on exit." + [translation-vector & body] + `(let [tr# ~translation-vector] + (push-matrix) + (translate tr#) + ~@body + (pop-matrix))) + +(defmacro with-rotation + "Berforms body with rotation, restores current transformation on exit. Accepts a vector [angle] or [angle x-axis y-axis z-axis]. - Example: - (with-rotation [angle] + Example: + (with-rotation [angle] (vertex 1 2))" - [rotation & body] - `(let [tr# ~rotation] - (push-matrix) - (apply rotate tr#) - ~@body - (pop-matrix))) + [rotation & body] + `(let [tr# ~rotation] + (push-matrix) + (apply rotate tr#) + ~@body + (pop-matrix))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; These CONSTANTS are located in the constants.clj file in the original +;; These CONSTANTS are located in the constants.clj file in the original ;; clj-processing library ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; @@ -1418,7 +1418,7 @@ ;; (def POINT PApplet/POINT) ; shared with shape feature (def SPOT PApplet/SPOT) -;; keys +;; keys (def BACKSPACE PApplet/BACKSPACE) (def TAB PApplet/TAB) @@ -1455,19 +1455,19 @@ Examples: -(use '(incanter core processing)) +(use '(incanter core processing)) (let [map-image (ref nil) sktch (sketch - + ;; define the setup function (setup [] - (dosync (ref-set map-image + (dosync (ref-set map-image (load-image this \"examples/images/map.png\"))) (size this 640 400)) - - ;; define the draw function + + ;; define the draw function (draw [] (doto this ;(background-float 255) @@ -1493,8 +1493,8 @@ title (or (:title opts) "Processing Sketch") width (or (:width opts) (.width (.getSize sketch))) height (or (:height opts) (.height (.getSize sketch))) - [width height] (or (:size opts) - [(.width (.getSize sketch)) + [width height] (or (:size opts) + [(.width (.getSize sketch)) (.height (.getSize sketch))]) frame (javax.swing.JFrame. title) ] @@ -1503,5 +1503,3 @@ (.setDefaultCloseOperation javax.swing.WindowConstants/EXIT_ON_CLOSE) (.setSize width height) (.setVisible true))))) - - src/incanter/processing.clj @@ -1,4 +1,3 @@ - ;;; som.clj -- Self-Organizing-Map Neural Network Library ;; by David Edgar Liebke http://incanter.org @@ -14,13 +13,13 @@ -(ns incanter.som +(ns incanter.som (:use [incanter.core :only (sel ncol nrow mult div plus minus trans to-vect sqrt sum pow)] [incanter.stats :only (mean principal-components covariance)])) -(defn- som-dimensions +(defn- som-dimensions ([pc1-sd pc2-sd] (let [dim-1 (mult 5 pc1-sd) dim-2 (mult (div pc2-sd pc1-sd) dim-1)] @@ -36,17 +35,17 @@ pc2 (sel (:rotation pc) :cols 1) [dim-1 dim-2] (map #(Math/ceil %) (som-dimensions pc1-sd pc2-sd)) data-mean (map mean (trans data)) - weight-fn (fn [i j data-mean pc1-sd dim-1 dim-2 pc1 pc2] + weight-fn (fn [i j data-mean pc1-sd dim-1 dim-2 pc1 pc2] (to-vect (plus data-mean (mult (mult 5 pc1-sd) (plus (mult pc1 (minus i (div dim-1 2))) (mult pc2 (minus j (div dim-2 2)))))))) weights (reduce conj {} - (for [i (range dim-1) j (range dim-2)] + (for [i (range dim-1) j (range dim-2)] {[i j] (weight-fn i j data-mean pc1-sd dim-1 dim-2 pc1 pc2)} )) sets (reduce conj {} - (for [i (range dim-1) j (range dim-2)] + (for [i (range dim-1) j (range dim-2)] {[i j] #{}} ))] {:dims [dim-1 dim-2] :weights weights @@ -57,37 +56,37 @@ (defn- dist-euclidean [x y] (sqrt (sum (pow (minus x y) 2)))) -(defn- get-distances - ([x som] +(defn- get-distances + ([x som] (reduce conj {} - (pmap #(hash-map % (dist-euclidean x ((:weights som) %))) + (pmap #(hash-map % (dist-euclidean x ((:weights som) %))) (keys (:weights som)))))) (defn- get-min-dist ([x som] (let [distances (get-distances x som) - min-dist-key (reduce #(if (<= (distances %1) (distances %2)) %1 %2) + min-dist-key (reduce #(if (<= (distances %1) (distances %2)) %1 %2) (keys distances))] {:index min-dist-key :dist (distances min-dist-key)} ))) (defn- som-update-cells ([data som] - (let [sets (loop [i 0 sets {}] + (let [sets (loop [i 0 sets {}] (if (= i (nrow data)) sets (let [{idx :index min-dist :dist} (get-min-dist (trans (nth data i)) som)] (recur (inc i) (assoc sets idx (conj (sets idx) i))))))] (assoc som :sets sets)))) - -(defn- alpha-fn + +(defn- alpha-fn ([r total-cycles alpha-init] (max 0.01 (* alpha-init (- 1 (/ r total-cycles)))))) -(defn- beta-fn +(defn- beta-fn ([r beta-init] (max 0 (- beta-init r)))) @@ -95,11 +94,11 @@ (defn- som-neighborhoods ([r dim-1 dim-2 total-cycles beta0] (reduce conj {} - (for [i (range dim-1) j (range dim-2)] - [[i j] - (for [s1 (range (if (pos? (- i (beta-fn r beta0))) (- i (beta-fn r beta0)) 0) + (for [i (range dim-1) j (range dim-2)] + [[i j] + (for [s1 (range (if (pos? (- i (beta-fn r beta0))) (- i (beta-fn r beta0)) 0) (if (<= (+ i (beta-fn r beta0)) dim-1) (+ i (beta-fn r beta0) 1) dim-1)) - s2 (range (if (pos? (- j (beta-fn r beta0))) (- j (beta-fn r beta0)) 0) + s2 (range (if (pos? (- j (beta-fn r beta0))) (- j (beta-fn r beta0)) 0) (if (<= (+ j (beta-fn r beta0)) dim-2) (+ j (beta-fn r beta0) 1) dim-2))] [s1 s2])])))) @@ -107,8 +106,8 @@ (defn- som-update-weights [r total-cycles som alpha-init beta-init] (let [ - sets (:sets som) - weights (:weights som) + sets (:sets som) + weights (:weights som) indices (keys weights) dims (:dims som) neighborhoods (som-neighborhoods r (first dims) (second dims) total-cycles beta-init) @@ -117,7 +116,7 @@ (reduce conj {} (pmap (fn [indx] {indx - (plus (weights indx) + (plus (weights indx) (mult (alpha-fn r total-cycles alpha-init) (minus (if (pos? (count (sets indx))) (div (sum (sets indx)) @@ -126,17 +125,17 @@ (weights indx))))} ) indices))))) -(defn- som-fitness +(defn- som-fitness ([data som] (/ (sum (for [indx (keys (:weights som))] - (sum (map #(dist-euclidean ((:weights som) indx) (trans (nth data %))) + (sum (map #(dist-euclidean ((:weights som) indx) (trans (nth data %))) ((:sets som) indx))))) (nrow data)))) -(defn som-batch-train -" Performs BL-SOM (batch-learning self organizing map) learning on - the given data, returning a hashmap containing resulting BL-SOM +(defn som-batch-train +" Performs BL-SOM (batch-learning self organizing map) learning on + the given data, returning a hashmap containing resulting BL-SOM values. @@ -162,7 +161,7 @@ Examples: (use '(incanter core som stats charts datasets)) - (def data (to-matrix (sel (get-dataset :iris) + (def data (to-matrix (sel (get-dataset :iris) :cols [\"Sepal.Length\" \"Sepal.Width\" \"Petal.Length\" \"Petal.Width\"]))) (def som (som-batch-train data :cycles 10 :alpha 0.5 :beta 3)) @@ -172,8 +171,8 @@ ;; view indices of data items in each cell (:sets som) ;; view the species in each cell - (doseq [rws (vals (:sets som))] - (println (sel (get-dataset :iris) :cols \"Species\" :rows rws) \\newline)) + (doseq [rws (vals (:sets som))] + (println (sel (get-dataset :iris) :cols \"Species\" :rows rws) \\newline)) ;; plot the means of the data vectors in each cell/cluster (def cell-means (map #(map mean (trans (sel data :rows ((:sets som) %)))) (keys (:sets som)))) @@ -200,9 +199,3 @@ (let [new-som (som-update-weights r total-cycles (som-update-cells data som) alpha-init beta-init)] (recur (inc r) new-som (conj fit (som-fitness data new-som))))))))) - - - - - - src/incanter/som.clj @@ -16,7 +16,7 @@ -(ns incanter.stats +(ns incanter.stats (:import (cern.colt.list.tdouble DoubleArrayList) (cern.jet.random.tdouble Gamma Beta Binomial ChiSquare DoubleUniform Exponential NegativeBinomial Normal Poisson @@ -24,7 +24,7 @@ (cern.jet.random.tdouble.engine DoubleMersenneTwister) (cern.jet.stat.tdouble DoubleDescriptive Probability)) - (:use [incanter.core :only (abs plus minus div mult mmult to-list bind-columns + (:use [incanter.core :only (abs plus minus div mult mmult to-list bind-columns gamma pow sqrt diag trans regularized-beta ncol nrow identity-matrix decomp-cholesky decomp-svd matrix length sum sum-of-squares sel matrix? @@ -33,7 +33,7 @@ -(defn indicator +(defn indicator " Returns a sequence of ones and zeros, where ones are returned when the given predicate is true for @@ -42,7 +42,7 @@ Examples: (use 'incanter.stats) - + (indicator #(neg? %) (sample-normal 10)) ;; return the sum of the positive values in a normal sample @@ -50,40 +50,40 @@ (sum (mult x (indicator #(pos? %) x))) " - ([pred coll] + ([pred coll] (for [el coll] (if (pred el) 1 0)))) -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; CONTINOUS DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; F DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + - (defn pdf-f -" Returns the F pdf of the given value, x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's df function. +" Returns the F pdf of the given value, x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's df function. - Options: - :df1 (default 1) + Options: + :df1 (default 1) :df2 (default 1) - See also: + See also: cdf-f and quantile-f - References: + References: http://en.wikipedia.org/wiki/F_distribution http://mathworld.wolfram.com/F-Distribution.html http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-f 1.0 :df1 5 :df2 2) " ([x & options] @@ -91,7 +91,7 @@ df1 (or (:df1 opts) 0) df2 (or (:df2 opts) 1) pdf-fx (fn [x] - (* (/ (gamma (/ (+ df1 df2) 2)) + (* (/ (gamma (/ (+ df1 df2) 2)) (* (gamma (/ df1 2)) (gamma (/ df2 2)))) (pow (/ df1 df2) (/ df1 2)) (pow x (- (/ df1 2) 1)) @@ -104,23 +104,23 @@ -(defn cdf-f -" Returns the F-distribution cdf of the given value, x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's pf function. +(defn cdf-f +" Returns the F-distribution cdf of the given value, x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's pf function. - Options: - :df1 (default 1) + Options: + :df1 (default 1) :df2 (default 1) - See also: + See also: pdf-f and quantile-f - References: + References: http://en.wikipedia.org/wiki/F_distribution http://mathworld.wolfram.com/F-Distribution.html http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-f 1.0 :df1 5 :df2 2) " ([x & options] @@ -129,11 +129,11 @@ df1 (or (:df1 opts) 1) df2 (or (:df2 opts) 1) cdf-fx (if lower-tail? - (fn [x1] (regularized-beta + (fn [x1] (regularized-beta (/ (* df1 x1) (+ df2 (* df1 x1))) (/ df1 2) (/ df2 2))) - (fn [x1] (- 1 (regularized-beta + (fn [x1] (- 1 (regularized-beta (/ (* df1 x1) (+ df2 (* df1 x1))) (/ df1 2) (/ df2 2))))) @@ -144,27 +144,27 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; NORMAL DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-normal -" Returns the Normal pdf of the given value, x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's dnorm function. +" Returns the Normal pdf of the given value, x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's dnorm function. - Options: - :mean (default 0) + Options: + :mean (default 0) :sd (default 1) - See also: + See also: cdf-normal, quantile-normal, sample-normal - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Normal.html http://en.wikipedia.org/wiki/Normal_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-normal 1.96 :mean -2 :sd (sqrt 0.5)) " ([x & options] @@ -179,55 +179,55 @@ (defn cdf-normal -" Returns the Normal cdf of the given value, x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's pnorm function. +" Returns the Normal cdf of the given value, x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's pnorm function. - Options: - :mean (default 0) + Options: + :mean (default 0) :sd (default 1) - See also: + See also: pdf-normal, quantile-normal, sample-normal - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Normal.html http://en.wikipedia.org/wiki/Normal_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-normal 1.96 :mean -2 :sd (sqrt 0.5)) " ([x & options] (let [opts (when options (apply assoc {} options)) mean (or (:mean opts) 0) sd (or (:sd opts) 1) - dist (Normal. mean sd (DoubleMersenneTwister.))] + dist (Normal. mean sd (DoubleMersenneTwister.))] (if (coll? x) (map #(.cdf dist %) x) (.cdf dist x))))) -(defn quantile-normal -" Returns the inverse of the Normal CDF for the given probability. - It will return a sequence of values, if given a sequence of +(defn quantile-normal +" Returns the inverse of the Normal CDF for the given probability. + It will return a sequence of values, if given a sequence of probabilities. This is equivalent to R's qnorm function. - Options: - :mean (default 0) + Options: + :mean (default 0) :sd (default 1) - Returns: + Returns: a value x, where (cdf-normal x) = probability - See also: + See also: pdf-normal, cdf-normal, and sample-normal - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/stat/tdouble/Probability.html http://en.wikipedia.org/wiki/Normal_distribution http://en.wikipedia.org/wiki/Quantile - Example: + Example: (quantile-normal 0.975) (quantile-normal [0.025 0.975] :mean -2 :sd (sqrt 0.5)) " @@ -235,7 +235,7 @@ (let [opts (when options (apply assoc {} options)) mean (or (:mean opts) 0) sd (or (:sd opts) 1) - x (if (coll? probability) + x (if (coll? probability) (map #(Probability/normalInverse %) probability) (Probability/normalInverse probability))] (plus mean (mult sd x))))) @@ -244,20 +244,20 @@ (defn sample-normal " Returns a sample of the given size from a Normal distribution - This is equivalent to R's rnorm function. + This is equivalent to R's rnorm function. - Options: - :mean (default 0) + Options: + :mean (default 0) :sd (default 1) - See also: + See also: pdf-normal, cdf-normal, quantile-normal - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Normal.html http://en.wikipedia.org/wiki/Normal_distribution - Example: + Example: (sample-normal 1000 :mean -2 :sd (sqrt 0.5)) " ([size & options] @@ -274,15 +274,15 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn sample-mvn -" Returns a sample of the given size from a Multivariate Normal - distribution. This is equivalent to R's mvtnorm::rmvnorm function. +(defn sample-mvn +" Returns a sample of the given size from a Multivariate Normal + distribution. This is equivalent to R's mvtnorm::rmvnorm function. Arguments: size -- the size of the sample to return - Options: - :mean (default (repeat (ncol sigma) 0)) + Options: + :mean (default (repeat (ncol sigma) 0)) :sigma (default (identity-matrix (count mean))) @@ -312,17 +312,17 @@ " ([size & options] (let [opts (when options (apply assoc {} options)) - mean (or (:mean opts) - (if (:sigma opts) + mean (or (:mean opts) + (if (:sigma opts) (repeat (ncol (:sigma opts)) 0) [0])) - sigma (or (:sigma opts) + sigma (or (:sigma opts) (identity-matrix (count mean))) p (count mean) chol (decomp-cholesky sigma) norm-samp (mmult (matrix (sample-normal (* size p)) p) chol) ] - (if (> (nrow norm-samp) 1) + (if (> (nrow norm-samp) 1) (matrix (map #(plus % (trans mean)) norm-samp)) (matrix (plus norm-samp (trans mean))))))) @@ -334,23 +334,23 @@ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn pdf-uniform -" Returns the Uniform pdf of the given value of x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's dunif function. +(defn pdf-uniform +" Returns the Uniform pdf of the given value of x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's dunif function. - Options: - :min (default 0) + Options: + :min (default 0) :max (default 1) - See also: + See also: cdf-uniform and sample-uniform - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/DoubleUniform.html http://en.wikipedia.org/wiki/Uniform_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-uniform 5) (pdf-uniform 5 :min 1 :max 10) " @@ -364,23 +364,23 @@ (.pdf dist x))))) -(defn cdf-uniform -" Returns the Uniform cdf of the given value of x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's punif function. +(defn cdf-uniform +" Returns the Uniform cdf of the given value of x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's punif function. - Options: - :min (default 0) + Options: + :min (default 0) :max (default 1) - See also: + See also: pdf-uniform and sample-uniform - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/DoubleUniform.html http://en.wikipedia.org/wiki/Uniform_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-uniform 5) (cdf-uniform 5 :min 1 :max 10) " @@ -394,23 +394,23 @@ (.cdf dist x))))) -(defn sample-uniform +(defn sample-uniform " Returns a sample of the given size from a Uniform distribution. - This is equivalent to R's runif function. + This is equivalent to R's runif function. - Options: - :min (default 0) + Options: + :min (default 0) :max (default 1) :integers (default false) - See also: + See also: pdf-uniform and cdf-uniform - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/DoubleUniform.html http://en.wikipedia.org/wiki/Uniform_distribution - Example: + Example: (sample-uniform 1000) (sample-uniform 1000 :min 1 :max 10) " @@ -421,35 +421,35 @@ ints? (if (true? (:integers opts)) true false) dist (DoubleUniform. min-val max-val (DoubleMersenneTwister.))] (if (= size 1) - (if ints? + (if ints? (DoubleUniform/staticNextIntFromTo min-val max-val) (DoubleUniform/staticNextDoubleFromTo min-val max-val)) - (if ints? + (if ints? (for [_ (range size)] (DoubleUniform/staticNextIntFromTo min-val max-val)) (for [_ (range size)] (DoubleUniform/staticNextDoubleFromTo min-val max-val))))))) -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; BETA DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-beta -" Returns the Beta pdf of the given value of x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's dbeta function. +" Returns the Beta pdf of the given value of x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's dbeta function. - Options: - :alpha (default 1) + Options: + :alpha (default 1) :beta (default 1) - See also: + See also: cdf-beta and sample-beta - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Beta.html http://en.wikipedia.org/wiki/Beta_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-beta 0.5 :alpha 1 :beta 2) " ([x & options] @@ -463,23 +463,23 @@ (defn cdf-beta -" Returns the Beta cdf of the given value of x. It will return a sequence - of values, if x is a sequence. This is equivalent to R's pbeta function. +" Returns the Beta cdf of the given value of x. It will return a sequence + of values, if x is a sequence. This is equivalent to R's pbeta function. - Options: - :alpha (default 1) + Options: + :alpha (default 1) :beta (default 1) :lower-tail (default true) - See also: + See also: pdf-beta and sample-beta - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Beta.html http://en.wikipedia.org/wiki/Beta_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-beta 0.5 :alpha 1 :beta 2) (cdf-beta 0.5 :alpha 1 :beta 2 :lower-tail false) " @@ -503,19 +503,19 @@ " Returns a sample of the given size from a Beta distribution. This is equivalent to R's rbeta function. - Options: - :alpha (default 1) + Options: + :alpha (default 1) :beta (default 1) These default values produce a Uniform distribution. - See also: + See also: pdf-beta and cdf-beta - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Beta.html http://en.wikipedia.org/wiki/Beta_distribution - Example: + Example: (sample-beta 1000 :alpha 1 :beta 2) " ([size & options] @@ -527,31 +527,31 @@ (for [_ (range size)] (Beta/staticNextDouble alpha beta)))))) -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; GAMMA DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-gamma -" Returns the Gamma pdf for the given value of x. It will return a sequence +" Returns the Gamma pdf for the given value of x. It will return a sequence of values, if x is a sequence. This is equivalent to R's dgamma function. - Options: - :shape (default 1) + Options: + :shape (default 1) :rate (default 1) - See also: + See also: cdf-gamma and sample-gamma - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Gamma.html http://en.wikipedia.org/wiki/Gamma_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-gamma 10 :shape 1 :rate 2) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) shape (or (:shape opts) 1) rate (or (:rate opts) 1) dist (Gamma. shape rate (DoubleMersenneTwister.))] @@ -562,31 +562,31 @@ (defn cdf-gamma -" Returns the Gamma cdf for the given value of x. It will return a sequence +" Returns the Gamma cdf for the given value of x. It will return a sequence of values, if x is a sequence. This is equivalent to R's pgamma function. - Options: - :shape (default 1) + Options: + :shape (default 1) :rate (default 1) :lower-tail (default true) - See also: + See also: pdf-gamma and sample-gamma - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Gamma.html http://en.wikipedia.org/wiki/Gamma_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-gamma 10 :shape 1 :rate 2) (cdf-gamma 3 :shape 1 :lower-tail false) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) shape (or (:shape opts) 1) rate (or (:rate opts) 1) - lower-tail? (if (false? (:lower-tail opts)) false true) + lower-tail? (if (false? (:lower-tail opts)) false true) cdf-fx (if lower-tail? (fn [x1] (Probability/gamma rate shape x1)) (fn [x1] (Probability/gammaComplemented rate shape x1)))] @@ -600,22 +600,22 @@ " Returns a sample of the given size from a Gamma distribution. This is equivalent to R's rgamma function. - Options: - :shape (default 1) + Options: + :shape (default 1) :rate (default 1) - See also: + See also: pdf-gamma, cdf-gamma, and quantile-gamma - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Gamma.html http://en.wikipedia.org/wiki/Gamma_distribution - Example: + Example: (sample-gamma 1000 :shape 1 :rate 2) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) shape (or (:shape opts) 1) rate (or (:rate opts) 1)] (if (= size 1) @@ -625,30 +625,30 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; CHI SQUARE DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-chisq -" Returns the Chi Square pdf of the given value of x. It will return a sequence +" Returns the Chi Square pdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dchisq function. - Options: - :df (default 1) + Options: + :df (default 1) - See also: + See also: cdf-chisq and sample-chisq - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/ChiSquare.html http://en.wikipedia.org/wiki/Chi_square_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-chisq 5.0 :df 2) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1) dist (ChiSquare. df (DoubleMersenneTwister.))] (if (coll? x) @@ -658,27 +658,27 @@ (defn cdf-chisq -" Returns the Chi Square cdf of the given value of x. It will return a sequence +" Returns the Chi Square cdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's pchisq function. - Options: - :df (default 1) - :lower-tail (default true) + Options: + :df (default 1) + :lower-tail (default true) - See also: + See also: pdf-chisq and sample-chisq - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/ChiSquare.html http://en.wikipedia.org/wiki/Chi_square_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-chisq 5.0 :df 2) (cdf-chisq 5.0 :df 2 :lower-tail false) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1) lower-tail? (if (false? (:lower-tail opts)) false true) cdf-fx (if lower-tail? @@ -694,21 +694,21 @@ " Returns a sample of the given size from a Chi Square distribution Same as R's rchisq function. - Options: - :df (default 1) + Options: + :df (default 1) - See also: + See also: pdf-chisq and cdf-chisq - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/ChiSquare.html http://en.wikipedia.org/wiki/Chi_square_distribution - Example: + Example: (sample-chisq 1000 :df 2) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1)] (if (= size 1) (ChiSquare/staticNextDouble df) @@ -716,30 +716,30 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; STUDENT'S T DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-t -" Returns the Student's t pdf for the given value of x. It will return a sequence +" Returns the Student's t pdf for the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dt function. - Options: - :df (default 1) + Options: + :df (default 1) - See also: + See also: cdf-t, quantile-t, and sample-t - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/StudentT.html http://en.wikipedia.org/wiki/Student-t_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-t 1.2 :df 10) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1) dist (StudentT. df (DoubleMersenneTwister.))] (if (coll? x) @@ -748,25 +748,25 @@ (defn cdf-t -" Returns the Student's t cdf for the given value of x. It will return a sequence +" Returns the Student's t cdf for the given value of x. It will return a sequence of values, if x is a sequence. Same as R's pt function. - Options: - :df (default 1) + Options: + :df (default 1) - See also: + See also: pdf-t, quantile-t, and sample-t - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/StudentT.html http://en.wikipedia.org/wiki/Student-t_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-t 1.2 :df 10) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1) lower-tail? (if (false? (:lower-tail opts)) false true) cdf-fx (if lower-tail? @@ -777,26 +777,26 @@ (cdf-fx x))))) -(defn quantile-t +(defn quantile-t " Returns the inverse of the Student's t CDF for the given probability - (i.e. the quantile). It will return a sequence of values, if x is + (i.e. the quantile). It will return a sequence of values, if x is a sequence of probabilities. This is equivalent to R's qt function. - Options: - :df (default 1) + Options: + :df (default 1) - Returns: + Returns: a value x, where (cdf-t x) = probability - See also: - pdf-t, cdf-t, and sample-t + See also: + pdf-t, cdf-t, and sample-t - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/stat/tdouble/Probability.html http://en.wikipedia.org/wiki/Student-t_distribution http://en.wikipedia.org/wiki/Quantile - Example: + Example: (quantile-t 0.975) (quantile-t [0.025 0.975] :df 25) (def df [1 2 3 4 5 6 7 8 9 10 20 50 100 1000]) @@ -806,11 +806,11 @@ (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1) to-alpha (fn [prob] ;; need to convert the probability to an alpha value - (if (< prob 1/2) + (if (< prob 1/2) (* 2 prob) (* 2 (- 1 prob)))) sign-fx (fn [x1 prob] (if (< prob 1/2) (* -1 x1) x1)) - x (if (coll? probability) + x (if (coll? probability) (map #(sign-fx (Probability/studentTInverse (to-alpha %) df) %) probability) (sign-fx (Probability/studentTInverse (to-alpha probability) df) probability))] x))) @@ -821,21 +821,21 @@ " Returns a sample of the given size from a Student's t distribution. Same as R's rt function. - Options: - :df (default 1) + Options: + :df (default 1) - See also: + See also: pdf-t, cdf-t, and quantile-t - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/StudentT.html http://en.wikipedia.org/wiki/Student-t_distribution - Example: + Example: (cdf-t 1000 :df 10) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) df (or (:df opts) 1)] (if (= size 1) (StudentT/staticNextDouble df) @@ -845,30 +845,30 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; EXPONENTIAL DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-exp -" Returns the Exponential pdf of the given value of x. It will return a sequence +" Returns the Exponential pdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dexp - Options: - :rate (default 1) + Options: + :rate (default 1) - See also: + See also: cdf-exp and sample-exp - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Exponential.html http://en.wikipedia.org/wiki/Exponential_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-exp 2.0 :rate 1/2) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) rate (or (:rate opts) 1) dist (Exponential. rate (DoubleMersenneTwister.))] (if (coll? x) @@ -877,25 +877,25 @@ (defn cdf-exp -" Returns the Exponential cdf of the given value of x. It will return a sequence +" Returns the Exponential cdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's pexp - Options: - :rate (default 1) + Options: + :rate (default 1) - See also: + See also: pdf-exp and sample-exp - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Exponential.html http://en.wikipedia.org/wiki/Exponential_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-exp 2.0 :rate 1/2) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) rate (or (:rate opts) 1) dist (Exponential. rate (DoubleMersenneTwister.))] (if (coll? x) @@ -907,21 +907,21 @@ " Returns a sample of the given size from a Exponential distribution. Same as R's rexp - Options: - :rate (default 1) + Options: + :rate (default 1) - See also: + See also: pdf-exp, and cdf-exp - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Exponential.html http://en.wikipedia.org/wiki/Exponential_distribution - Example: + Example: (sample-exp 1000 :rate 1/2) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) rate (or (:rate opts) 1)] (if (= size 1) (Exponential/staticNextDouble rate) @@ -930,9 +930,9 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; WISHART DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn sample-wishart @@ -961,7 +961,7 @@ scale (or (:scale opts) (when (:p opts) (identity-matrix (:p opts)))) p (count scale) df (or (:df opts) p) - diagonal (for [i (range 1 (inc p))] + diagonal (for [i (range 1 (inc p))] (pow (sample-chisq 1 :df (inc (- df i))) 1/2)) mat (diag diagonal) indices (for [i (range p) j (range p) :when (< j i)] [i j]) @@ -1015,7 +1015,7 @@ ;; the counts y1, y2, and y3 are assumed to have a multinomial distribution ;; If a uniform prior distribution is assigned to the multinomial vector ;; theta = (th1, th2, th3), then the posterior distribution of theta is - ;; proportional to g(theta) = th1^y1 * th2^y2 * th3^y3, which is a + ;; proportional to g(theta) = th1^y1 * th2^y2 * th3^y3, which is a ;; dirichlet distribution with parameters (y1+1, y2+1, y3+1) (def theta (sample-dirichlet 1000 [(inc 727) (inc 583) (inc 137)])) ;; view means, 95% CI, and histograms of the proportion parameters @@ -1041,40 +1041,40 @@ (matrix (map #(div %1 %2) W T))))) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; DISCRETE DISTRIBUTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; BINOMIAL DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-binomial -" Returns the Bionomial pdf of the given value of x. It will return a sequence +" Returns the Bionomial pdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dbinom - Options: - :size (default 1) + Options: + :size (default 1) :prob (default 1/2) - See also: + See also: cdf-binomial and sample-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Binomial.html http://en.wikipedia.org/wiki/Binomial_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-binomial 10 :prob 1/4 :size 20) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) n (or (:size opts) 1) p (or (:prob opts) 1/2) dist (Binomial. n p (DoubleMersenneTwister.))] @@ -1085,27 +1085,27 @@ (defn cdf-binomial -" Returns the Bionomial cdf of the given value of x. It will return a sequence +" Returns the Bionomial cdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's pbinom - Options: - :size (default 1) + Options: + :size (default 1) :prob (default 1/2) :lower-tail (default true) - See also: + See also: pdf-binomial and sample-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Binomial.html http://en.wikipedia.org/wiki/Binomial_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-binomial 10 :prob 1/4 :size 20) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) n (or (:size opts) 1) p (or (:prob opts) 1/2) lower-tail? (if (false? (:lower-tail opts)) false true) @@ -1122,22 +1122,22 @@ " Returns a sample of the given size from a Binomial distribution. Same as R's rbinom - Options: - :size (default 1) + Options: + :size (default 1) :prob (default 1/2) - See also: + See also: cdf-binomial and sample-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Binomial.html http://en.wikipedia.org/wiki/Binomial_distribution - Example: + Example: (sample-binomial 1000 :prob 1/4 :size 20) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) n (or (:size opts) 1) p (or (:prob opts) 1/2)] (if (= size 1) @@ -1146,30 +1146,30 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; POISSON DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-poisson -" Returns the Poisson pdf of the given value of x. It will return a sequence +" Returns the Poisson pdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dpois - Options: - :lambda (default 1) + Options: + :lambda (default 1) - See also: + See also: cdf-poisson and sample-poisson - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Poisson.html http://en.wikipedia.org/wiki/Poisson_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-poisson 5 :lambda 10) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) lambda (or (:lambda opts) 1) dist (Poisson. lambda (DoubleMersenneTwister.))] (if (coll? x) @@ -1179,26 +1179,26 @@ (defn cdf-poisson -" Returns the Poisson cdf of the given value of x. It will return a sequence +" Returns the Poisson cdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's ppois - Options: - :lambda (default 1) + Options: + :lambda (default 1) :lower-tail (default true) - See also: + See also: cdf-poisson and sample-poisson - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Poisson.html http://en.wikipedia.org/wiki/Poisson_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-poisson 5 :lambda 10) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) lambda (or (:lambda opts) 1) lower-tail? (if (false? (:lower-tail opts)) false true) cdf-fx (if lower-tail? @@ -1214,21 +1214,21 @@ " Returns a sample of the given size from a Poisson distribution. Same as R's rpois - Options: - :lambda (default 1) + Options: + :lambda (default 1) - See also: + See also: pdf-poisson and cdf-poisson - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/Poisson.html http://en.wikipedia.org/wiki/Poisson_distribution - Example: + Example: (sample-poisson 1000 :lambda 10) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) lambda (or (:lambda opts) 1)] (if (= size 1) (Poisson/staticNextInt lambda) @@ -1236,31 +1236,31 @@ -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; NEGATIVE BINOMIAL DISTRIBUTION FUNCTIONS -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defn pdf-neg-binomial -" Returns the Negative Binomial pdf of the given value of x. It will return a sequence +" Returns the Negative Binomial pdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dnbinom - Options: - :size (default 10) + Options: + :size (default 10) :prob (default 1/2) - See also: + See also: cdf-neg-binomial and sample-neg-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/NegativeBinomial.html http://en.wikipedia.org/wiki/Negative_binomial_distribution http://en.wikipedia.org/wiki/Probability_density_function - Example: + Example: (pdf-neg-binomial 10 :prob 1/2 :size 20) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) size (or (:size opts) 10) prob (or (:prob opts) 1/2) dist (NegativeBinomial. size prob (DoubleMersenneTwister.))] @@ -1272,26 +1272,26 @@ (defn cdf-neg-binomial -" Returns the Negative Binomial cdf of the given value of x. It will return a sequence +" Returns the Negative Binomial cdf of the given value of x. It will return a sequence of values, if x is a sequence. Same as R's dnbinom - Options: - :size (default 10) + Options: + :size (default 10) :prob (default 1/2) - See also: + See also: cdf-neg-binomial and sample-neg-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/NegativeBinomial.html http://en.wikipedia.org/wiki/Negative_binomial_distribution http://en.wikipedia.org/wiki/Cumulative_distribution_function - Example: + Example: (cdf-neg-binomial 10 :prob 1/2 :size 20) " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) size (or (:size opts) 10) prob (or (:prob opts) 1/2) lower-tail? (if (false? (:lower-tail opts)) false true) @@ -1308,22 +1308,22 @@ " Returns a sample of the given size from a Negative Binomial distribution. Same as R's rnbinom - Options: - :size (default 10) + Options: + :size (default 10) :prob (default 1/2) - See also: + See also: pdf-neg-binomial and cdf-neg-binomial - References: + References: http://incanter.org/docs/parallelcolt/api/cern/jet/random/tdouble/NegativeBinomial.html http://en.wikipedia.org/wiki/Negative_binomial_distribution - Example: + Example: (sample-neg-binomial 1000 :prob 1/2 :size 20) " ([size & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) size (or (:size opts) 10) prob (or (:prob opts) 1/2)] (if (= size 1) @@ -1336,13 +1336,13 @@ ;; EMPIRICAL DISTRIBUTION ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -(defn cdf-empirical +(defn cdf-empirical " Returns a step-function representing the empirical cdf of the given data. Equivalent to R's ecdf function. - The following description is from the ecdf help in R: The e.c.d.f. - (empirical cumulative distribution function) Fn is a step function - with jumps i/n at observation values, where i is the number of tied + The following description is from the ecdf help in R: The e.c.d.f. + (empirical cumulative distribution function) Fn is a step function + with jumps i/n at observation values, where i is the number of tied observations at that value. Missing values are ignored. For observations 'x'= (x1,x2, ... xn), Fn is the fraction of @@ -1354,8 +1354,8 @@ Examples: (use '(incanter core stats charts)) - (def exam1 [192 160 183 136 162 165 181 188 150 163 192 164 184 - 189 183 181 188 191 190 184 171 177 125 192 149 188 + (def exam1 [192 160 183 136 162 165 181 188 150 163 192 164 184 + 189 183 181 188 191 190 184 171 177 125 192 149 188 154 151 159 141 171 153 169 168 168 157 160 190 166 150]) ;; the ecdf function returns an empirical cdf function for the given data @@ -1365,7 +1365,7 @@ (view (scatter-plot exam1 (map ecdf exam1))) - + " ([x] (fn [t] (div (sum (indicator #(<= % t) x)) (count x))))) @@ -1378,7 +1378,7 @@ (defn mean " - Returns the mean of the data, x. + Returns the mean of the data, x. Examples: (mean (sample-normal 100)) @@ -1387,14 +1387,14 @@ http://incanter.org/docs/parallelcolt/api/cern/jet/stat/tdouble/DoubleDescriptive.html http://en.wikipedia.org/wiki/Mean -" +" ([x] (let [xx (to-list x)] (DoubleDescriptive/mean (DoubleArrayList. (double-array xx)))))) -(defn variance +(defn variance " Returns the sample variance of the data, x. Equivalent to R's var function. @@ -1410,9 +1410,9 @@ -(defn covariance +(defn covariance " - Returns the sample covariance of x and y. + Returns the sample covariance of x and y. Examples: ;; create some data that covaries @@ -1431,20 +1431,20 @@ xx (to-list x) yy (to-list y) ] - (DoubleDescriptive/covariance + (DoubleDescriptive/covariance (DoubleArrayList. (double-array xx)) (DoubleArrayList. (double-array yy))))) ([mat] (let [n (ncol mat)] - (matrix - (for [i (range n) j (range n)] + (matrix + (for [i (range n) j (range n)] (covariance (sel mat true i) (sel mat true j))) n)))) -(defn sd +(defn sd " - Returns the sample standard deviation of the data, x. Equivalent to + Returns the sample standard deviation of the data, x. Equivalent to R's sd function. Examples: @@ -1453,7 +1453,7 @@ References: http://incanter.org/docs/parallelcolt/api/cern/jet/stat/tdouble/DoubleDescriptive.html http://en.wikipedia.org/wiki/Standard_deviation -" +" ([x] ;; population sd, not the sample sd ;(DoubleDescriptive/sampleStandardDeviation (length x) (variance x)))) @@ -1464,8 +1464,8 @@ (defn correlation " - Returns the sample correlation of x and y, or the correlation - matrix of the given matrix. + Returns the sample correlation of x and y, or the correlation + matrix of the given matrix. Examples: @@ -1478,18 +1478,18 @@ xx (to-list x) yy (to-list y) ] - (DoubleDescriptive/correlation + (DoubleDescriptive/correlation (DoubleArrayList. (double-array xx)) (sd x) (DoubleArrayList. (double-array yy)) (sd y)))) ([mat] - (div (covariance mat) + (div (covariance mat) (sqrt (mmult (diag (covariance mat)) (trans (diag (covariance mat)))))))) -(defn median +(defn median " - Returns the median of the data, x. + Returns the median of the data, x. Examples: (median (sample-normal 100)) @@ -1498,18 +1498,18 @@ http://incanter.org/docs/parallelcolt/api/cern/jet/stat/tdouble/DoubleDescriptive.html http://en.wikipedia.org/wiki/Median -" +" ([x] (let [xx (sort (to-list x))] (DoubleDescriptive/median (DoubleArrayList. (double-array xx)))))) - -(defn kurtosis + +(defn kurtosis " - Returns the kurtosis of the data, x. \"Kurtosis is a measure of the \"peakedness\" - of the probability distribution of a real-valued random variable. Higher kurtosis - means more of the variance is due to infrequent extreme deviations, as opposed to + Returns the kurtosis of the data, x. \"Kurtosis is a measure of the \"peakedness\" + of the probability distribution of a real-valued random variable. Higher kurtosis + means more of the variance is due to infrequent extreme deviations, as opposed to frequent modestly-sized deviations.\" (Wikipedia) Examples: @@ -1526,9 +1526,9 @@ -(defn skewness +(defn skewness " - Returns the skewness of the data, x. \"Skewness is a measure of the asymmetry + Returns the skewness of the data, x. \"Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable.\" (Wikipedia) Examples: @@ -1545,9 +1545,9 @@ -(defn quantile +(defn quantile " - Returns the quantiles of the data, x. By default it returns the min, + Returns the quantiles of the data, x. By default it returns the min, 25th-percentile, 50th-percentile, 75th-percentile, and max value. Options: @@ -1564,10 +1564,10 @@ " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) _x (if (matrix? x) (to-list x) x) data (DoubleArrayList. (double-array (sort _x))) - probs (cond + probs (cond (number? (:probs opts)) (:probs opts) (coll? (:probs opts)) @@ -1580,8 +1580,8 @@ -(defn sample -" Returns a sample of the given size from the given collection. If replacement +(defn sample +" Returns a sample of the given size from the given collection. If replacement is set to false it returns a set, otherwise it returns a list. Arguments: @@ -1599,7 +1599,7 @@ " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) size (or (:size opts) (count x)) replacement? (if (false? (:replacement opts)) false true) max-idx (dec (count x))] @@ -1629,32 +1629,32 @@ Options: :size -- the number of bootstrap samples to return - :smooth -- (default false) smoothing option - :smooth-sd -- (default (/ (sqrt (count data)))) determines the standard + :smooth -- (default false) smoothing option + :smooth-sd -- (default (/ (sqrt (count data)))) determines the standard deviation of the noise to use for smoothing - :replacement -- (default true) determines if sampling of the data + :replacement -- (default true) determines if sampling of the data should be done with replacement References: 1. Clifford E. Lunneborg, Data Analysis by Resampling Concepts and Applications, 2000, pages 105-117 2. http://en.wikipedia.org/wiki/Bootstrapping_(statistics) - + Examples: - + ;; example from Data Analysis by Resampling Concepts and Applications ;; Clifford E. Lunneborg (pages 119-122) (use '(incanter core stats charts)) ;; weights (in grams) of 50 randomly sampled bags of preztels - (def weights [464 447 446 454 450 457 450 442 - 433 452 449 454 450 438 448 449 - 457 451 456 452 450 463 464 453 - 452 447 439 449 468 443 433 460 - 452 447 447 446 450 459 466 433 - 445 453 454 446 464 450 456 456 + (def weights [464 447 446 454 450 457 450 442 + 433 452 449 454 450 438 448 449 + 457 451 456 452 450 463 464 453 + 452 447 439 449 468 443 433 460 + 452 447 447 446 450 459 466 433 + 445 453 454 446 464 450 456 456 447 469]) ;; calculate the sample median, 450 @@ -1689,17 +1689,17 @@ ;; Newcomb's speed of light data (use '(incanter core stats charts)) - + ;; A numeric vector giving the Third Series of measurements of the ;; passage time of light recorded by Newcomb in 1882. The given ;; values divided by 1000 plus 24 give the time in millionths of a ;; second for light to traverse a known distance. The 'true' value is ;; now considered to be 33.02. - (def speed-of-light [28 -44 29 30 24 28 37 32 36 27 26 28 29 - 26 27 22 23 20 25 25 36 23 31 32 24 27 - 33 16 24 29 36 21 28 26 27 27 32 25 28 - 24 40 21 31 32 28 26 30 27 26 24 32 29 + (def speed-of-light [28 -44 29 30 24 28 37 32 36 27 26 28 29 + 26 27 22 23 20 25 25 36 23 31 32 24 27 + 33 16 24 29 36 21 28 26 27 27 32 25 28 + 24 40 21 31 32 28 26 30 27 26 24 32 29 34 -2 25 19 36 29 30 22 28 33 39 25 16 23]) ;; view histogram of data to see outlier observations @@ -1714,11 +1714,11 @@ (view (histogram smooth-samp :density true :nbins 30)) (mean smooth-samp) (quantile smooth-samp :probs [0.025 0.975]) - + " ([data statistic & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) size (:size opts) replacement (if (false? (:replacement opts)) false true) n (or (:n opts) (count data)) @@ -1742,13 +1742,13 @@ (if smooth? (plus samp (sample-normal samp-size :sd smooth-sd)) samp)))) - -(defn cumulative-mean - " Returns a sequence of cumulative means for the given collection. For instance + +(defn cumulative-mean + " Returns a sequence of cumulative means for the given collection. For instance The first value equals the first value of the argument, the second value is the mean of the first two arguments, the third is the mean of the first three arguments, etc. @@ -1773,29 +1773,29 @@ :fun (defaul minus) the function used to sweep the stat out Example: - + (use '(incanter core stats)) - + (def x (sample-normal 30 :mean 10 :sd 5)) (sweep x) ;; center the data around mean (sweep x :stat sd :fun div) ;; divide data by its sd " ([x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) stat-fn (or (:stat opts) mean) fun (or (:fun opts) minus) stat (stat-fn x)] (fun x stat)))) - + (defn permute -" If provided a single argument, returns a permuted version of the +" If provided a single argument, returns a permuted version of the given collection. (perm x) is the same as (sample x). - If provided two arguments, returns two lists that are permutations + If provided two arguments, returns two lists that are permutations across the given collections. In other words, each of the new collections will contain elements from both of the given collections. Useful for permutation tests or randomization tests. @@ -1804,7 +1804,7 @@ (permute (range 10)) (permute (range 10) (range 10 20)) - + " ([x] @@ -1819,7 +1819,7 @@ -(defn sample-permutations +(defn sample-permutations " If provided a two arguments (n x), it returns a list of n permutations of x. If provided three (n x y) arguments, returns a list with two with n permutations of each arguments, where each permutation is drawn from the pooled arguments. @@ -1838,42 +1838,42 @@ ;; extended example with plant-growth data (use '(incanter core stats datasets charts)) - + ;; load the plant-growth dataset (def data (to-matrix (get-dataset :plant-growth))) - + ;; break the first column of the data into groups based on treatment (second column). (def groups (group-by data 1 :cols 0)) - + ;; define a function for the statistic of interest (defn means-diff [x y] (minus (mean x) (mean y))) - + ;; calculate the difference in sample means between the two groups (def samp-mean-diff (means-diff (first groups) (second groups))) ;; 0.371 - + ;; create 500 permuted versions of the original two groups (def permuted-groups (sample-permutations 1000 (first groups) (second groups))) - + ;; calculate the difference of means of the 500 samples (def permuted-means-diffs1 (map means-diff (first permuted-groups) (second permuted-groups))) - + ;; use an indicator function that returns 1 when the randomized means diff is greater ;; than the original sample mean, and zero otherwise. Then take the mean of this sequence ;; of ones and zeros. That is the proportion of times you would see a value more extreme ;; than the sample mean (i.e. the p-value). (mean (indicator #(> % samp-mean-diff) permuted-means-diffs1)) ;; 0.088 - + ;; calculate the 95% confidence interval of the null hypothesis. If the - ;; sample difference in means is outside of this range, that is evidence + ;; sample difference in means is outside of this range, that is evidence ;; that the two means are statistically significantly different. (quantile permuted-means-diffs1 :probs [0.025 0.975]) ;; (-0.606 0.595) ;; Plot a histogram of the permuted-means-diffs using the density option, - ;; instead of the default frequency, and then add a normal pdf curve with + ;; instead of the default frequency, and then add a normal pdf curve with ;; the mean and sd of permuted-means-diffs data for a visual comparison. (doto (histogram permuted-means-diffs1 :density true) - (add-lines (range -1 1 0.01) (pdf-normal (range -1 1 0.01) - :mean (mean permuted-means-diffs1) + (add-lines (range -1 1 0.01) (pdf-normal (range -1 1 0.01) + :mean (mean permuted-means-diffs1) :sd (sd permuted-means-diffs1))) view) @@ -1883,19 +1883,19 @@ (def samp-mean-diff (means-diff (first groups) (last groups))) ;; -0.4939 (mean (indicator #(< % samp-mean-diff) permuted-means-diffs2)) ;; 0.022 (quantile permuted-means-diffs2 :probs [0.025 0.975]) ;; (-0.478 0.466) - + ;; compare the means of treatment 1 and treatment 2 (def permuted-groups (sample-permutations 1000 (second groups) (last groups))) (def permuted-means-diffs3 (map means-diff (first permuted-groups) (second permuted-groups))) (def samp-mean-diff (means-diff (second groups) (last groups))) ;; -0.865 (mean (indicator #(< % samp-mean-diff) permuted-means-diffs3)) ;; 0.002 (quantile permuted-means-diffs3 :probs [0.025 0.975]) ;; (-0.676 0.646) - + (doto (box-plot permuted-means-diffs1) (add-box-plot permuted-means-diffs2) (add-box-plot permuted-means-diffs3) view) - + Further Reading: http://en.wikipedia.org/wiki/Resampling_(statistics) @@ -1905,7 +1905,7 @@ (loop [samp '() i 0] (if (= i n) samp - (recur + (recur (conj samp (sample x)) (inc i))))) ([n x y] @@ -1926,7 +1926,7 @@ (defn linear-model -" +" Returns the results of performing a OLS linear regression of y on x. Arguments: @@ -1970,7 +1970,7 @@ (def x1 (range 0.0 3 0.1)) (view (xy-plot x1 (cdf-f x1 :df1 4 :df2 144))) - + References: http://en.wikipedia.org/wiki/OLS_Regression @@ -1978,7 +1978,7 @@ " ([y x & options] - (let [opts (when options (apply assoc {} options)) + (let [opts (when options (apply assoc {} options)) intercept? (if (false? (:intercept opts)) false true) _x (if intercept? (bind-columns (replicate (nrow x) 1) x) x) xtx (mmult (trans _x) _x) @@ -1986,7 +1986,7 @@ xty (mmult (trans _x) y) coefs (if (and (number? xtxi) (number? xty)) (* xtxi xty) - (to-list (if (or (number? xtxi) (number? xty)) + (to-list (if (or (number? xtxi) (number? xty)) (mult xtxi xty) (mmult xtxi xty)))) fitted (to-list (if (number? coefs) @@ -2013,11 +2013,11 @@ t-probs (mult 2 (cdf-t (abs t-tests) :df df2 :lower-tail false)) t-95 (mult (quantile-t 0.975 :df df2) std-errors) coefs-ci (if (number? std-errors) - [(plus coefs t-95) + [(plus coefs t-95) (minus coefs t-95)] - (partition 2 - (interleave - (minus coefs t-95) + (partition 2 + (interleave + (minus coefs t-95) (plus coefs t-95)))) ] (with-meta @@ -2042,7 +2042,7 @@ :coef-var coef-var :r-square r-square :adj-r-square adj-r-square - } + } {:type ::linear-model})))) @@ -2061,8 +2061,8 @@ :conf-level (default 0.95) for returned confidence interval Examples: - - (t-test (range 1 11) :mu 0) + + (t-test (range 1 11) :mu 0) (t-test (range 1 11) :mu 0 :alternative :less) (t-test (range 1 11) :mu 0 :alternative :greater) @@ -2078,9 +2078,9 @@ " ([x & options] (let [opts (when options (apply assoc {} options)) - y (or (:y opts) nil) + y (or (:y opts) nil) one-sample? (nil? y) - mu (or (:mu opts) + mu (or (:mu opts) (if y (mean y) 0)) paired? (if (true? (:paired opts)) true false) var-equal? (if (true? (:var-equal opts)) true false) @@ -2129,14 +2129,14 @@ [(if (= alternative :less) Double/NEGATIVE_INFINITY (+ x-mean (/ (* qt (sqrt x-var)) (sqrt n1)))) - (if (= alternative :greater) + (if (= alternative :greater) Double/POSITIVE_INFINITY (- x-mean (/ (* qt (sqrt x-var)) (sqrt n1))))] ;; two-sample confidence interval [(if (= alternative :less) - Double/NEGATIVE_INFINITY + Double/NEGATIVE_INFINITY (+ (- x-mean y-mean) (* qt (sqrt (+ (/ x-var n1) (/ y-var n2)))))) - (if (= alternative :greater) + (if (= alternative :greater) Double/POSITIVE_INFINITY (- (- x-mean y-mean) (* qt (sqrt (+ (/ x-var n1) (/ y-var n2))))))]) }))) @@ -2152,14 +2152,14 @@ :levels -- a sequence of sequences, where each sequence list the levels (possible values) of the corresponding column of x. - :margins -- a sequence of sequences, where each sequence - represents the marginal total for each level + :margins -- a sequence of sequences, where each sequence + represents the marginal total for each level of the corresponding column of x. :counts -- a hash-map, where vectors of unique combinations of the cross-tabulated levels are the keys and the values are the total count of each combination. :N -- the grand-total for the contingency table - + Examples: @@ -2172,26 +2172,26 @@ (tabulate (sel math-prog :cols [1 2])) - (def data (matrix [[1 0 1] - [1 1 1] - [1 1 1] - [1 0 1] - [0 0 0] - [1 1 1] - [1 1 1] - [1 0 1] + (def data (matrix [[1 0 1] + [1 1 1] + [1 1 1] + [1 0 1] + [0 0 0] + [1 1 1] + [1 1 1] + [1 0 1] [1 1 0]])) (tabulate data) - (def data (matrix [[1 0] - [1 1] - [1 1] - [1 0] - [0 0] - [1 1] - [1 1] - [1 0] + (def data (matrix [[1 0] + [1 1] + [1 1] + [1 0] + [0 0] + [1 1] + [1 1] + [1 0] [1 1]])) (tabulate data) @@ -2208,7 +2208,7 @@ marg (let [lvl (sel _x :rows i :cols j)] (recur (let [cnt (get marg lvl)] - (if cnt + (if cnt (assoc marg lvl (inc cnt)) (assoc marg lvl 1))) (inc i)))))) @@ -2228,40 +2228,40 @@ {:counts counts :margins margins :table (when (= p 2) - (matrix (for [r (first levels) - c (second levels)] - (let [c (get counts (to-list (trans [r c])))] - (if c c 0))) + (matrix (for [r (first levels) + c (second levels)] + (let [c (get counts (to-list (trans [r c])))] + (if c c 0))) (second n-levels))) :n-vars p :N (reduce + (vals (first margins))) :n-levels n-levels :levels levels}))) - + (defn chisq-test -" +" Performs chi-squared contingency table tests and goodness-of-fit tests. - + If the optional argument :y is not provided then a goodness-of-fit test - is performed. In this case, the hypothesis tested is whether the - population probabilities equal those in :probs, or are all equal if + is performed. In this case, the hypothesis tested is whether the + population probabilities equal those in :probs, or are all equal if :probs is not given. If :y is provided, it must be a sequence of integers that is the - same length as x. A contingency table is computed from x and :y. + same length as x. A contingency table is computed from x and :y. Then, Pearson's chi-squared test of the null hypothesis that the joint distribution of the cell counts in a 2-dimensional contingency table is the product of the row and column marginals is performed. By default the Yates' continuity correction for 2x2 contingency - tables is performed, this can be disabled by setting the :correct + tables is performed, this can be disabled by setting the :correct option to false. - + Options: :x -- a sequence of numbers. @@ -2287,7 +2287,7 @@ (chisq-test :table (trans table)) ;; throws exception (chisq-test :x [1 0 0 0 1 1 1 0 0 1 0 0 1 1 1 1]) ;; 0.25 - + (use '(incanter core stats datasets)) (def math-prog (to-matrix (get-dataset :math-prog))) (def x (sel math-prog :cols 1)) @@ -2316,7 +2316,7 @@ (def probs [0.40 0.20 0.20 0.19 0.01]) (chisq-test :table table :probs probs) ;; X-sq = 5.7947, df = 4, p-value = 0.215 - ;; use frequencies instead of probabilities + ;; use frequencies instead of probabilities (def freq [40 20 20 15 5]) (chisq-test :table table :freq freq) ;; X-sq = 9.9901, df = 4, p-value = 0.04059 @@ -2335,47 +2335,47 @@ y (:y opts) table? (if (:table opts) true false) xtab (when (or x y) - (if y - (tabulate (bind-columns x y)) + (if y + (tabulate (bind-columns x y)) (tabulate x))) - table (cond - table? - (:table opts) + table (cond + table? + (:table opts) (and x y) (:table xtab)) - two-samp? (if (or (and x y) - (and table? - (and (> (nrow table) 1) (> (ncol table) 1)))) + two-samp? (if (or (and x y) + (and table? + (and (> (nrow table) 1) (> (ncol table) 1)))) true false) r-levels (if table? - (range (nrow table)) + (range (nrow table)) (first (:levels xtab))) - c-levels (if table? - (range (ncol table)) + c-levels (if table? + (range (ncol table)) (second (:levels xtab))) r-margins (if table? (if two-samp? (apply hash-map (interleave r-levels (map sum (trans table)))) - (if (> (nrow table) 1) - (to-list table) + (if (> (nrow table) 1) + (to-list table) (throw (Exception. "One dimensional tables must have only a single column")))) (second (:margins xtab))) - c-margins (if table? + c-margins (if table? (if two-samp? (apply hash-map (interleave c-levels (map sum table))) 0) (first (:margins xtab))) - + counts (if two-samp? (vectorize table) table) - N (if table? - (sum counts) + N (if table? + (sum counts) (:N xtab)) n (when (not two-samp?) (count r-levels)) df (if two-samp? (* (dec (nrow table)) (dec (ncol table))) (dec n)) probs (when (not two-samp?) - (cond - (:probs opts) - (:probs opts) + (cond + (:probs opts) + (:probs opts) (:freq opts) (div (:freq opts) (sum (:freq opts))) :else @@ -2404,11 +2404,11 @@ -(defn principal-components +(defn principal-components " Performs a principal components analysis on the given data matrix. Equivalent to R's prcomp function. - + Returns: A map with the following fields: :std-dev -- the standard deviations of the principal compoenents @@ -2420,7 +2420,7 @@ Examples: - + (use '(incanter core stats charts datasets)) ;; load the iris dataset (def iris (to-matrix (get-dataset :iris))) @@ -2432,9 +2432,9 @@ ;; project the first four dimension of the iris data onto the first ;; two principal components - (def x1 (mmult (sel iris :cols (range 4)) pc1)) - (def x2 (mmult (sel iris :cols (range 4)) pc2)) - + (def x1 (mmult (sel iris :cols (range 4)) pc1)) + (def x2 (mmult (sel iris :cols (range 4)) pc2)) + ;; now plot the transformed data, coloring each species a different color (doto (scatter-plot (sel x1 :rows (range 50)) (sel x2 :rows (range 50)) :x-label \"PC1\" :y-label \"PC2\" :title \"Iris PCA\") @@ -2444,7 +2444,7 @@ ;; alternatively, the :group-by option can be used in scatter-plot - (view (scatter-plot x1 x2 + (view (scatter-plot x1 x2 :group-by (sel iris :cols 4) :x-label \"PC1\" :y-label \"PC2\" :title \"Iris PCA\")) @@ -2476,17 +2476,17 @@ (tabulate (detabulate :table table)) ;; example 2 - (def data (matrix [[1 0] - [1 1] - [1 1] - [1 0] - [0 0] - [1 1] - [1 1] - [1 0] + (def data (matrix [[1 0] + [1 1] + [1 1] + [1 0] + [0 0] + [1 1] + [1 1] + [1 0] [1 1]])) (tabulate data) - + (tabulate (detabulate :table (:table (tabulate data)))) " @@ -2495,12 +2495,8 @@ table (:table opts) row-labels (when table (or (:row-labels opts) (range (nrow table)))) col-labels (when table (or (:col-labels opts) (range (ncol table)))) - data (apply bind-rows + data (apply bind-rows (apply concat (for [r row-labels c col-labels] (repeat (sel table :rows r :cols c) [r c]))))] data))) - - - - src/incanter/stats.clj @@ -1,7 +1,7 @@ #!/usr/bin/env bin/clj -;;; test-incanter.clj -- Unit tests of Incanter functions +;;; test-incanter.clj -- Unit tests of Incanter functions ;; by David Edgar Liebke http://incanter.org ;; March 12, 2009 @@ -19,7 +19,7 @@ (ns tests.runtests - (:use [clojure.test :only (run-tests)] + (:use [clojure.test :only (run-tests)] (tests test-cases) (incanter core io stats charts bayes))) tests/runtests.clj 652fadad4777fbcd37026e3ba61dabeddd7317e1 cd4a3601fbbfe91f7b62133d7173ad0aafe657e9 David Edgar Liebke liebke@gmail.com http://github.com/liebke/incanter/commit/08f68bc55b83f22053a8b733ef9bebdc7793af6f 08f68bc55b83f22053a8b733ef9bebdc7793af6f 2009-08-13T06:09:19-07:00 2009-08-13T06:09:19-07:00 fixed .gitignore merge 8f7858e85be7191cb54ee187104eebf3dab4a6ba David Edgar Liebke liebke@gmail.com