Refactor RandomVariable and merge with Distribution

benchmarks/distribution.py

@@ -50,9 +50,9 @@ def get_randvar(distribution_name):
     elif distribution_name == "symmatrixvar_normal":
         distribution = prob.Normal(mean=mean_2d_mat, cov=cov_2d_symkron)
     elif distribution_name == "operatorvar_normal":
-        distribution = prob.Normal(mean=mean_2d_linop)
+        distribution = prob.Normal(mean=mean_2d_linop, cov=cov_2d_symkron)
 
-    return prob.RandomVariable(distribution=distribution)
+    return distribution
 
 
 class Functions:
@@ -72,13 +72,13 @@ def time_distr_functions(self, dist, property):
         """Times evaluation of the pdf, logpdf, cdf and logcdf."""
         try:
             if property == "pdf":
-                self.randvar.distribution.pdf(x=self.eval_point)
+                self.randvar.pdf(x=self.eval_point)
             elif property == "logpdf":
-                self.randvar.distribution.pdf(x=self.eval_point)
+                self.randvar.pdf(x=self.eval_point)
             elif property == "cdf":
-                self.randvar.distribution.pdf(x=self.quantile)
+                self.randvar.pdf(x=self.quantile)
             elif property == "logcdf":
-                self.randvar.distribution.pdf(x=self.quantile)
+                self.randvar.pdf(x=self.quantile)
         except NotImplementedError:
             pass
 

docs/source/introduction/quickstart.ipynb

