diff --git a/examples/brian_drawert/param_search.m b/examples/brian_drawert/param_search.m
deleted file mode 100644
index 9afef42b43b59..0000000000000
--- a/examples/brian_drawert/param_search.m
+++ /dev/null
@@ -1,223 +0,0 @@
-function out=param_search()
-close all;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Smax=10000;
-V=1;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Nfrogs_range = [5 10 15 20 25 30 40 50];
-gamma_range = [ 0.0005 0.001 0.0025 0.005 0.01];
-eta_range = [5 7.5 10 12.5 15 17.5];
-nu_range = [0.15 0.3 0.5 0.7];
-f_range = [ 0.05 0.1 0.25 0.5];
-sigma_range = [0.1 0.2 0.275 0.4 0.5];
-mu_range = [0.01 0.1 1 10 100];
-%%%
-nu1_range = [1.3863e-04 6.9315e-05 2.8768e-05];  % = -log( [ .25 .5 .75 ] )/Smax
-sigma1_range = [5e-06 5e-05 5e-04];  % sigma0/(4*Smax) = 6.8750e-06, iff sigma0=0.275
-phi_range = [ 0.1 0.5 1];
-Ar_range = [ 1e-5 1e-4 1e-3 1e-2];
-Ac_range = [ 0.001 0.01 0.1];
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-num_runs = length(Nfrogs_range)*length(gamma_range)*length(eta_range)*length(nu_range)*length(f_range)*length(sigma_range)*length(mu_range);
-num_runs2=length(Nfrogs_range)*length(gamma_range)*length(sigma_range)*length(nu1_range)*length(sigma1_range)*length(phi_range)*length(Ar_range)*length(Ac_range)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Nfrogs = {};
-gamma = {};
-eta = {};
-nu= {};
-f = {};
-sigma = {};
-mu = {};
-%%%
-nu1=0;
-sigma1=0;
-phi=0;
-Ar=0;
-Ac=0;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-for mu_ndx=1:length(mu_range)
-    for sigma_ndx=1:length(sigma_range)
-        for f_ndx=1:length(f_range)
-            for nu_ndx=1:length(nu_range)
-                for eta_ndx=1:length(eta_range)
-                    for gamma_ndx=1:length(gamma_range)
-                        for Nfrogs_ndx=1:length(Nfrogs_range)
-                            run_param_set();
-                            %if this parameter set works, check non-linear
-                            if (eta_ndx==6 && nu_ndx==1 && f_ndx==2 && mu_ndx==3)
-                            for nu1_ndx=0:length(nu1_range)
-                                for sigma1_ndx=0:length(sigma1_range)
-                                    for phi_ndx = 0:length(phi_range)
-                                        if(phi_ndx>0)
-                                        for Ar_ndx = 1:length(Ar_range)
-                                            for Ac_ndx=1:length(Ac_range)
-                            run_param_set_w_antibodies();
-                                            end
-                                        end
-                                        end
-                                    end
-                                end
-                            end
-                            end
-                        end
-                    end
-                end
-            end
-        end
-    end
-end
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    function out = run_param_set_w_antibodies()
-        nu1=0;
-        if(nu1_ndx>0),nu1=nu1_range(nu1_ndx);end
-        sigma1=0;
-        if(sigma1_ndx>0),sigma1=sigma1_range(sigma1_ndx);end
-        phi=0;
-        Ar=0;
-        Ac=0;
-        Ar_ndx2=0;
-        Ac_ndx2=0;
-        if(phi_ndx>0)
-            phi=phi_range(phi_ndx);
-            Ar=Ar_range(Ar_ndx);
-            Ac=Ac_range(Ac_ndx);
-            Ar_ndx2=Ar_ndx;
-            Ac_ndx2=Ac_ndx;
-        end
-        %%%
-        filename=sprintf('data/param_search_data_%i_%i_%i_%i_%i_%i_%i--%i_%i_%i_%i_%i.mat',...
-            mu_ndx,sigma_ndx,f_ndx,nu_ndx,eta_ndx,gamma_ndx,Nfrogs_ndx,...
-            nu1_ndx,sigma1_ndx,phi_ndx,Ar_ndx2,Ac_ndx2);
-        fid=fopen(filename);
-        if(fid~=-1)
-            fprintf('skipping: %s\n',filename);
-            fclose(fid);
-            return;
-        end
-        fprintf('%s\n',filename);
-                            %out = run_param_set();
-                            %return;
-        %%%%%%
-        Nfrogs = Nfrogs_range(Nfrogs_ndx);
-        gamma = gamma_range(gamma_ndx);
-        eta = eta_range(eta_ndx);
-        nu=nu_range(nu_ndx);
-        f = f_range(f_ndx);
-        sigma = sigma_range(sigma_ndx);
-        mu = mu_range(mu_ndx);
-        %%%%%%
-        out.param_set={Nfrogs,gamma,eta,nu,f,sigma,mu,nu1,sigma1,phi,Ar,Ac};
-        %%%%%%
-        end_time=1000;
-        %%%%%%
-        % Find: days to reach_Smax
-        if(phi>0)
-            y0 = zeros(2*Nfrogs+1,1);
-        else
-            y0 = zeros(Nfrogs+1,1);
-        end
-        y0(1)=100;  %first frog is infected
-        options = odeset('Events',@events,'Refine',4);
-        [t,y,te,ye,ie] = ode23(@dydt,[0 end_time],y0,options);
-        out.days_to_reach_Smax = t(end);
-        %%%%%%
-        % Find: frac_surving_frogs
-        [t2,y2] = ode23(@dydt,[0 end_time],y0);
-        out.frac_surviving_frogs = sum( y2(end,1:Nfrogs) < Smax )/Nfrogs;
-        out.y=y;
-        out.t=t;
-        out.y2=y2;
-        out.t2=t2;
-        %%%%%%
-        save(filename,'-STRUCT','out')
-    end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    function out = run_param_set()
-        nu1=0;
-        sigma1=0;
-        phi=0;
-        Ar=0;
-        Ac=0;
-        filename=sprintf('data/param_search_data_%i_%i_%i_%i_%i_%i_%i.mat',mu_ndx,sigma_ndx,f_ndx,nu_ndx,eta_ndx,gamma_ndx,Nfrogs_ndx);
-        fid=fopen(filename);
-        if(fid~=-1)
-            fprintf('skipping: %s\n',filename);
-            fclose(fid);
-            return;
-        end
-        fprintf('%s\n',filename);
-                            %out = run_param_set();
-                            %return;
-        %%%%%%
-        Nfrogs = Nfrogs_range(Nfrogs_ndx);
-        gamma = gamma_range(gamma_ndx);
-        eta = eta_range(eta_ndx);
-        nu=nu_range(nu_ndx);
-        f = f_range(f_ndx);
-        sigma = sigma_range(sigma_ndx);
-        mu = mu_range(mu_ndx);
-        %%%%%%
-        out.param_set={Nfrogs,gamma,eta,nu,f,sigma,mu};
-        %%%%%%
-        end_time=1000;
-        %%%%%%
-        % Find: days to reach_Smax
-        if(phi>0)
-            y0 = zeros(2*Nfrogs+1,1);
-        else
-            y0 = zeros(Nfrogs+1,1);
-        end
-        y0(1)=100;  %first frog is infected
-        options = odeset('Events',@events,'Refine',4);
-        [t,y,te,ye,ie] = ode23(@dydt,[0 end_time],y0,options);
-        out.days_to_reach_Smax = t(end);
-        %%%%%%
-        % Find: frac_surving_frogs
-        [t2,y2] = ode23(@dydt,[0 end_time],y0);
-        out.frac_surviving_frogs = sum( y2(end,1:Nfrogs) < Smax )/Nfrogs;
-        out.y=y;
-        out.t=t;
-        out.y2=y2;
-        out.t2=t2;
-        %%%%%%
-        save(filename,'-STRUCT','out')
-    end
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    function d = dydt(t,y)
-        z = y;
-        for i=1:Nfrogs
-            if(y(i)>Smax)
-                z(i)=0;
-                d(i) = 0;
-            else
-                nu0 = nu*exp(-nu1*y(i));
-                if(phi>0)
-                    nu0=nu*exp(-phi*y(i+Nfrogs+1));
-                end
-                sigma0=sigma+sigma1*y(i);
-                beta = gamma*nu0;
-                d(i) = beta/V*y(Nfrogs+1) + (eta*nu0*f-sigma0)*y(i);
-            end
-        end
-        d(Nfrogs+1) = sum(z(1:Nfrogs)*eta*(1-f)-(gamma/V)*y(Nfrogs+1)) - mu*y(Nfrogs+1);
-        if(phi>0)
-            for i=1:Nfrogs
-                d(i+Nfrogs+1) = Ar*y(i) - Ac*y(i+Nfrogs+1);
-            end
-        end
-        d=d';
-     end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    function [value,isterminal,direction]= events(t,y)
-        value(1:Nfrogs) = Smax-y(1:Nfrogs);
-        value(Nfrogs+1:end)=1;
-        isterminal=ones(size(value));
-        direction=zeros(size(value));
-    end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-end%function param_search()
\ No newline at end of file
diff --git a/examples/select.j b/examples/select.j
deleted file mode 100644
index b3d860407e195..0000000000000
--- a/examples/select.j
+++ /dev/null
@@ -1,33 +0,0 @@
-# http://en.wikipedia.org/wiki/Selection_algorithm#Partition-based_general_selection_algorithm
-
-function _jl_quickselect(a::AbstractVector, k::Integer, lo::Integer, hi::Integer)
-    while true
-        #partition(a, lo, hi, m)
-        i, j = lo, hi
-        pivot = a[div(lo+hi, 2)]
-        while i <= j
-            while a[i] < pivot; i += 1; end
-            while pivot < a[j]; j -= 1; end
-            if i <= j
-                a[i], a[j] = a[j], a[i]
-                i += 1
-                j -= 1
-            end
-        end
-
-        pivot_new = j + 1
-        pivot_dist = pivot_new - lo + 1
-        if pivot_dist == k
-            return a[pivot_new]
-        elseif k < pivot_dist
-            hi = pivot_new - 1
-        else
-            k = k - pivot_dist
-            lo = pivot_new + 1
-        end
-
-    end
-end
-
-select(a::AbstractVector, k::Integer) = _jl_quickselect(copy(a), k, 1, length(a))
-select!(a::AbstractVector, k::Integer) = _jl_quickselect(a, k, 1, length(a))
diff --git a/j/linalg.j b/j/linalg.j
index 173e60431cbf9..c11ea4ddbc46a 100644
--- a/j/linalg.j
+++ b/j/linalg.j
@@ -237,18 +237,6 @@ function trace{T}(A::Matrix{T})
     return t
 end
 
