From be0a32ec17d1f5b04ed8f2e52604c70c69f416b2 Mon Sep 17 00:00:00 2001
From: Spencer Lyon <spencerlyon2@gmail.com>
Date: Tue, 20 Mar 2018 18:04:30 -0400
Subject: [PATCH] RFC: Sl/qnwdist (#202)

* WIP: quadrature weights and nodes for distributions

* DOC: add docstring to qnwdist

* ENH: export qnwdist
---
 src/QuantEcon.jl |  9 +++---
 src/quad.jl      | 72 ++++++++++++++++++++++++++++++++++++++++++++----
 2 files changed, 71 insertions(+), 10 deletions(-)

diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl
index 1abfcf4b..bd5c8f5b 100644
--- a/src/QuantEcon.jl
+++ b/src/QuantEcon.jl
@@ -53,7 +53,8 @@ export
     discrete_var, Even, Quantile, Quadrature,
 
 # modeltools
-    AbstractUtility, LogUtility, CRRAUtility, CFEUtility, EllipticalUtility, derivative,
+    AbstractUtility, LogUtility, CRRAUtility, CFEUtility, EllipticalUtility,
+    derivative,
 
 # gth_solve
     gth_solve,
@@ -96,8 +97,8 @@ export
     smooth, periodogram, ar_periodogram,
 
 # util
-    meshgrid, gridmake, gridmake!, ckron, is_stable, num_compositions, simplex_grid, simplex_index,
-    next_k_array!, k_array_rank,
+    meshgrid, gridmake, gridmake!, ckron, is_stable, num_compositions,
+    simplex_grid, simplex_index, next_k_array!, k_array_rank,
 
 # robustlq
     RBLQ,
@@ -109,7 +110,7 @@ export
 
 # quad
     qnwlege, qnwcheb, qnwsimp, qnwtrap, qnwbeta, qnwgamma, qnwequi, qnwnorm,
-    qnwunif, qnwlogn, qnwmonomial1, qnwmonomial2,
+    qnwunif, qnwlogn, qnwmonomial1, qnwmonomial2, qnwdist,
     quadrect,
     do_quad,
 
diff --git a/src/quad.jl b/src/quad.jl
index a0548695..feeea4f2 100644
--- a/src/quad.jl
+++ b/src/quad.jl
@@ -431,13 +431,13 @@ function qnwgamma(n::Int, a::Real=1.0, b::Real=1.0)
         end
 
         # rootfinding iterations
-        it = 0
         pp = 0.0
         p2 = 0.0
-        for it = 1:maxit
+        err = 100.0
+        for it in 1:maxit
             p1 = 1.0
             p2 = 0.0
-            for j=1:n
+            for j in 1:n
                 # Recurrance relation for Laguerre polynomials
                 p3 = p2
                 p2 = p1
@@ -446,12 +446,13 @@ function qnwgamma(n::Int, a::Real=1.0, b::Real=1.0)
             pp = (n*p1-(n+a)*p2)./z
             z1 = z
             z = z1-p1 ./ pp
-            if abs(z - z1) < 3e-14
+            err = abs(z - z1)
+            if err < 3e-14
                 break
             end
-        end
 
-        if it >= maxit
+        end
+        if err > 3e-14
             error("failure to converge.")
         end
 
@@ -870,3 +871,62 @@ function qnwmonomial2(vcv::AbstractMatrix)
                1/(n+2)^2 * ones(size(z2, 1)))
     return ϵj, ωj
 end
+
+
+function _quadnodes(
+        d::Distributions.ContinuousUnivariateDistribution, N::Int,
+        q0::Real, qN::Real, ::Union{Even,Type{Even}}
+    )
+    collect(linspace(quantile(d, q0), quantile(d, qN), N))
+end
+
+function _quadnodes(
+        d::Distributions.ContinuousUnivariateDistribution, N::Int,
+        q0::Real, qN::Real, ::Union{Quantile,Type{Quantile}}
+    )
+    quantiles = linspace(q0, qN, N)
+    z = quantile.(d, quantiles)
+end
+
+"""
+    qnwdist(
+        d::Distributions.ContinuousUnivariateDistribution, N::Int,
+        q0::Real=0.001, qN::Real=0.999, method::Union{T,Type{T}}=Quantile
+    ) where T
+
+Construct `N` quadrature weights and nodes for distribution `d` from the
+quantile `q0` to the quantile `qN`. `method` can be one of:
+
+- `Even`: nodes will be evenly spaced between the quantiles
+- `Quantile`: nodes will be placed at evenly spaced quantile values
+
+To construct the weights, consider splitting the nodes into cells centered at
+each node. Specifically, let notation `z_i` mean the `i`th node and let
+`z_{i-1/2}` be 1/2 between nodes `z_{i-1}` and `z_i`. Then, weights are
+determined as follows:
+
+- `weights[1] = cdf(d, z_{1+1/2})`
+- `weights[N] = 1 - cdf(d, z_{N-1/2})`
+- `weights[i] = cdf(d, z_{i+1/2}) - cdf(d, z_{i-1/2})` for all i in 2:N-1
+
+In effect, this strategy assigns node `i` all the probability associated with a
+random variable occuring within the node `i`s cell.
+
+The weights always sum to 1, so they can be used as a proper probability
+distribution. This means that `E[f(x) | x ~ d] ≈ dot(f.(nodes), weights)`.
+"""
+function qnwdist(
+        d::Distributions.ContinuousUnivariateDistribution, N::Int,
+        q0::Real=0.001, qN::Real=0.999, method::Union{T,Type{T}}=Quantile
+    ) where T
+
+    z = _quadnodes(d, N, q0, qN, method)
+    zprob = zeros(N)
+
+    for i in 2:N-1
+        zprob[i] = cdf(d, (z[i] + z[i+1])/2) - cdf(d, (z[i] + z[i-1])/2)
+    end
+    zprob[1] = cdf(d, (z[1] + z[2])/2)
+    zprob[end] = 1 - cdf(d, (z[end-1] + z[end])/2)
+    return z, zprob
+end