Add cdf for poisson

notebooks/models/poisson.md

@@ -33,18 +33,36 @@ The characteristic exponent is given by
 
 ```{code-cell} ipython3
 from quantflow.sp.poisson import PoissonProcess
-pr = PoissonProcess(intensity=2)
+pr = PoissonProcess(intensity=1)
 pr
 ```
 
+```{code-cell} ipython3
+m = pr.marginal(0.1)
+import numpy as np
+cdf = m.cdf(np.arange(5))
+cdf
+```
+
+```{code-cell} ipython3
+n = 128*8
+m.cdf_from_characteristic(5, frequency_n=n).y
+```
+
+```{code-cell} ipython3
+cdf1 = m.cdf_from_characteristic(5, frequency_n=n).y
+cdf2 = m.cdf_from_characteristic(5, frequency_n=n, simpson_rule=False).y
+10000*np.max(np.abs(cdf-cdf1)), 10000*np.max(np.abs(cdf-cdf2))
+```
+
 ### Marginal
 
 ```{code-cell} ipython3
 import numpy as np
 from quantflow.utils import plot
 
 m = pr.marginal(1)
-plot.plot_marginal_pdf(m, 16)
+plot.plot_marginal_pdf(m, frequency_n=128*8)
 ```
 
 ```{code-cell} ipython3
@@ -206,8 +224,8 @@ The intensity function of a DSPP is given by:
 
 ```{code-cell} ipython3
 from quantflow.sp.dsp import DSP, PoissonProcess, CIR
-pr = DSP(intensity=CIR(sigma=1, kappa=1), poisson=PoissonProcess(intensity=3))
-pr2 = DSP(intensity=CIR(sigma=0.1, kappa=10), poisson=PoissonProcess(intensity=3))
+pr = DSP(intensity=CIR(sigma=1, kappa=1), poisson=PoissonProcess(intensity=1))
+pr2 = DSP(intensity=CIR(sigma=0.1, kappa=10), poisson=PoissonProcess(intensity=1))
 pr
 ```
 
@@ -216,10 +234,11 @@ import numpy as np
 from quantflow.utils import plot
 import plotly.graph_objects as go
 
-m = pr.marginal(1)
-pdf = m.pdf_from_characteristic(16)
-fig = plot.plot_marginal_pdf(m, 16, analytical=False, label=f"sigma={pr.intensity.sigma}")
-plot.plot_marginal_pdf(pr2.marginal(1), 16, analytical=False, fig=fig, marker_color="yellow", label=f"sigma={pr2.intensity.sigma}")
+n=16
+m = pr.marginal(2)
+pdf = m.pdf_from_characteristic(n)
+fig = plot.plot_marginal_pdf(m, n, analytical=False, label=f"sigma={pr.intensity.sigma}")
+plot.plot_marginal_pdf(pr2.marginal(1), n, analytical=False, fig=fig, marker_color="yellow", label=f"sigma={pr2.intensity.sigma}")
 fig.add_trace(go.Scatter(x=pdf.x, y=pr.poisson.marginal(1).pdf(pdf.x), name="Poisson", mode="markers", marker_color="blue"))
 ```
 
@@ -233,7 +252,7 @@ pr2.marginal(1).mean(), pr2.marginal(1).variance()
 
 ```{code-cell} ipython3
 from quantflow.utils.plot import plot_characteristic
-m = pr.marginal(1)
+m = pr.marginal(2)
 plot_characteristic(m)
 ```
 

pyproject.toml

@@ -24,7 +24,7 @@ black = "^23.3.0"
 pytest-cov = "^4.0.0"
 mypy = "^1.4.0"
 ghp-import = "^2.0.2"
-ruff = "^0.0.277"
+ruff = "^0.0.285"
 pytest-asyncio = "^0.21.1"
 
 

quantflow/sp/dsp.py

@@ -5,8 +5,9 @@
 from scipy.optimize import Bounds
 
 from ..utils.types import FloatArray, FloatArrayLike, Vector
