diff --git a/src/array/_globals.kwargs.ts b/src/array/_globals.kwargs.ts
index 8bd03f6..1110006 100644
--- a/src/array/_globals.kwargs.ts
+++ b/src/array/_globals.kwargs.ts
@@ -1,12 +1,12 @@
 //@ts-check
 import NDArray from "../NDArray";
-import { Arr, ArrJS, kwargs_decorator, frequently_used_parsers } from "../array/kwargs";
+import { Arr, ArrOrAny, kwargs_decorator, frequently_used_parsers } from "../array/kwargs";
 
 
 export namespace Func_clip {
   export type Implementation = (a: Arr, a_min: Arr, a_max: Arr, out: Arr) => Arr;
-  export type Kwargs = { a?: ArrJS, a_min?: ArrJS, a_max?: ArrJS, out?: NDArray<any> | null };
-  export type Wrapper = (a: ArrJS, a_min: ArrJS, a_max: ArrJS, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, a_min?: ArrOrAny, a_max?: ArrOrAny, out?: NDArray<any> | null };
+  export type Wrapper = (a: ArrOrAny, a_min: ArrOrAny, a_max: ArrOrAny, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [
     ["a", undefined],
diff --git a/src/array/_globals.ts b/src/array/_globals.ts
index a590b78..f0e29dd 100644
--- a/src/array/_globals.ts
+++ b/src/array/_globals.ts
@@ -5,6 +5,17 @@ import { TypedArrayConstructor } from "../dtypes";
 const { np, NDArray: __NDArray } = GLOBALS;
 if (!__NDArray) throw new Error(`Programming error: NDArray not defined`);
 
+
+// Types used everywhere
+export type Arr = NDArray<any>;
+export type ArrOrAny = NDArray<any> | number | boolean | any[];
+export type ArrOrConst = NDArray<any> | number | boolean;
+
+export type AxisArg = number | null;
+export type ArrOrNull = NDArray<any> | null;
+
+
+
 // Functions to avoid importing NDArray (because if I import NDArray, I can't use it as a type annotation in the same file)
 export const _NDArray = __NDArray;
 
diff --git a/src/array/kwargs.ts b/src/array/kwargs.ts
index e4a763b..bce95ae 100644
--- a/src/array/kwargs.ts
+++ b/src/array/kwargs.ts
@@ -1,16 +1,12 @@
 import { NDArray } from '../NDArray';
 import { TypedArrayConstructor } from '../dtypes';
-import { asarray, broadcast_n_shapes, isarray } from "./_globals";
+import { broadcast_n_shapes, isarray } from "./_globals";
+import { asarray, Arr, ArrOrAny, ArrOrConst, AxisArg, ArrOrNull, } from './_globals';
+export { asarray, Arr, ArrOrAny, ArrOrConst, AxisArg, ArrOrNull } from './_globals';
 import { new_from } from './basic';
 import { Where } from './indexes';
 
 
-export type Arr = NDArray<any>;
-export type ArrJS = NDArray<any> | number | boolean;
-export type ArrNum = NDArray<any> | number | boolean;
-
-export type AxisArg = number | null;
-export type OutArg = NDArray<any> | null;
 
 
 
@@ -22,8 +18,8 @@ export type OutArg = NDArray<any> | null;
 
 
 export function kwargs_decorator<
-  Wrapper extends (...args: any[]) => ArrJS,
-  Implementation extends (...args: any[]) => NDArray,
+  Wrapper extends (...args: any[]) => ArrOrConst,
+  Implementation extends (...args: any[]) => Arr,
 >({ defaults, implementation, parsers, this_as_first_arg }: {
   defaults: [string, any][],
   implementation: Implementation,
@@ -165,29 +161,11 @@ export const frequently_used_parsers = {
   }
 }
 
-type NDArray_non_0D = NDArray<any> | number[];
-
-// Used in statistics.ts:
-export namespace Func_a_q_axis {
-  export type Implementation = (a: NDArray<any>, q: number, axis: number) => NDArray<any>;
-  export type Kwargs = { a?: NDArray_non_0D, q?: number, axis?: AxisArg };
-  export type Wrapper = (a: NDArray_non_0D | Kwargs, q: number | Kwargs, axis?: AxisArg | Kwargs) => NDArray<any>;
-  export const decorator = kwargs_decorator<Wrapper, Implementation>;
-  export const defaults: [string, any][] = [["a", undefined], ["q", undefined], ["axis", null]];
-  export const parsers = [
-    frequently_used_parsers.a_axis_flatten,
-    (kwargs) => { kwargs.q = asarray(kwargs.q); },
-  ];
-  export const defaultDecorator = (implementation: Implementation) => decorator({
-    defaults, implementation, parsers
-  })
-}
-
 // Used by sort:
 export namespace Func_a_lastAxis {
-  export type Implementation = (a: ArrJS, axis: number) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, axis?: AxisArg };
-  export type Wrapper = (a: ArrJS | Kwargs, axis?: AxisArg | Kwargs) => NDArray<any>;
+  export type Implementation = (a: ArrOrAny, axis: number) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, axis?: AxisArg };
+  export type Wrapper = (a: ArrOrAny | Kwargs, axis?: AxisArg | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["axis", -1]];
   export const parsers = [
@@ -202,8 +180,8 @@ export namespace Func_a_lastAxis {
 // Used by elementwise operators and methods:
 export namespace Func_a_out {
   export type Implementation = (a: NDArray<any>, out: NDArray<any> | null) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, out?: NDArray<any> | null };
-  export type Wrapper = (a: ArrJS | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, out?: NDArray<any> | null };
+  export type Wrapper = (a: ArrOrAny | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["out", null]];
   export const parsers = [
@@ -236,8 +214,8 @@ export namespace Method_out {
 // Used by round
 export namespace Func_a_decimals_out {
   export type Implementation = (a: NDArray<any>, decimals: number, out: NDArray<any> | null) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, decimals?: number, out?: NDArray<any> | null };
-  export type Wrapper = (a: ArrJS | Kwargs, decimals: number | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, decimals?: number, out?: NDArray<any> | null };
+  export type Wrapper = (a: ArrOrAny | Kwargs, decimals: number | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["decimals", 0], ["out", null]];
   export const parsers = [
@@ -270,8 +248,8 @@ export namespace Method_a_decimals_out {
 // Used by binary operators and methods:
 export namespace Func_a_other_out {
   export type Implementation = (a: NDArray<any>, other: NDArray<any>, out: NDArray<any> | null) => NDArray<any>;
-  export type Kwargs = { a: ArrJS, other?: ArrJS, out?: NDArray<any> | null };
-  export type Wrapper = (a: ArrJS | Kwargs, other: ArrJS | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
+  export type Kwargs = { a: ArrOrAny, other?: ArrOrAny, out?: NDArray<any> | null };
+  export type Wrapper = (a: ArrOrAny | Kwargs, other: ArrOrAny | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["other", undefined], ["out", null]];
   export const parsers = [
@@ -285,10 +263,10 @@ export namespace Func_a_other_out {
 }
 export namespace Method_other_out {
   export type Implementation = Func_a_other_out.Implementation;
-  export type Kwargs = { other?: ArrJS, out?: NDArray<any> | null };
+  export type Kwargs = { other?: ArrOrAny, out?: NDArray<any> | null };
   export type Wrapper<
     T extends TypedArrayConstructor = Float64ArrayConstructor,
-  > = (other: ArrJS | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<T>;
+  > = (other: ArrOrAny | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<T>;
   export const defaults: [string, any][] = [["a", undefined], ["other", undefined], ["out", null]];
   export const parsers = [
     frequently_used_parsers.isarray('a'),
@@ -307,8 +285,8 @@ export namespace Method_other_out {
 // Used by assign operators and methods:
 export namespace Func_a_values_where {
   export type Implementation = (a: NDArray<any>, values: NDArray<any>, where: Where) => NDArray<any>;
-  export type Kwargs = { a: NDArray<any>, values?: ArrJS, where: Where };
-  export type Wrapper = (a: NDArray<any> | Kwargs, values: ArrJS | Kwargs, ...where: Where) => NDArray<any>;
+  export type Kwargs = { a: NDArray<any>, values?: ArrOrAny, where: Where };
+  export type Wrapper = (a: NDArray<any> | Kwargs, values: ArrOrAny | Kwargs, ...where: Where) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["values", undefined], ["...where", null]];
   export const parsers = [
@@ -321,8 +299,8 @@ export namespace Func_a_values_where {
 }
 export namespace Method_values_where {
   export type Implementation = Func_a_values_where.Implementation;
-  export type Kwargs = { values?: ArrJS, where: Where };
-  export type Wrapper = (values: ArrJS | Kwargs, ...where: Where) => NDArray<any>;
+  export type Kwargs = { values?: ArrOrAny, where: Where };
+  export type Wrapper = (values: ArrOrAny | Kwargs, ...where: Where) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["values", undefined], ["...where", null]];
   export const parsers = [
@@ -337,10 +315,10 @@ export namespace Method_values_where {
 // Used by reduce functions and methods:
 export namespace Func_a_axis_keepdims {
   export type Implementation = (a: NDArray<any>, axis: number, keepdims: boolean) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, axis?: AxisArg, keepdims?: boolean };
+  export type Kwargs = { a?: ArrOrAny, axis?: AxisArg, keepdims?: boolean };
   export type Wrapper<
     T extends TypedArrayConstructor = Float64ArrayConstructor,
-  > = (a: ArrJS | Kwargs, axis?: AxisArg | Kwargs, keepdims?: boolean | Kwargs) => NDArray<T>;
+  > = (a: ArrOrAny | Kwargs, axis?: AxisArg | Kwargs, keepdims?: boolean | Kwargs) => NDArray<T>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["axis", null], ["keepdims", false]];
   export const parsers = [
@@ -354,6 +332,7 @@ export namespace Func_a_axis_keepdims {
     });
   }
 }
+
 export namespace Method_a_axis_keepdims {
   export type Implementation = Func_a_axis_keepdims.Implementation;
   export type Kwargs = { axis?: AxisArg, keepdims?: boolean };
@@ -378,8 +357,8 @@ export namespace Method_a_axis_keepdims {
 // For norm:
 export namespace Func_a_ord_axis_keepdims {
   export type Implementation = (a: NDArray<any>, ord: number, axis: number, keepdims: boolean) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, ord?: number, axis?: AxisArg, keepdims?: boolean };
-  export type Wrapper = (a: ArrJS | Kwargs, ord?: number | Kwargs, axis?: AxisArg | Kwargs, keepdims?: boolean | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, ord?: number, axis?: AxisArg, keepdims?: boolean };
+  export type Wrapper = (a: ArrOrAny | Kwargs, ord?: number | Kwargs, axis?: AxisArg | Kwargs, keepdims?: boolean | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["ord", 2], ["axis", null], ["keepdims", false]];
   export const parsers = [
@@ -408,8 +387,8 @@ export namespace Method_a_ord_axis_keepdims {
 // For std:
 export namespace Func_a_axis_ddof_keepdims {
   export type Implementation = (a: NDArray<any>, axis: number, ddof: number, keepdims: boolean) => NDArray<any>;
-  export type Kwargs = { a?: ArrJS, axis?: AxisArg, ddof?: number, keepdims?: boolean };
-  export type Wrapper = (a: ArrJS | Kwargs, axis?: AxisArg | Kwargs, ddof?: number | Kwargs, keepdims?: boolean | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, axis?: AxisArg, ddof?: number, keepdims?: boolean };
+  export type Wrapper = (a: ArrOrAny | Kwargs, axis?: AxisArg | Kwargs, ddof?: number | Kwargs, keepdims?: boolean | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["a", undefined], ["ddof", 0], ["axis", null], ["keepdims", false]];
   export const parsers = [
@@ -439,9 +418,9 @@ export namespace Method_a_axis_ddof_keepdims {
 
 // Used by atan2:
 export namespace Func_y_x_out {
-  export type Implementation = (y: ArrJS, x: ArrJS, out?: NDArray | null) => NDArray<any>;
-  export type Kwargs = { y?: ArrJS, x?: ArrJS, out?: NDArray<any> | null };
-  export type Wrapper = (y: ArrJS | Kwargs, x: ArrJS | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
+  export type Implementation = (y: ArrOrAny, x: ArrOrAny, out?: NDArray | null) => NDArray<any>;
+  export type Kwargs = { y?: ArrOrAny, x?: ArrOrAny, out?: NDArray<any> | null };
+  export type Wrapper = (y: ArrOrAny | Kwargs, x: ArrOrAny | Kwargs, out?: NDArray<any> | null | Kwargs) => NDArray<any>;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [["x", undefined], ["y", undefined], ["out", null]];
   export const parsers = [
diff --git a/src/index.ts b/src/index.ts
index 33f95a1..063b166 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -250,8 +250,8 @@ np.geomspace = np.modules.constructors.geomspace;
 np.take = np.modules.indexing.take;
 np.where = np.modules.indexing.where;
 np.nonzero = np.modules.indexing.nonzero;
-np.quantile = np.modules.statistics.kw_exported.quantile;
-np.nanquantile = np.modules.statistics.kw_exported.nanquantile;
+np.quantile = np.modules.statistics.quantile;
+np.nanquantile = np.modules.statistics.nanquantile;
 // np.percentile = np.modules.statistics.kw_exported.percentile;
 // np.median = np.modules.statistics.kw_exported.median;
 // np.average = np.modules.statistics.kw_exported.average;
diff --git a/src/modules/math-functions.kwargs.ts b/src/modules/math-functions.kwargs.ts
index e0a0a19..77070c2 100644
--- a/src/modules/math-functions.kwargs.ts
+++ b/src/modules/math-functions.kwargs.ts
@@ -1,12 +1,11 @@
 //@ts-check
-import NDArray from "../NDArray";
-import { Arr, ArrJS, OutArg, kwargs_decorator, frequently_used_parsers } from "../array/kwargs";
-
+import { Arr, ArrOrAny, ArrOrNull, kwargs_decorator, frequently_used_parsers, ArrOrConst } from "../array/kwargs";
+import * as math from "./math-functions";
 
 export namespace Func_clip {
   export type Implementation = (a: Arr, a_min: Arr, a_max: Arr, out: Arr) => Arr;
-  export type Kwargs = { a?: ArrJS, a_min?: ArrJS, a_max?: ArrJS, out?: OutArg };
-  export type Wrapper = (a: ArrJS, a_min: ArrJS, a_max: ArrJS, out?: OutArg | Kwargs) => NDArray<any>;
+  export type Kwargs = { a?: ArrOrAny, a_min?: ArrOrAny, a_max?: ArrOrAny, out?: ArrOrNull };
+  export type Wrapper = (a: ArrOrAny, a_min: ArrOrAny, a_max: ArrOrAny, out?: ArrOrNull | Kwargs) => ArrOrConst;
   export const decorator = kwargs_decorator<Wrapper, Implementation>;
   export const defaults: [string, any][] = [
     ["a", undefined],
@@ -23,4 +22,6 @@ export namespace Func_clip {
   export const defaultDecorator = (implementation: Implementation) => decorator({
     defaults, implementation, parsers
   });
-}
\ No newline at end of file
+}
+
+export const clip = Func_clip.defaultDecorator(math.clip);
diff --git a/src/modules/math-functions.ts b/src/modules/math-functions.ts
index 288f7aa..3651f7c 100644
--- a/src/modules/math-functions.ts
+++ b/src/modules/math-functions.ts
@@ -1,6 +1,5 @@
 //@ts-check
-import { Arr } from "../array/kwargs";
-import { Func_clip } from "./math-functions.kwargs";
+import { Arr } from "../array/_globals";
 import { op_binary } from "../array/operators";
 import { apply_along_axis, concatenate } from "../array/transform";
 
@@ -76,9 +75,37 @@ export function cumtrapz(y: Arr, x: Arr | null, dx: number, axis: number, initia
   }, y.dtype);
 }
 
+export function trapz(y: Arr, x: Arr | null, dx: number, axis: number) {
+  return apply_along_axis(y, axis, (arr: any[]) => {
+    let sum = 0;
+    for (let i = 1; i < arr.length; i++) {
+      const dxi = x ? x[i] - x[i - 1] : dx;
+      sum += (arr[i] + arr[i - 1]) / 2 * dxi;
+    }
+    return sum;
+  }, y.dtype);
+}
 
+export function interp(x: Arr, xp: Arr, fp: Arr, left: number | null, right: number | null) {
+  const x_flat = x.flat;
+  const xp_flat: number[] = xp.flat;
+  const fp_flat = fp.flat;
+  // Use binary search instead of findIndex:
+  return x_flat.map((v) => {
+    let lo = 0, hi = xp_flat.length;
+    while (lo < hi) {
+      const mid = lo + hi >> 1;
+      if (xp_flat[mid] < v) lo = mid + 1;
+      else hi = mid;
+    }
+    if (lo == 0) return left;
+    if (lo == xp_flat.length) return right;
+    const x0 = xp_flat[lo - 1];
+    const x1 = xp_flat[lo];
+    const f0 = fp_flat[lo - 1];
+    const f1 = fp_flat[lo];
+    return f0 + (f1 - f0) * (v - x0) / (x1 - x0);
+  });
+}
 
 
-export const kw_exported = {
-  clip: Func_clip.defaultDecorator(clip),
-};
\ No newline at end of file
diff --git a/src/modules/statistics.kwargs.ts b/src/modules/statistics.kwargs.ts
new file mode 100644
index 0000000..9703596
--- /dev/null
+++ b/src/modules/statistics.kwargs.ts
@@ -0,0 +1,29 @@
+//@ts-check
+import { Arr, AxisArg, kwargs_decorator, frequently_used_parsers, asarray, ArrOrConst, Func_a_axis_keepdims } from "../array/kwargs";
+import * as stats from "./statistics";
+
+
+type NDArray_non_0D = Arr | number[];
+
+// Used in statistics.ts:
+export namespace Func_a_q_axis {
+  export type Implementation = (a: Arr, q: number, axis: number) => Arr;
+  export type Kwargs = { a?: NDArray_non_0D, q?: number, axis?: AxisArg };
+  export type Wrapper = (a: NDArray_non_0D | Kwargs, q: number | Kwargs, axis?: AxisArg | Kwargs) => ArrOrConst;
+  export const decorator = kwargs_decorator<Wrapper, Implementation>;
+  export const defaults: [string, any][] = [["a", undefined], ["q", undefined], ["axis", null]];
+  export const parsers = [
+    frequently_used_parsers.a_axis_flatten,
+    (kwargs) => { kwargs.q = asarray(kwargs.q); },
+  ];
+  export const defaultDecorator = (implementation: Implementation) => decorator({
+    defaults, implementation, parsers
+  })
+}
+
+export const quantile = Func_a_q_axis.defaultDecorator(stats.quantile);
+export const nanquantile = Func_a_q_axis.defaultDecorator(stats.nanquantile);
+export const percentile = Func_a_q_axis.defaultDecorator(stats.percentile);
+export const nanpercentile = Func_a_q_axis.defaultDecorator(stats.nanpercentile);
+export const median = Func_a_axis_keepdims.defaultDecorator(stats.median);
+export const nanmedian = Func_a_axis_keepdims.defaultDecorator(stats.nanmedian);
\ No newline at end of file
diff --git a/src/modules/statistics.ts b/src/modules/statistics.ts
index ada7e7d..be65109 100644
--- a/src/modules/statistics.ts
+++ b/src/modules/statistics.ts
@@ -1,43 +1,121 @@
 //@ts-check
+import assert = require("assert");
 import NDArray from "../NDArray";
-import { asarray, new_NDArray } from "../array/_globals";
-import { AxisArg, Func_a_q_axis, Func_a_axis_keepdims } from "../array/kwargs";
-import { n_ary_operation, op_binary } from "../array/operators";
+import { ArrOrAny, asarray, new_NDArray } from "../array/_globals";
+import { as_number } from "../array/basic";
+import { Arr, AxisArg } from "../array/kwargs";
+
 import { cmp_nan_at_the_end, swapAxes } from "../array/transform";
 import { np } from "./_globals";
 
 
-/**
- * Compute the q-th percentile of the data along the specified axis.
- */
-export function percentile(a: NDArray<any>, q: NDArray<any> | number, axis: AxisArg | null) {
-  q = np.divide(q, 100);
-  return quantile(a, q, axis);
+
+
+export const sum_list = (arr: number[] | boolean[]) => arr.reduce((a, b) => (a as number) + (b as number), 0) as number;
+export const product_list = (arr: number[] | boolean[]) => arr.reduce((a, b) => (a as number) + (b as number), 1) as number;
+export function argmax_list(arr: number[] | boolean[]) {
+  if (arr.length == 0) throw new Error("argmax of an empty array");
+  let i = 0, v = arr[0];
+  for (let j = 1; j < arr.length; j++) if (arr[j] > v) v = arr[i = j];
+  return i;
+}
+export function argmin_list(arr: number[] | boolean[]) {
+  if (arr.length == 0) throw new Error("argmin of an empty array");
+  let i = 0, v = arr[0];
+  for (let j = 1; j < arr.length; j++) if (arr[j] < v) v = arr[i = j];
+  return i;
+}
+export function max_list<T>(arr: T[]): T {
+  if (arr.length == 0) throw new Error("max of an empty array");
+  if (typeof arr[0] == "boolean") return any_list(arr as boolean[]) as T;
+  return Math.max(...arr as number[]) as T;
+}
+export function min_list<T>(arr: T[]): T {
+  if (arr.length == 0) throw new Error("min of an empty array");
+  if (typeof arr[0] == "boolean") return all_list(arr as boolean[]) as T;
+  return Math.min(...arr as number[]) as T;
+}
+export const any_list = (arr: number[] | boolean[]) => {
+  for (let x of arr) if (x) return true;
+  return false;
+};
+export const all_list = (arr: number[] | boolean[]) => {
+  for (let x of arr) if (!x) return false;
+  return true;
+}
+export const mean_list = (arr: number[] | boolean[]) => sum_list(arr) / arr.length;
+export const var_list = (arr: number[] | boolean[], ddof: number = 0) => {
+  const mean = mean_list(arr);
+  const n = arr.length;
+  return sum_list(arr.map(x => (x - mean) ** 2)) / (n - ddof);
+}
+export const std_list = (arr: number[] | boolean[], ddof: number = 0) => {
+  return Math.pow(var_list(arr, ddof), 0.5);
 }
 
 
-function js_quantile_a_1D_q_01D(
-  a_flat: number[],
-  q_flat: number[],
+/**
+ * a must be a list.
+ * q can be a number in [0,1] or a list.
+ * the output is shaped as q.
+ * @param a
+ * @param q
+ * @param _nan_handling 
+ * @returns 
+ */
+function quantile_list<T extends number | number[] = number>(
+  a: number[] | boolean[],
+  q: T,
   _nan_handling = false,
-): number[] {
-  a_flat = [...a_flat].sort(cmp_nan_at_the_end);
-  let n = a_flat.length;
-  if (_nan_handling) n -= a_flat.reduce((cum, x) => cum + (isNaN(x) ? 1 : 0), 0);
-  return q_flat.map(q => {
+): T {
+  const x = a.map(v => +v).sort(cmp_nan_at_the_end);
+  let n = x.length;
+  if (_nan_handling) n -= x.reduce((cum, x) => cum + (isNaN(x) ? 1 : 0), 0);
+  const f = (q: number) => {
     const nq = q * (n - 1);
     const lo = Math.floor(nq);
-    if (nq == lo) return a_flat[lo];
+    if (nq == lo) return x[lo];
     const hi = lo + 1;
-    if (nq == hi) return a_flat[hi];
-    return a_flat[lo] * (hi - nq) + a_flat[hi] * (nq - lo);
-  });
+    if (nq == hi) return x[hi];
+    return x[lo] * (hi - nq) + x[hi] * (nq - lo);
+  }
+  return (Array.isArray(q) ? q.map(f) : f(q)) as T;
+}
+
+// const norm_list = (arr: number[], ord: number = 2) => {
+//   if (ord % 2 != 0) arr = abs_list(arr);
+//   if (ord == Infinity) return max_list(arr);
+//   if (ord == 1) return sum_list(arr);
+//   return Math.pow(sum_list(pow_list(arr, ord)), 1 / ord);
+// }
+
+
+// MOVE THIS!!!
+export function linspace_list(start: number, stop: number, steps = 50, endpoint = true) {
+  let n = (steps - (endpoint ? 1 : 0));
+  const shift = (stop - start) / n;
+  let arr = new Array(steps);
+  for (let i = 0; i < steps; i++) arr[i] = start + i * shift;
+  return arr;
+}
+export const median_list = (a: number[] | boolean[]) => quantile_list(a, 0.5);
+
+
+
+
+
+/**
+ * Compute the q-th percentile of the data along the specified axis.
+ */
+export function percentile(a: Arr, q: Arr | number, axis: AxisArg | null) {
+  q = np.divide(q, 100);
+  return quantile(a, q, axis);
 }
 
 /**
  * Compute the q-th quantile of the data along the specified axis.
  */
-export function quantile(a: NDArray<any>, q: NDArray<any> | number, axis: number, _nan_handling = false) {
+export function quantile(a: Arr, q: Arr | number, axis: number, _nan_handling = false): Arr {
   q = asarray(q);
   if (axis != a.shape.length - 1) { a = swapAxes(a, axis, -1); }
   const outer_shape = q.shape;
@@ -50,7 +128,7 @@ export function quantile(a: NDArray<any>, q: NDArray<any> | number, axis: number
   const out = new Float64Array(q.size * nrows);
   for (let i = 0; i < nrows; i++) {
     const row = a_flat.slice(i * ncols, (i + 1) * ncols);
-    const values = js_quantile_a_1D_q_01D(row, q_flat, _nan_handling);
+    const values = quantile_list(row, q_flat, _nan_handling);
     let j = i;
     for (let k in values) {
       out[j] = values[k];
@@ -75,7 +153,7 @@ export function quantile(a: NDArray<any>, q: NDArray<any> | number, axis: number
 /**
  * Compute the median along the specified axis.
  */
-export function median(a: NDArray, axis: AxisArg | null, keepdims: boolean) {
+export function median(a: Arr, axis: AxisArg | null, keepdims: boolean) {
   let out = quantile(a, 0.5, axis);
   if (keepdims) {
     const shape = a.shape.slice();
@@ -88,7 +166,7 @@ export function median(a: NDArray, axis: AxisArg | null, keepdims: boolean) {
 /**
  * Compute the weighted average along the specified axis.
  */
-export function average(a: NDArray, axis: AxisArg | null, weights: NDArray | null, keepdims: boolean) {
+export function average(a: Arr, axis: AxisArg | null, weights: NDArray | null, keepdims: boolean) {
   if (weights === null) return np.mean(a, axis);
   // MISSING: assert weights is 1D
   const denominator = np.sum(weights);
@@ -103,7 +181,7 @@ export function average(a: NDArray, axis: AxisArg | null, weights: NDArray | nul
 /**
  * Compute the q-th percentile of the data along the specified axis, while ignoring nan values.
 */
-export function nanpercentile(a: NDArray, q: NDArray | number, axis: AxisArg | null) {
+export function nanpercentile(a: Arr, q: Arr | number, axis: AxisArg | null) {
   q = np.divide(q, 100);
   return nanquantile(a, q, axis);
 }
@@ -112,7 +190,7 @@ export function nanpercentile(a: NDArray, q: NDArray | number, axis: AxisArg | n
 /**
  * Compute the q-th quantile of the data along the specified axis, while ignoring nan values.
  */
-export function nanquantile(a: NDArray, q: NDArray | number, axis: AxisArg | null) {
+export function nanquantile(a: Arr, q: Arr | number, axis: AxisArg | null) {
   return quantile(a, q, axis, true)
 }
 
@@ -120,7 +198,7 @@ export function nanquantile(a: NDArray, q: NDArray | number, axis: AxisArg | nul
 /**
  * Compute the median along the specified axis, while ignoring NaNs.
  */
-export function nanmedian(a: NDArray, axis: number, keepdims: boolean) {
+export function nanmedian(a: Arr, axis: number, keepdims: boolean) {
   let out = nanquantile(a, 0.5, axis);
   if (keepdims) {
     const shape = a.shape.slice();
@@ -137,7 +215,7 @@ export function nanmedian(a: NDArray, axis: number, keepdims: boolean) {
 // /**
 //  * Return Pearson product-moment correlation coefficients.
 //  */
-// export function corrcoef(x: NDArray, y: NDArray, rowvar = true, bias = false, ddof = 0, dtype = Number) {
+// export function corrcoef(x: Arr, y: Arr, rowvar = true, bias = false, ddof = 0, dtype = Number) {
 //   const x_shape = x.shape;
 //   const y_shape = y.shape;
 //   if (x_shape.length !== 1 || y_shape.length !== 1) throw new Error("Expected 1D arrays");
@@ -155,14 +233,14 @@ export function nanmedian(a: NDArray, axis: number, keepdims: boolean) {
 // /**
 //  * Cross-correlation of two 1-dimensional sequences.
 //  */
-// export function correlate(a: NDArray, v: NDArray, mode) {
+// export function correlate(a: Arr, v: Arr, mode) {
 //   throw new Error("Not implemented");
 // }
 
 // /**
 //  * Estimate a covariance matrix, given data and weights.
 //  */
-// export function cov(m: NDArray, y: NDArray, rowvar = true, bias = false, ddof = 0, fweights = null, aweights = null) {
+// export function cov(m: Arr, y: Arr, rowvar = true, bias = false, ddof = 0, fweights = null, aweights = null) {
 //   const m_shape = m.shape;
 //   const y_shape = y.shape;
 //   if (m_shape.length !== 2 || y_shape.length !== 2) throw new Error("Expected 2D arrays");
@@ -185,44 +263,101 @@ export function nanmedian(a: NDArray, axis: number, keepdims: boolean) {
 // /**
 //  * Compute the histogram of a dataset.
 //  */
-// export function histogram(a: NDArray, bins, range, density, weights) {
+// export function histogram(a: Arr, bins, range, density, weights) {
 //   throw new Error("Not implemented");
 // }
 
 // /**
 //  * Compute the bi-dimensional histogram of two data samples.
 //  */
-// export function histogram2d(x: NDArray, y: NDArray, bins, range, density, weights) {
+// export function histogram2d(x: Arr, y: Arr, bins, range, density, weights) {
 //   throw new Error("Not implemented");
 // }
 
 // /**
 //  * Compute the multidimensional histogram of some data.
 //  */
-// export function histogramdd(sample: NDArray, bins, range, density, weights) {
+// export function histogramdd(sample: Arr, bins, range, density, weights) {
 //   throw new Error("Not implemented");
 // }
 
 // /**
 //  * Count number of occurrences of each value in array of non-negative ints.
 //  */
-// export function bincount(x: NDArray, weights, minlength) {
+// export function bincount(x: Arr, weights, minlength) {
 //   throw new Error("Not implemented");
 // }
 
-// /**
-//  * Function to calculate only the edges of the bins used by the histogram function.
-//  */
-// export function histogram_bin_edges(a: NDArray, bins, range, weights) {
-// }
+/**
+ * Function to calculate only the edges of the bins used by the histogram function.
+ */
+
+export function _n_bins(a: number[], bins: number | number[] | 'auto' | 'rice' | 'fd' | 'scott' | 'sqrt' | 'sturges', range: [number, number]) {
+  const n = a.length;
+  if (bins === 'auto') {
+    return Math.max(_n_bins(a, 'fd', range), _n_bins(a, 'sturges', range));
+  }
+  if (bins == 'sturges') {
+    return Math.ceil(Math.log2(n) + 1);
+  }
+  if (bins == 'fd') {
+    const IQR = quantile_list(a, 0.75) - quantile_list(a, 0.25);
+    return Math.ceil(2 * IQR / Math.pow(n, 1 / 3));
+  }
+  if (bins == 'rice') {
+    const nh = 2 * Math.pow(n, 1 / 3);
+    return Math.ceil(n / nh);
+  }
+  if (bins == 'scott') {
+    const std = std_list(a);
+    const h = std * Math.pow(24 * Math.sqrt(Math.PI) / n, 1 / 3);
+    return Math.ceil((max_list(a) - min_list(a)) / h);
+  }
+  if (bins == 'sqrt') {
+    return Math.ceil(Math.sqrt(n));
+  }
+  throw new Error("Unsupported bins value");
+}
+
+export function histogram_bin_edges_list(a: number[], bins: number | number[] | 'auto' | 'rice' | 'fd' | 'scott' | 'sqrt' = 'auto', range: [number, number] | null = null): number[] {
+  if (range === null) range = [min_list(a), max_list(a)];
+  else {
+    // Exclude values outside the range
+    a = a.filter(x => x >= range[0] && x <= range[1]);
+  }
+  if (typeof bins === 'string') bins = _n_bins(a, bins, range);
+  if (typeof bins === 'number') bins = linspace_list(min_list(a), max_list(a), bins + 1);
+  assert(Array.isArray(bins), "bins must have at least 2 elements");
+  return bins;
+}
+
+
+export function histogram_list(a: number[], bins: number | number[] | 'auto' | 'rice' | 'fd' | 'scott' | 'sqrt', range: [number, number] | null = null, weights: number[] | null = null): [number[], number[]] {
+  const edges = histogram_bin_edges_list(a, bins, range);
+  const counts = new Array(edges.length - 1).fill(0);
+  for (let i = 0; i < a.length; i++) {
+    const x = a[i];
+    let j = 0;
+    while (j < edges.length - 1 && x >= edges[j + 1]) j++;
+    counts[j] += weights ? weights[i] : 1;
+  }
+  return [edges, counts];
+}
+
 
 // /**
 //  * Return the indices of the bins to which each value in input array belongs.
 //  */
-// export function digitize(x: NDArray, bins, right) {
-//   throw new Error("Not implemented");
-// }
-
+export function digitize_list(x: number[], bins: number[], right = false): number[] {
+  const out = new Array(x.length);
+  for (let i = 0; i < x.length; i++) {
+    const xi = x[i];
+    let j = 0;
+    while (j < bins.length && (right ? xi > bins[j] : xi >= bins[j])) j++;
+    out[i] = j;
+  }
+  return out;
+}
 
 
 
@@ -254,12 +389,3 @@ Averages and variances
   - nanvar(a[, axis, dtype, out, ddof, ...]) Compute the variance along the specified axis, while ignoring NaNs.
  */
 
-
-export const kw_exported = {
-  quantile: Func_a_q_axis.defaultDecorator(quantile),
-  nanquantile: Func_a_q_axis.defaultDecorator(nanquantile),
-  percentile: Func_a_q_axis.defaultDecorator(percentile),
-  nanpercentile: Func_a_q_axis.defaultDecorator(nanpercentile),
-  median: Func_a_axis_keepdims.defaultDecorator(median),
-  nanmedian: Func_a_axis_keepdims.defaultDecorator(nanmedian),
-};
\ No newline at end of file