@@ -121,7 +121,7 @@
    ],
    "source": [
     "y = 2 * x + 1\n",
-    "print(f\"Mean: {y.mean()} and covariance: {y.cov()} of y.\")"
+    "print(f\"Mean: {y.mean} and covariance: {y.cov} of y.\")"
    ]
   },
   {
@@ -164,7 +164,7 @@
     "\n",
     "# Solve using ProbNum\n",
     "x_rv, _, _, info = probnum.linalg.problinsolve(A, b)\n",
-    "print(x_rv.mean())"
+    "print(x_rv.mean)"
    ]
   },
   {
@@ -212,7 +212,7 @@
    ],
    "source": [
     "# Covariance of solution representing uncertainty\n",
-    "print(f\"Covariance matrix: \\n{x_rv.cov().todense()}\")"
+    "print(f\"Covariance matrix: \\n{x_rv.cov.todense()}\")"
    ]
   },
   {
@@ -591,7 +591,7 @@
    ],
    "source": [
     "# Plot of true solution, mean and samples\n",
-    "rvdict = {\"$x_*$\" : x, \"$\\mathbb{E}(\\mathsf{x})$\" : x_rv.mean(), \n",
+    "rvdict = {\"$x_*$\" : x, \"$\\mathbb{E}(\\mathsf{x})$\" : x_rv.mean, \n",
     "          \"$\\mathsf{x}_1$\" : x_samples[0], \"$\\mathsf{x}_2$\" : x_samples[1], \"$\\mathsf{x}_3$\" : x_samples[2]}\n",
     "vmin = np.min([np.min(mat) for mat in list(rvdict.values())])\n",
     "vmax = np.max([np.max(mat) for mat in list(rvdict.values())])\n",
@@ -731,7 +731,7 @@
    ],
    "source": [
     "# Covariance of solution representing uncertainty\n",
-    "print(f\"Covariance matrix: \\n{x_rv.cov().todense()}\")"
+    "print(f\"Covariance matrix: \\n{x_rv.cov.todense()}\")"
    ]
   },
   {
@@ -1110,7 +1110,7 @@
    ],
    "source": [
     "# Plot of true solution, mean and samples\n",
-    "rvdict = {\"$x_*$\" : x, \"$\\mathbb{E}(\\mathsf{x})$\" : x_rv.mean(), \n",
+    "rvdict = {\"$x_*$\" : x, \"$\\mathbb{E}(\\mathsf{x})$\" : x_rv.mean, \n",
     "          \"$\\mathsf{x}_1$\" : x_samples[0], \"$\\mathsf{x}_2$\" : x_samples[1], \"$\\mathsf{x}_3$\" : x_samples[2]}\n",
     "\n",
     "fig, axes = plt.subplots(nrows=1, ncols=2 + n_samples, figsize=(8, 2.5), sharey=True)\n",

docs/source/tutorials/adaptive_steps_odefilter.ipynb

@@ -982,7 +982,7 @@
     "ivp = lotkavolterra([0.0, 20.0], initrv)\n",
     "sol = probsolve_ivp(ivp, tol=1e-1)\n",
     "\n",
-    "means = sol.y.mean()\n",
+    "means = sol.y.mean\n",
     "\n",
     "plt.plot(sol.t, means[:, 0], label=\"Prey\")\n",
     "plt.plot(sol.t, means[:, 1], label=\"Predators\")\n",
@@ -19320,7 +19320,7 @@
     "            which_prior=priors[idx],\n",
     "            method=filters[jdx],\n",
     "        )\n",
-    "        ts, means, stds = sol.t, sol.y.mean(), sol.y.std()\n",
+    "        ts, means, stds = sol.t, sol.y.mean, sol.y.std()\n",
     "        ax[idx][jdx].plot(ts, means)\n",
     "        ax[idx][jdx].fill_between(\n",
     "            ts, means[:, 0] - 3 * stds[:, 0], means[:, 0] + 3 * stds[:, 0], alpha=0.25\n",
@@ -65242,7 +65242,7 @@
     "            which_prior=priors[idx],\n",
     "            method=filters[jdx],\n",
     "        )\n",
-    "        ts, means, stds = sol.t, sol.y.mean(), sol.y.std()\n",
+    "        ts, means, stds = sol.t, sol.y.mean, sol.y.std()\n",
     "        ax[idx][jdx].plot(ts, means)\n",
     "        ax[idx][jdx].fill_between(\n",
     "            ts, means[:, 0] - 3 * stds[:, 0], means[:, 0] + 3 * stds[:, 0], alpha=0.25\n",
@@ -83739,7 +83739,7 @@
     "            which_prior=priors[idx],\n",
     "            method=filters[jdx],\n",
     "        )\n",
-    "        ts, means, stds = sol.t, sol.y.mean(), sol.y.std()\n",
+    "        ts, means, stds = sol.t, sol.y.mean, sol.y.std()\n",
     "        ax[idx][jdx].plot(ts, means)\n",
     "        ax[idx][jdx].fill_between(\n",
     "            ts, means[:, 0] - 3 * stds[:, 0], means[:, 0] + 3 * stds[:, 0], alpha=0.25\n",

docs/source/tutorials/galerkin_method.ipynb

@@ -690,7 +690,7 @@
     "print(info)\n",
     "\n",
     "# Save solution to file\n",
-    "np.save(data_path + \"solution_res{}.npy\".format(mesh_resolution_coarse), uhat.mean())"
+    "np.save(data_path + \"solution_res{}.npy\".format(mesh_resolution_coarse), uhat.mean)"
    ]
   },
   {
@@ -3350,8 +3350,8 @@
     "\n",
     "# Interpolation levels and color scaling\n",
     "interpol_levels = None\n",
-    "vmin_tricont = np.min(uhat.mean().ravel())\n",
-    "vmax_tricont = np.max(uhat.mean().ravel())\n",
+    "vmin_tricont = np.min(uhat.mean.ravel())\n",
+    "vmax_tricont = np.max(uhat.mean.ravel())\n",
     "\n",
     "# Sample from Posterior\n",
     "np.random.seed(42)\n",
@@ -3368,7 +3368,7 @@
     "\n",
     "# Estimated solution\n",
     "fig.add_subplot(gs00[0, 0])\n",
-    "plt.tricontourf(triang, uhat.mean().ravel(), levels=interpol_levels, vmin=vmin_tricont, vmax=vmax_tricont)\n",
+    "plt.tricontourf(triang, uhat.mean.ravel(), levels=interpol_levels, vmin=vmin_tricont, vmax=vmax_tricont)\n",
     "plt.triplot(triang, linewidth=2.5, color=\"white\", alpha=.25)\n",
     "plt.title(\"$\\mathbb{E}[\\\\bm{\\\\mathsf{u}}]$\", fontsize=24)\n",
     "plt.axis('scaled')\n",
@@ -3732,8 +3732,8 @@
    ],
    "source": [
     "# Plot linear operators\n",
-    "matdict = {\"$A$\" : A.todense(), \"$\\mathbb{E}(\\mathsf{A})$\" : Ahat.mean().todense(), \n",
-    "           \"$\\mathbb{E}(\\mathsf{H})$\" : Ainvhat.mean().todense()}\n",
+    "matdict = {\"$A$\" : A.todense(), \"$\\mathbb{E}(\\mathsf{A})$\" : Ahat.mean.todense(), \n",
+    "           \"$\\mathbb{E}(\\mathsf{H})$\" : Ainvhat.mean.todense()}\n",
     "\n",
     "# Set imshow limits uniform across matrices\n",
     "vmax = np.max([np.max(np.abs(mat)) for mat in list(matdict.values())])\n",
@@ -3792,7 +3792,7 @@
     "rel_error_u = error_u / np.linalg.norm(u)\n",
     "\n",
     "# Scaled error\n",
-    "Sigma = uhat.cov().todense()\n",
+    "Sigma = uhat.cov.todense()\n",
     "Sigma_pseudoinv = np.linalg.pinv(Sigma, hermitian=True)\n",
     "scaled_error_u = np.real_if_close(scipy.linalg.sqrtm(Sigma_pseudoinv) @ error_u, tol=10**8)"
    ]