+from .base import StochasticProcess1DMarginal
 from .cir import CIR, IntensityProcess
-from .poisson import PoissonBase, PoissonProcess, poisson_arrivals
+from .poisson import MarginalDiscrete1D, PoissonBase, PoissonProcess, poisson_arrivals
 
 
 class DSP(PoissonBase):
@@ -23,6 +24,9 @@ class DSP(PoissonBase):
     )
     poisson: PoissonProcess = Field(default_factory=PoissonProcess, exclude=True)
 
+    def marginal(self, t: FloatArrayLike) -> StochasticProcess1DMarginal:
+        return MarginalDiscrete1D(process=self, t=t)
+
     def frequency_range(self, std: float, max_frequency: float | None = None) -> Bounds:
         """Frequency range of the process"""
         return Bounds(0, np.pi)

quantflow/sp/poisson.py

@@ -1,16 +1,18 @@
 from abc import abstractmethod
-from typing import Generic, TypeVar
+from typing import Any, Generic, TypeVar
 
 import numpy as np
 from pydantic import Field
+from scipy.integrate import simpson
 from scipy.optimize import Bounds
 from scipy.stats import poisson
 
 from ..utils.distributions import Distribution1D
 from ..utils.functions import factorial
 from ..utils.paths import Paths
+from ..utils.transforms import TransformResult
 from ..utils.types import FloatArray, FloatArrayLike, Vector
-from .base import Im, StochasticProcess1D
+from .base import Im, StochasticProcess1D, StochasticProcess1DMarginal
 
 D = TypeVar("D", bound=Distribution1D)
 
@@ -63,6 +65,9 @@ def poisson_arrivals(intensity: float, time_horizon: float = 1) -> list[float]:
 class PoissonProcess(PoissonBase):
     intensity: float = Field(default=1.0, ge=0, description="intensity rate")
 
+    def marginal(self, t: FloatArrayLike) -> StochasticProcess1DMarginal:
+        return MarginalDiscrete1D(process=self, t=t)
+
     def characteristic_exponent(self, t: Vector, u: Vector) -> Vector:
         return t * self.intensity * (1 - np.exp(Im * u))
 
@@ -159,3 +164,44 @@ def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
     def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
         """Expected variance at a time horizon"""
         return self.intensity * t * (self.jumps.variance() + self.jumps.mean() ** 2)
+
+
+class MarginalDiscrete1D(StochasticProcess1DMarginal):
+    def pdf_from_characteristic(
+        self,
+        n: int | None = None,
+        **kwargs: Any,
+    ) -> TransformResult:
+        cdf = self.cdf_from_characteristic(n, **kwargs)
+        return cdf._replace(y=np.diff(cdf.y, prepend=0))
+
+    def cdf_from_characteristic(
+        self,
+        n: int | None = None,
+        *,
+        frequency_n: int | None = None,
+        simpson_rule: bool = True,
+        **kwargs: Any,
+    ) -> TransformResult:
+        n = n or 10
+        transform = self.get_transform(frequency_n, self.support, **kwargs)
+        frequency = transform.frequency_domain
+        c = self.characteristic(frequency)
+        a = 1 / np.pi
+        result = []
+        x = np.arange(n or 10)
+        for m in x:
+            d = np.sin(0.5 * frequency)
+            d[0] = 1.0
+            f = (
+                np.sin(0.5 * (m + 1) * frequency)
+                * (c * np.exp(-0.5 * Im * m * frequency)).real
+                / d
+            )
+            f[0] = c[0].real  # type: ignore[index]
+            if simpson_rule:
+                result.append(a * simpson(f, frequency))
+            else:
+                result.append(a * np.trapz(f, frequency))
+        pdf = np.maximum(np.diff(result, prepend=0), 0)
+        return TransformResult(x=x, y=np.cumsum(pdf))  # type: ignore[arg-type]

quantflow/utils/marginal.py

@@ -79,6 +79,15 @@ def variance_from_characteristic(self) -> FloatArrayLike:
         s = -(c1 - 2 * c0 + c2) / (d * d) - m * m
         return s.real
 