-mean(V::Vector) = sum(V) / length(V)
-
-function std(V::Vector)
-    n = numel(V)
-    m = mean(V)
-    s = 0.0
-    for i=1:n
-        s += (V[i] - m)^2
-    end
-    return sqrt(s/(n-1))
-end
-
 kron(a::Vector, b::Vector) = [ a[i]*b[j] | i=1:length(a), j=1:length(b) ]
 
 function kron{T,S}(a::Matrix{T}, b::Matrix{S})
diff --git a/j/sort.j b/j/sort.j
index b214dd8f1fd34..107a968c4f983 100644
--- a/j/sort.j
+++ b/j/sort.j
@@ -313,3 +313,41 @@ function issorted(v::AbstractVector)
   end
   return true
 end
+
+function _jl_quickselect(a::AbstractVector, k::Integer, lo::Integer, hi::Integer)
+    if k < lo || k > hi; error("k is out of bounds"); end
+
+    while true
+
+        if lo == hi; return a[lo]; end
+
+        # pivot_ind = partition(a, lo, hi)
+        i, j = lo, hi
+        pivot = a[randival(lo,hi)]
+        while i < j
+            while a[i] < pivot; i += 1; end
+            while a[j] > pivot; j -= 1; end
+            if a[i] == a[j]
+                i += 1
+            elseif i < j
+                a[i], a[j] = a[j], a[i]
+            end
+        end
+        pivot_ind = j
+
+        length = pivot_ind - lo + 1
+        if k == length
+            return a[pivot_ind]
+        elseif k <  length
+            hi = pivot_ind - 1
+        else
+            lo = pivot_ind + 1
+            k = k - length
+        end
+
+    end # while true...
+
+end
+
+select(a::AbstractVector, k::Integer) = _jl_quickselect(copy(a), k, 1, length(a))
+select!(a::AbstractVector, k::Integer) = _jl_quickselect(a, k, 1, length(a))
diff --git a/j/statistics.j b/j/statistics.j
new file mode 100644
index 0000000000000..b72f5ee4b4838
--- /dev/null
+++ b/j/statistics.j
@@ -0,0 +1,13 @@
+mean(V::AbstractVector) = sum(V) / length(V)
+
+function std(V::AbstractVector)
+    n = numel(V)
+    m = mean(V)
+    s = 0.0
+    for i=1:n
+        s += (V[i] - m)^2
+    end
+    return sqrt(s/(n-1))
+end
+
+median(V::AbstractVector) = select(V, div(numel(V),2))
diff --git a/j/sysimg.j b/j/sysimg.j
index 36ff0ba0ab9a7..c8019670a1cc4 100644
--- a/j/sysimg.j
+++ b/j/sysimg.j
@@ -98,6 +98,7 @@ load("math.j")
 load("math_libm.j")
 load("sort.j")
 load("combinatorics.j")
+load("statistics.j")
 
 # linear algebra
 load("linalg.j")