docs/source/tutorials/linear_operators.ipynb

@@ -1045,7 +1045,7 @@
    ],
    "source": [
     "# Plot\n",
-    "rvdict = {\"$\\mathbb{E}(\\mathsf{Y})$\" : Y.mean().todense(), \n",
+    "rvdict = {\"$\\mathbb{E}(\\mathsf{Y})$\" : Y.mean.todense(), \n",
     "          \"$\\mathsf{Y}_1$\" : Ysamples[2], \"$\\mathsf{Y}_2$\" : Ysamples[1], \"$\\mathsf{Y}_3$\" : Ysamples[2]}\n",
     "vmin = np.min([np.min(mat) for mat in list(rvdict.values())])\n",
     "vmax = np.max([np.max(mat) for mat in list(rvdict.values())])\n",

docs/source/tutorials/linear_systems.ipynb

@@ -758,7 +758,7 @@
     "Ainv_samples = Ainv.sample(3)\n",
     "\n",
     "# Plot\n",
-    "rvdict = {\"$A$\" : A, \"$\\mathbb{E}(\\mathsf{A})$\" : Ahat.mean().todense(),\n",
+    "rvdict = {\"$A$\" : A, \"$\\mathbb{E}(\\mathsf{A})$\" : Ahat.mean.todense(),\n",
     "          \"$\\mathsf{A}_1$\" : Ahat_samples[0], \"$\\mathsf{A}_2$\" : Ahat_samples[1]}\n",
     "\n",
     "fig, axes = plt.subplots(nrows=1, ncols=len(rvdict), figsize=(10, 3), sharey=True)\n",
@@ -1140,7 +1140,7 @@
    ],
    "source": [
     "# Plot\n",
-    "rvdict = {\"$A^{-1}\": np.linalg.inv(A), \"$\\mathbb{E}(\\mathsf{H})\" : Ainv.mean().todense(), \n",
+    "rvdict = {\"$A^{-1}\": np.linalg.inv(A), \"$\\mathbb{E}(\\mathsf{H})\" : Ainv.mean.todense(), \n",
     "          \"$\\mathsf{H}_1\" : Ainv_samples[0], \"$\\mathsf{H}_2\" : Ainv_samples[1]}\n",
     "\n",
     "fig, axes = plt.subplots(nrows=2, ncols=len(rvdict), figsize=(10, 6), sharey=True)\n",

docs/source/tutorials/uncertainties_odefilter.ipynb

@@ -201417,7 +201417,7 @@
     }
    ],
    "source": [
-    "means, stds = sol.y.mean(), sol.y.std()\n",
+    "means, stds = sol.y.mean, sol.y.std()\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True)\n",
     "ax1.plot(sol.t, means[:, 0], label=\"x1\")\n",
@@ -402823,7 +402823,7 @@
     }
    ],
    "source": [
-    "means, stds = sol.y.mean(), sol.y.std()\n",
+    "means, stds = sol.y.mean, sol.y.std()\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True)\n",
     "ax1.plot(sol.t, means[:, 0], label=\"x1\")\n",

src/probnum/diffeq/ode/ivp.py

@@ -61,7 +61,7 @@ def hess(t, y):
         return log_hess(t, y, params)
 
     def sol(t):
-        return log_sol(t, params, initrv.mean())
+        return log_sol(t, params, initrv.mean)
 
     return IVP(timespan, initrv, rhs, jac, hess, sol)
 
@@ -264,7 +264,7 @@ class IVP(ODE):
     >>> from probnum.diffeq import IVP
     >>> rhsfun = lambda t, y, **kwargs: 2.0*y
     >>> from probnum.prob import Dirac, Normal, RandomVariable
-    >>> initrv = RandomVariable(distribution=Dirac(0.1))
+    >>> initrv = Dirac(0.1)
     >>> timespan = (0, 10)
     >>> ivp = IVP(timespan, initrv, rhsfun)
     >>> print(ivp.rhs(0., 2.))
@@ -274,7 +274,7 @@ class IVP(ODE):
     >>> print(ivp.t0)
     0
 
-    >>> initrv = RandomVariable(distribution=Normal(0.1, 1.0))
+    >>> initrv = Normal(0.1, 1.0)
     >>> ivp = IVP(timespan, initrv, rhsfun)
     >>> jac = lambda t, y, **kwargs: 2.0
     >>> ivp = IVP(timespan, initrv, rhs=rhsfun, jac=jac)
@@ -310,7 +310,7 @@ def ndim(self):
 
         Depends on the mean of the initial random variable.
         """
-        if np.isscalar(self.initrv.mean()):
+        if np.isscalar(self.initrv.mean):
             return 1
         else:
-            return len(self.initrv.mean())
+            return len(self.initrv.mean)

src/probnum/diffeq/odefiltsmooth/ivp2filter.py

@@ -7,8 +7,7 @@
 
 from probnum.filtsmooth import ExtendedKalman, UnscentedKalman
 from probnum.filtsmooth.statespace.discrete import DiscreteGaussianModel
-from probnum.prob import RandomVariable
-from probnum.prob.distributions import Normal, Dirac
+from probnum.prob import Normal, Dirac
 
 
 def ivp2ekf0(ivp, prior, evlvar):
@@ -167,10 +166,10 @@ def _initialdistribution(ivp, prior):
     Note that the projection matrices :math:`H_0` and :math:`H_1`
     become :math:`H_0 P^{-1}` and :math:`H_1 P^{-1}`.
     """
-    if not issubclass(type(ivp.initrv.distribution), Normal):
-        if not issubclass(type(ivp.initrv.distribution), Dirac):
+    if not issubclass(type(ivp.initrv), Normal):
+        if not issubclass(type(ivp.initrv), Dirac):
             raise RuntimeError("Initial distribution not Normal nor Dirac")
-    x0 = ivp.initialdistribution.mean()
+    x0 = ivp.initialdistribution.mean
     dx0 = ivp.rhs(ivp.t0, x0)
     ddx0 = _ddx(ivp.t0, x0, ivp)
     dddx0 = _dddx(ivp.t0, x0, ivp)
@@ -194,17 +193,17 @@ def _initialdistribution(ivp, prior):
         projmat = np.hstack((h0.T, h1.T, h2.T, h3.T)).T
         data = np.hstack((x0, dx0, ddx0, dddx0))
         _size = 4
-    largecov = np.kron(np.eye(_size), ivp.initialdistribution.cov())
+    largecov = np.kron(np.eye(_size), ivp.initialdistribution.cov)
     s = projmat @ initcov @ projmat.T + largecov
     crosscov = initcov @ projmat.T
     newmean = crosscov @ np.linalg.solve(s, data)
     newcov = initcov - (crosscov @ np.linalg.solve(s.T, crosscov.T)).T
-    return RandomVariable(distribution=Normal(newmean, newcov))
+    return Normal(newmean, newcov)
 
 
 def _initialdistribution_no_precond(ivp, prior):
-
-    x0 = ivp.initialdistribution.mean()
+    """ """
+    x0 = ivp.initialdistribution.mean
     dx0 = ivp.rhs(ivp.t0, x0)
     ddx0 = _ddx(ivp.t0, x0, ivp)
     h0 = prior.proj2coord(coord=0)
@@ -221,7 +220,7 @@ def _initialdistribution_no_precond(ivp, prior):
     crosscov = initcov @ projmat.T  # @ np.linalg.inv(s)
     newmean = crosscov @ np.linalg.solve(s, data)
     newcov = initcov - (crosscov @ np.linalg.solve(s, crosscov)).T
-    return RandomVariable(distribution=Normal(newmean, newcov))
+    return Normal(newmean, newcov)
 
 
 def _ddx(t, x, ivp):

src/probnum/diffeq/odefiltsmooth/ivpfiltsmooth.py

@@ -1,8 +1,7 @@
 import warnings
 import numpy as np
 
-from probnum.prob import RandomVariable
-from probnum.prob.distributions import Normal
+from probnum.prob import Normal
 from probnum.diffeq import odesolver
 from probnum.diffeq.odefiltsmooth.prior import ODEPrior
 from probnum.diffeq.odesolution import ODESolution
@@ -67,7 +66,7 @@ def solve(self, firststep, **kwargs):
             )
 
             errorest, ssq = self._estimate_error(
-                filt_rv.mean(), crosscov, meas_cov, meas_mean
+                filt_rv.mean, crosscov, meas_cov, meas_mean
             )
 
             if self.steprule.is_accepted(step, errorest):