+    def cdf(self, x: FloatArrayLike) -> FloatArrayLike:
+        """
+        Compute the cumulative distribution function
+
+        :param n: Location in the stochastic process domain space. If a numpy array,
+            the output should have the same shape as the input.
+        """
+        raise NotImplementedError("Analytical CFD not available")
+
     def pdf(self, x: FloatArrayLike) -> FloatArrayLike:
         """
         Computes the probability density (or mass) function of the process.
@@ -90,7 +99,7 @@ def pdf(self, x: FloatArrayLike) -> FloatArrayLike:
         :param n: Location in the stochastic process domain space. If a numpy array,
             the output should have the same shape as the input.
         """
-        return self.cdf(x) - self.cdf(x - 1)
+        raise NotImplementedError("Analytical PFD not available")
 
     def pdf_from_characteristic(
         self,
@@ -99,6 +108,7 @@ def pdf_from_characteristic(
         max_frequency: float | None = None,
         simpson_rule: bool = False,
         use_fft: bool = False,
+        frequency_n: int | None = None,
     ) -> TransformResult:
         """
         Compute the probability density function from the characteristic function.
@@ -112,7 +122,7 @@ def pdf_from_characteristic(
         :param use_fft: Use FFT for the transform. Default is False.
         """
         transform = self.get_transform(
-            n,
+            frequency_n or n,
             self.support,
             max_frequency=max_frequency,
             simpson_rule=simpson_rule,
@@ -121,6 +131,17 @@ def pdf_from_characteristic(
         psi = cast(np.ndarray, self.characteristic(transform.frequency_domain))
         return transform(psi, use_fft=use_fft)
 
+    def cdf_from_characteristic(
+        self,
+        n: int | None = None,
+        *,
+        max_frequency: float | None = None,
+        simpson_rule: bool = False,
+        use_fft: bool = False,
+        frequency_n: int | None = None,
+    ) -> TransformResult:
+        raise NotImplementedError("CFD not available")
+
     def call_option(
         self,
         n: int | None = None,
@@ -214,15 +235,6 @@ def get_transform(
             simpson_rule=simpson_rule,
         )
 
-    def cdf(self, x: FloatArrayLike) -> FloatArrayLike:
-        """
-        Compute the cumulative distribution function
-
-        :param n: Location in the stochastic process domain space. If a numpy array,
-            the output should have the same shape as the input.
-        """
-        raise NotImplementedError("Analytical CFD not available")
-
     def pdf_jacobian(self, x: FloatArrayLike) -> FloatArrayLike:
         """
         Jacobian of the pdf with respect to the parameters of the process.

tests/test_poisson.py

@@ -36,6 +36,16 @@ def test_characteristic(poisson: PoissonProcess) -> None:
     assert pytest.approx(m2.variance_from_characteristic(), 0.001) == 4
 
 
+def test_poisson_cdf_from_characteristic(poisson: PoissonProcess) -> None:
+    m = poisson.marginal(0.1)
+    # m.pdf(x)
+    cdf1 = m.cdf(1.0 * np.arange(10))
+    cdf2 = m.cdf_from_characteristic(10, frequency_n=128 * 8).y
+    np.testing.assert_almost_equal(cdf1, cdf2, decimal=3)
+    # TODO: fix this
+    # np.testing.assert_almost_equal(pdf, c_pdf.y[:10])
+
+
 def test_poisson_pdf(poisson: PoissonProcess) -> None:
     m = poisson.marginal(1)
     analytical_tests(poisson)
@@ -68,3 +78,10 @@ def test_dsp_sample(dsp: DSP):
     assert mean[0] == 0
     std = paths.std()
     assert std[0] == 0
+
+
+def test_dsp_pdf(dsp: DSP):
+    m = dsp.marginal(1)
+    pdf1 = m.pdf_from_characteristic(32).y
+    pdf2 = m.pdf_from_characteristic(64).y
+    np.testing.assert_almost_equal(pdf2[:32], pdf1)