From acc4fe99aecb05aee9f40d08df6d54306a1435ad Mon Sep 17 00:00:00 2001 From: mjhajharia Date: Tue, 25 Jan 2022 20:30:38 +0530 Subject: [PATCH 1/4] add RV for AR --- pymc/distributions/timeseries.py | 76 ++++++++++++++------------------ 1 file changed, 34 insertions(+), 42 deletions(-) diff --git a/pymc/distributions/timeseries.py b/pymc/distributions/timeseries.py index dbdd591ba5..4e47dfe99d 100644 --- a/pymc/distributions/timeseries.py +++ b/pymc/distributions/timeseries.py @@ -22,6 +22,8 @@ from pymc.distributions.continuous import Flat, Normal, get_tau_sigma from pymc.distributions.shape_utils import to_tuple +from aesara.tensor.random.op import RandomVariable + __all__ = [ "AR1", "AR", @@ -32,49 +34,39 @@ "MvStudentTRandomWalk", ] +class ARRV(RandomVariable): + name = "AR" + ndim_supp = 0 + ndims_params = [1, 0, 0] + dtype = "floatX" + _print_name = ("AR", "\\operatorname{AR}") + + #set default values for parameters + def __call__(self, phi, mu=0.0, sigma=1.0, init=None **kwargs) -> TensorVariable: + return super().__call__(phi, mu, sigma,init **kwargs) + + @classmethod + def rng_fn( + cls, + rng: np.random.default_rng(), + phi: np.ndarray, + mu: np.ndarray, + sigma: np.ndarray, + size: Tuple[int, ...], + init: float) -> np.ndarray: + + # size parameter *should* be necessary for time series generation + if size is None: + raise ValueError('Specify size') + + # sign conventions used in signal.lfilter (or signal processing) + phi = np.r_[1, -phi] # add zero-lag and negate + + #ifilter convolutes x with the coefficient theta to create a linear recurrence relation and generate the AR process + if init is None: + init = rng.normal(loc=mu, scale=sigma,size=size) + return signal.lfilter(phi, init, axis=-1) -class AR1(distribution.Continuous): - """ - Autoregressive process with 1 lag. - - Parameters - ---------- - k: tensor - effect of lagged value on current value - tau_e: tensor - precision for innovations - """ - - def __init__(self, k, tau_e, *args, **kwargs): - super().__init__(*args, **kwargs) - self.k = k = at.as_tensor_variable(k) - self.tau_e = tau_e = at.as_tensor_variable(tau_e) - self.tau = tau_e * (1 - k ** 2) - self.mode = at.as_tensor_variable(0.0) - - def logp(self, x): - """ - Calculate log-probability of AR1 distribution at specified value. - - Parameters - ---------- - x: numeric - Value for which log-probability is calculated. - - Returns - ------- - TensorVariable - """ - k = self.k - tau_e = self.tau_e # innovation precision - tau = tau_e * (1 - k ** 2) # ar1 precision - - x_im1 = x[:-1] - x_i = x[1:] - boundary = Normal.dist(0.0, tau=tau).logp - - innov_like = Normal.dist(k * x_im1, tau=tau_e).logp(x_i) - return boundary(x[0]) + at.sum(innov_like) class AR(distribution.Continuous): From ced3de958c9df080884acc07b3a5660dff7c5293 Mon Sep 17 00:00:00 2001 From: mjhajharia Date: Wed, 26 Jan 2022 00:32:26 +0530 Subject: [PATCH 2/4] Revert "add RV for AR" This reverts commit acc4fe99aecb05aee9f40d08df6d54306a1435ad. --- pymc/distributions/timeseries.py | 76 ++++++++++++++++++-------------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/pymc/distributions/timeseries.py b/pymc/distributions/timeseries.py index 4e47dfe99d..dbdd591ba5 100644 --- a/pymc/distributions/timeseries.py +++ b/pymc/distributions/timeseries.py @@ -22,8 +22,6 @@ from pymc.distributions.continuous import Flat, Normal, get_tau_sigma from pymc.distributions.shape_utils import to_tuple -from aesara.tensor.random.op import RandomVariable - __all__ = [ "AR1", "AR", @@ -34,39 +32,49 @@ "MvStudentTRandomWalk", ] -class ARRV(RandomVariable): - name = "AR" - ndim_supp = 0 - ndims_params = [1, 0, 0] - dtype = "floatX" - _print_name = ("AR", "\\operatorname{AR}") - - #set default values for parameters - def __call__(self, phi, mu=0.0, sigma=1.0, init=None **kwargs) -> TensorVariable: - return super().__call__(phi, mu, sigma,init **kwargs) - - @classmethod - def rng_fn( - cls, - rng: np.random.default_rng(), - phi: np.ndarray, - mu: np.ndarray, - sigma: np.ndarray, - size: Tuple[int, ...], - init: float) -> np.ndarray: - - # size parameter *should* be necessary for time series generation - if size is None: - raise ValueError('Specify size') - - # sign conventions used in signal.lfilter (or signal processing) - phi = np.r_[1, -phi] # add zero-lag and negate - - #ifilter convolutes x with the coefficient theta to create a linear recurrence relation and generate the AR process - if init is None: - init = rng.normal(loc=mu, scale=sigma,size=size) - return signal.lfilter(phi, init, axis=-1) +class AR1(distribution.Continuous): + """ + Autoregressive process with 1 lag. + + Parameters + ---------- + k: tensor + effect of lagged value on current value + tau_e: tensor + precision for innovations + """ + + def __init__(self, k, tau_e, *args, **kwargs): + super().__init__(*args, **kwargs) + self.k = k = at.as_tensor_variable(k) + self.tau_e = tau_e = at.as_tensor_variable(tau_e) + self.tau = tau_e * (1 - k ** 2) + self.mode = at.as_tensor_variable(0.0) + + def logp(self, x): + """ + Calculate log-probability of AR1 distribution at specified value. + + Parameters + ---------- + x: numeric + Value for which log-probability is calculated. + + Returns + ------- + TensorVariable + """ + k = self.k + tau_e = self.tau_e # innovation precision + tau = tau_e * (1 - k ** 2) # ar1 precision + + x_im1 = x[:-1] + x_i = x[1:] + boundary = Normal.dist(0.0, tau=tau).logp + + innov_like = Normal.dist(k * x_im1, tau=tau_e).logp(x_i) + return boundary(x[0]) + at.sum(innov_like) class AR(distribution.Continuous): From d16370c0d7de9550a79bae6ee4ba0202953220a0 Mon Sep 17 00:00:00 2001 From: mjhajharia Date: Wed, 26 Jan 2022 00:33:49 +0530 Subject: [PATCH 3/4] AR RV --- pymc/distributions/timeseries.py | 79 +++++++++++++++----------------- 1 file changed, 37 insertions(+), 42 deletions(-) diff --git a/pymc/distributions/timeseries.py b/pymc/distributions/timeseries.py index dbdd591ba5..cba965ad1a 100644 --- a/pymc/distributions/timeseries.py +++ b/pymc/distributions/timeseries.py @@ -22,6 +22,11 @@ from pymc.distributions.continuous import Flat, Normal, get_tau_sigma from pymc.distributions.shape_utils import to_tuple +from aesara.tensor.random.op import RandomVariable +from pymc.aesaraf import floatX, intX +from aesara.tensor.var import TensorVariable +from typing import Tuple + __all__ = [ "AR1", "AR", @@ -33,49 +38,39 @@ ] -class AR1(distribution.Continuous): - """ - Autoregressive process with 1 lag. - - Parameters - ---------- - k: tensor - effect of lagged value on current value - tau_e: tensor - precision for innovations - """ - - def __init__(self, k, tau_e, *args, **kwargs): - super().__init__(*args, **kwargs) - self.k = k = at.as_tensor_variable(k) - self.tau_e = tau_e = at.as_tensor_variable(tau_e) - self.tau = tau_e * (1 - k ** 2) - self.mode = at.as_tensor_variable(0.0) - - def logp(self, x): - """ - Calculate log-probability of AR1 distribution at specified value. - - Parameters - ---------- - x: numeric - Value for which log-probability is calculated. - - Returns - ------- - TensorVariable - """ - k = self.k - tau_e = self.tau_e # innovation precision - tau = tau_e * (1 - k ** 2) # ar1 precision - - x_im1 = x[:-1] - x_i = x[1:] - boundary = Normal.dist(0.0, tau=tau).logp - - innov_like = Normal.dist(k * x_im1, tau=tau_e).logp(x_i) - return boundary(x[0]) + at.sum(innov_like) +class ARrv(RandomVariable): + name = "AR" + ndim_supp = 0 + ndims_params = [0, 0, 0, 0, 0] + dtype = "floatX" + _print_name = ("AR", "\\operatorname{AR}") + + + def __call__(self, phi, init=0.0, mu=0.0, sigma=1.0, **kwargs) -> TensorVariable: + return super().__call__(phi, init, mu, sigma, **kwargs) + + + @classmethod + def rng_fn( + cls, + rng: np.random.default_rng(), + phi: float, + init: float, + mu: np.ndarray, + sigma: np.ndarray, + size: Tuple[int, ...], + ) -> np.ndarray: + + phi = np.asarray(phi) + if init!=0.0: + y = np.asarray([init]) + else: + y = rng.normal(0.0, 1.0,size=1) + for i in range(1, size[0]): + y = np.append(y, mu + phi*y[-1] + rng.normal(0.0, sigma)) + + return y class AR(distribution.Continuous): r""" From 75f865ac50ce1e45f8c284ab4fc2f84409d6ddc7 Mon Sep 17 00:00:00 2001 From: mjhajharia Date: Wed, 26 Jan 2022 01:23:09 +0530 Subject: [PATCH 4/4] trying stuff --- pymc/distributions/AR1.ipynb | 249 +++++++++++++++++++++++++++++++ pymc/distributions/timeseries.py | 69 ++++++++- 2 files changed, 311 insertions(+), 7 deletions(-) create mode 100644 pymc/distributions/AR1.ipynb diff --git a/pymc/distributions/AR1.ipynb b/pymc/distributions/AR1.ipynb new file mode 100644 index 0000000000..3cd44d8e97 --- /dev/null +++ b/pymc/distributions/AR1.ipynb @@ -0,0 +1,249 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "0df772d5", + "metadata": {}, + "outputs": [], + "source": [ + "from aesara.tensor.random.op import RandomVariable\n", + "from pymc.aesaraf import floatX, intX\n", + "from aesara.tensor.var import TensorVariable\n", + "from typing import Tuple\n", + "\n", + "import numpy as np\n", + "import aesara.tensor as at\n", + "import pymc as pm\n", + "from pymc.distributions import distribution, multivariate\n", + "from pymc.distributions.continuous import Flat, Normal, HalfNormal, Uniform, get_tau_sigma" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e863a50e", + "metadata": {}, + "outputs": [], + "source": [ + "class AR1rv(RandomVariable):\n", + " name = \"AR1\"\n", + " ndim_supp = 0\n", + " ndims_params = [0, 0, 0, 0, 0]\n", + " dtype = \"floatX\"\n", + " _print_name = (\"AR1\", \"\\\\operatorname{AR1}\")\n", + " \n", + " \n", + " def __call__(self, phi, init=0.0, mu=0.0, sigma=1.0,*args,**kwargs) -> TensorVariable:\n", + " return super().__call__(phi, init, mu, sigma,*args, **kwargs)\n", + "\n", + " \n", + " @classmethod\n", + " def rng_fn(\n", + " cls,\n", + " rng: np.random.default_rng(),\n", + " phi: float,\n", + " init: float,\n", + " mu: np.ndarray,\n", + " sigma: np.ndarray,\n", + " size: Tuple[int, ...],\n", + " ) -> np.ndarray:\n", + "\n", + " phi = np.asarray(phi)\n", + " \n", + " if init!=0.0:\n", + " y = np.asarray([init])\n", + " else:\n", + " y = rng.normal(0.0, 1.0,size=1)\n", + " \n", + " y[-1]\n", + " for i in range(1, size[0]):\n", + " y = np.append(y, mu + phi*y[-1] + rng.normal(0.0, sigma))\n", + " \n", + " return y" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ec0e61e7", + "metadata": {}, + "outputs": [], + "source": [ + "class AR1(distribution.Continuous):\n", + " rv_op = AR1rv()\n", + " r\"\"\"\n", + " Autoregressive with order 1.\n", + "\n", + " Parameters\n", + " ----------\n", + " phi: tensor\n", + " Autoregressive coefficient.\n", + " sigma: float\n", + " Standard deviation of noise (sigma > 0, default: 1.0).\n", + " mu: float\n", + " Mean or constant (default: 0.0)\n", + " size: int\n", + " Number of observations in time series.\n", + " \n", + " Examples\n", + " --------\n", + " \n", + " .. code-block:: python\n", + " \n", + " #Generate an AR1 process\n", + " Y = AR1.dist(phi=.50, size=500)\n", + " \n", + " #Estimate AR1 parameter\n", + " with pm.Model() as ar1:\n", + " \n", + " #priors\n", + " phi = Uniform(\"phi\",-0.99,0.99,shape=1)\n", + " \n", + " #sampling and results\n", + " ar1 = AR1(\"AR1\", phi=phi, observed = Y)\n", + " trace = pm.sample(1000, tune=500, target_accept=0.99,return_inferencedata=True)\n", + " az.summary(trace)\n", + " \n", + " \"\"\"\n", + " @classmethod\n", + " def dist(cls, phi,*args,**kwargs):\n", + " return super().dist([phi],*args,**kwargs)\n", + "\n", + " def logp(obs, phi, mu=0.0, sigma=1.0, *args,**kwargs):\n", + " \n", + " n = obs.shape[0]\n", + " epsilon = at.ivector()\n", + " \n", + " for i in range(1, n):\n", + " epsilon = at.sub(Y, mu + phi*obs[-1])\n", + " \n", + " logp = at.sum(pm.logp(Normal.dist(mu=nu, sigma=sigma,size=n), epsilon))\n", + "\n", + " return logp\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "b193b123", + "metadata": {}, + "outputs": [], + "source": [ + "Y = AR1.dist(phi=[.50], size=500, init=0.3, mu=0.0, sigma=1.0).eval()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "7785438b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " SARIMAX Results \n", + "==============================================================================\n", + "Dep. Variable: y No. Observations: 500\n", + "Model: SARIMAX(1, 0, 0) Log Likelihood -694.258\n", + "Date: Wed, 26 Jan 2022 AIC 1392.516\n", + "Time: 01:18:13 BIC 1400.945\n", + "Sample: 0 HQIC 1395.823\n", + " - 500 \n", + "Covariance Type: opg \n", + "==============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "------------------------------------------------------------------------------\n", + "ar.L1 0.5082 0.041 12.388 0.000 0.428 0.589\n", + "sigma2 0.9404 0.061 15.390 0.000 0.821 1.060\n", + "===================================================================================\n", + "Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 0.22\n", + "Prob(Q): 0.89 Prob(JB): 0.89\n", + "Heteroskedasticity (H): 1.04 Skew: -0.01\n", + "Prob(H) (two-sided): 0.81 Kurtosis: 2.90\n", + "===================================================================================\n", + "\n", + "Warnings:\n", + "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n" + ] + } + ], + "source": [ + "import statsmodels.api as sm\n", + "\n", + "mod = sm.tsa.statespace.SARIMAX(Y, order=(1, 0, 0))\n", + "\n", + "res_mle = mod.fit(disp=False)\n", + "print(res_mle.summary())" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "dc317dae", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n" + ] + }, + { + "ename": "TypeError", + "evalue": "Cannot convert Type TensorType(float64, vector) (of Variable TensorConstant{[ 3.000000..07616e-01]}) into Type TensorType(float64, (True,)). You can try to manually convert TensorConstant{[ 3.000000..07616e-01]} into a TensorType(float64, (True,)).", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/var/folders/sd/pc07b3wn65nflgx8wpkkr7wr0000gn/T/ipykernel_4857/805995269.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mphi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mUniform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"phi\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m0.99\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0.99\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mAR1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAR1\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"AR1\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mphi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mphi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobserved\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtune\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtarget_accept\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.99\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mreturn_inferencedata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/sampling.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(draws, step, init, n_init, initvals, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, **kwargs)\u001b[0m\n\u001b[1;32m 492\u001b[0m \u001b[0;31m# By default, try to use NUTS\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 493\u001b[0m \u001b[0m_log\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Auto-assigning NUTS sampler...\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 494\u001b[0;31m initial_points, step = init_nuts(\n\u001b[0m\u001b[1;32m 495\u001b[0m \u001b[0minit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minit\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 496\u001b[0m \u001b[0mchains\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mchains\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/sampling.py\u001b[0m in \u001b[0;36minit_nuts\u001b[0;34m(init, chains, n_init, model, seeds, progressbar, jitter_max_retries, tune, initvals, **kwargs)\u001b[0m\n\u001b[1;32m 2252\u001b[0m ]\n\u001b[1;32m 2253\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2254\u001b[0;31m initial_points = _init_jitter(\n\u001b[0m\u001b[1;32m 2255\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2256\u001b[0m \u001b[0minitvals\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/sampling.py\u001b[0m in \u001b[0;36m_init_jitter\u001b[0;34m(model, initvals, seeds, jitter, jitter_max_retries)\u001b[0m\n\u001b[1;32m 2138\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mjitter_max_retries\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2139\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2140\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_start_vals\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2141\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mSamplingError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2142\u001b[0m \u001b[0;31m# Retry with a new seed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/model.py\u001b[0m in \u001b[0;36mcheck_start_vals\u001b[0;34m(self, start)\u001b[0m\n\u001b[1;32m 1667\u001b[0m )\n\u001b[1;32m 1668\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1669\u001b[0;31m \u001b[0minitial_eval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpoint_logps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0melem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1670\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1671\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitial_eval\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/model.py\u001b[0m in \u001b[0;36mpoint_logps\u001b[0;34m(self, point, round_vals)\u001b[0m\n\u001b[1;32m 1707\u001b[0m for factor, factor_logp in zip(\n\u001b[1;32m 1708\u001b[0m \u001b[0mfactors\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1709\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfactor\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mfactor\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogp_elemwiset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfactors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1710\u001b[0m )\n\u001b[1;32m 1711\u001b[0m },\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/model.py\u001b[0m in \u001b[0;36mlogp_elemwiset\u001b[0;34m(self, vars, jacobian)\u001b[0m\n\u001b[1;32m 774\u001b[0m \u001b[0mrv_logps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 775\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mrv_values\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 776\u001b[0;31m \u001b[0mrv_logps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlogpt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrv_values\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrv_values\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjacobian\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mjacobian\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 777\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrv_logps\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 778\u001b[0m \u001b[0mrv_logps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mrv_logps\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/pymc/distributions/logprob.py\u001b[0m in \u001b[0;36mlogpt\u001b[0;34m(var, rv_values, jacobian, scaling, transformed, sum, **kwargs)\u001b[0m\n\u001b[1;32m 223\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 224\u001b[0m \u001b[0mtransform_opt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTransformValuesOpt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtransform_map\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 225\u001b[0;31m temp_logp_var_dict = factorized_joint_logprob(\n\u001b[0m\u001b[1;32m 226\u001b[0m \u001b[0mtmp_rvs_to_values\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextra_rewrites\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtransform_opt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0muse_jacobian\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mjacobian\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 227\u001b[0m )\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aeppl/joint_logprob.py\u001b[0m in \u001b[0;36mfactorized_joint_logprob\u001b[0;34m(rv_values, warn_missing_rvs, extra_rewrites, **kwargs)\u001b[0m\n\u001b[1;32m 177\u001b[0m \u001b[0;31m# Also, store the results in the `replacements` map for the nodes\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 178\u001b[0m \u001b[0;31m# that follow.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 179\u001b[0;31m remapped_vars, _ = rvs_to_value_vars(\n\u001b[0m\u001b[1;32m 180\u001b[0m \u001b[0mq_rv_value_vars\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 181\u001b[0m \u001b[0minitial_replacements\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreplacements\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aeppl/utils.py\u001b[0m in \u001b[0;36mrvs_to_value_vars\u001b[0;34m(graphs, initial_replacements, **kwargs)\u001b[0m\n\u001b[1;32m 222\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 223\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 224\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mreplace_rvs_in_graphs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgraphs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreplace_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_replacements\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 225\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 226\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aeppl/utils.py\u001b[0m in \u001b[0;36mreplace_rvs_in_graphs\u001b[0;34m(graphs, replacement_fn, initial_replacements, **kwargs)\u001b[0m\n\u001b[1;32m 188\u001b[0m )\n\u001b[1;32m 189\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 190\u001b[0;31m \u001b[0mfg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreplace_all\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreplacements\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mimport_missing\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 191\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 192\u001b[0m \u001b[0mgraphs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aesara/graph/fg.py\u001b[0m in \u001b[0;36mreplace_all\u001b[0;34m(self, pairs, **kwargs)\u001b[0m\n\u001b[1;32m 552\u001b[0m \u001b[0;34m\"\"\"Replace variables in the `FunctionGraph` according to ``(var, new_var)`` pairs in a list.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 553\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mvar\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_var\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpairs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 554\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreplace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_var\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 555\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 556\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mattach_feature\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeature\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mFeature\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aesara/graph/fg.py\u001b[0m in \u001b[0;36mreplace\u001b[0;34m(self, var, new_var, reason, verbose, import_missing)\u001b[0m\n\u001b[1;32m 513\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"optimizer: rewrite {reason} replaces {var} with {new_var}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 514\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 515\u001b[0;31m \u001b[0mnew_var\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvar\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilter_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_var\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallow_convert\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 516\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 517\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mvar\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/anaconda3/envs/pymc4/lib/python3.9/site-packages/aesara/tensor/type.py\u001b[0m in \u001b[0;36mfilter_variable\u001b[0;34m(self, other, allow_convert)\u001b[0m\n\u001b[1;32m 245\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mother2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 247\u001b[0;31m raise TypeError(\n\u001b[0m\u001b[1;32m 248\u001b[0m \u001b[0;34mf\"Cannot convert Type {other.type} \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 249\u001b[0m \u001b[0;34mf\"(of Variable {other}) into Type {self}. \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mTypeError\u001b[0m: Cannot convert Type TensorType(float64, vector) (of Variable TensorConstant{[ 3.000000..07616e-01]}) into Type TensorType(float64, (True,)). You can try to manually convert TensorConstant{[ 3.000000..07616e-01]} into a TensorType(float64, (True,))." + ] + } + ], + "source": [ + "with pm.Model() as AR:\n", + " phi = Uniform(\"phi\",-0.99,0.99,shape=1)\n", + " AR1 = AR1(\"AR1\", phi=phi, observed = Y)\n", + " trace = pm.sample(1000, tune=500, target_accept=0.99,return_inferencedata=True)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pymc/distributions/timeseries.py b/pymc/distributions/timeseries.py index cba965ad1a..48f569e8f0 100644 --- a/pymc/distributions/timeseries.py +++ b/pymc/distributions/timeseries.py @@ -37,17 +37,16 @@ "MvStudentTRandomWalk", ] - -class ARrv(RandomVariable): - name = "AR" +class AR1rv(RandomVariable): + name = "AR1" ndim_supp = 0 ndims_params = [0, 0, 0, 0, 0] dtype = "floatX" - _print_name = ("AR", "\\operatorname{AR}") + _print_name = ("AR1", "\\operatorname{AR1}") - def __call__(self, phi, init=0.0, mu=0.0, sigma=1.0, **kwargs) -> TensorVariable: - return super().__call__(phi, init, mu, sigma, **kwargs) + def __call__(self, phi, init=0.0, mu=0.0, sigma=1.0,*args,**kwargs) -> TensorVariable: + return super().__call__(phi, init, mu, sigma,*args, **kwargs) @classmethod @@ -62,16 +61,72 @@ def rng_fn( ) -> np.ndarray: phi = np.asarray(phi) + if init!=0.0: y = np.asarray([init]) else: y = rng.normal(0.0, 1.0,size=1) - + + y[-1] for i in range(1, size[0]): y = np.append(y, mu + phi*y[-1] + rng.normal(0.0, sigma)) return y +class AR1(distribution.Continuous): + rv_op = AR1rv() + r""" + Autoregressive with order 1. + + Parameters + ---------- + phi: tensor + Autoregressive coefficient. + sigma: float + Standard deviation of noise (sigma > 0, default: 1.0). + mu: float + Mean or constant (default: 0.0) + size: int + Number of observations in time series. + + Examples + -------- + + .. code-block:: python + + #Generate an AR1 process + Y = AR1.dist(phi=.50, size=500) + + #Estimate AR1 parameter + with pm.Model() as ar1: + + #priors + phi = Uniform("phi",-0.99,0.99,shape=1) + + #sampling and results + ar1 = AR1("AR1", phi=phi, observed = Y) + trace = pm.sample(1000, tune=500, target_accept=0.99,return_inferencedata=True) + az.summary(trace) + + """ + @classmethod + def dist(cls, phi,*args,**kwargs): + return super().dist([phi],*args,**kwargs) + + def logp(obs, phi, mu=0.0, sigma=1.0, *args,**kwargs): + + n = obs.shape[0] + epsilon = at.ivector() + + for i in range(1, n): + epsilon = at.sub(Y, mu + phi*obs[-1]) + + logp = at.sum(pm.logp(Normal.dist(mu=nu, sigma=sigma,size=n), epsilon)) + + return logp + + + class AR(distribution.Continuous): r""" Autoregressive process with p lags.