@@ -81,10 +80,7 @@ def solve(self, firststep, **kwargs):
             step = self._suggest_step(step, errorest)
             step = min(step, self.ivp.tmax - t)
 
-        rvs = [
-            RandomVariable(distribution=Normal(rv.mean(), ssqest * rv.cov()))
-            for rv in rvs
-        ]
+        rvs = [Normal(rv.mean, ssqest * rv.cov) for rv in rvs]
 
         return ODESolution(times, rvs, self)
 
@@ -116,9 +112,9 @@ def odesmooth(self, filter_solution, **kwargs):
 
     def undo_preconditioning(self, rv):
         ipre = self.gfilt.dynamicmodel.invprecond
-        newmean = ipre @ rv.mean()
-        newcov = ipre @ rv.cov() @ ipre.T
-        newrv = RandomVariable(distribution=Normal(newmean, newcov))
+        newmean = ipre @ rv.mean
+        newcov = ipre @ rv.cov @ ipre.T
+        newrv = Normal(newmean, newcov)
         return newrv
 
     def _estimate_error(self, currmn, ccest, covest, mnest):

src/probnum/diffeq/odefiltsmooth/odefiltsmooth.py

@@ -153,8 +153,8 @@ def probsolve_ivp(
         y : :obj:`list` of :obj:`RandomVariable`, length=N
             Discrete-time solution at times :math:`t_1, ..., t_N`,
             as a list of random variables.
-            The means and covariances can be accessed with ``solution.y.mean()``
-            and ``solution.y.cov()``.
+            The means and covariances can be accessed with ``solution.y.mean``
+            and ``solution.y.cov``.
 
     See Also
     --------
@@ -182,8 +182,8 @@ def probsolve_ivp(
     Examples
     --------
     >>> from probnum.diffeq import logistic, probsolve_ivp
-    >>> from probnum.prob import RandomVariable, Dirac, Normal
-    >>> initrv = RandomVariable(distribution=Dirac(0.15))
+    >>> from probnum.prob import Dirac, Normal
+    >>> initrv = Dirac(0.15)
     >>> ivp = logistic(timespan=[0., 1.5], initrv=initrv, params=(4, 1))
     >>> solution = probsolve_ivp(ivp, method="ekf0", step=0.1)
     >>> print(solution.y.mean())
@@ -204,7 +204,7 @@ def probsolve_ivp(
      [0.97577054]
      [0.9831919 ]]
 
-    >>> initrv = RandomVariable(distribution=Dirac(0.15))
+    >>> initrv = Dirac(0.15)
     >>> ivp = logistic(timespan=[0., 1.5], initrv=initrv, params=(4, 1))
     >>> solution = probsolve_ivp(ivp, method="eks1", which_prior="ioup3", step=0.1)
     >>> print(solution.y.mean())