2
$\begingroup$

I'm following a paper described in the link which models rainfall. Due to a cyclical nature the authors use a model equation (53):

$$ log(\mu_{i}) = a_{0} + a_{1} sin (\frac{2\pi i}{365}) + a_{2} cos (\frac{2\pi i}{365}) $$

where $$\mu_{i}$$ is the mean of the distribution. My code in Python is below.

import pandas as pd
import numpy as np
import math 
import statsmodels
import statsmodels.api as sm 

from patsy import dmatrices
from statsmodels.formula.api import glm


def log_y(y):
    return np.log(y)

def sin_x(x):
    return np.sin(2*math.pi*x/365.0)

def cos_x(x):
    return np.cos(2*math.pi*x/365.0)

def model():
    N = 365*10
    x = np.arange(0, N).astype(int)
    y = 0.1653 + 0.9049*np.sin(2*math.pi*x/365.0) + 2.0326*np.cos(2*math.pi*x/365.0)

    formula = 'np.log(y) ~ sin_x(x) + cos_x(x)'
    
    data = {'x': x, 'y': y}
    model = glm(formula, data, family=sm.families.Tweedie())
    result = model.fit()
    print(result.summary())

# Main function 
if __name__== "__main__": 
    model()

I have two questions:

  1. How to represent the constant term a0
  2. When I run the example above I received the error with model.fit()

ValueError: The first guess on the deviance function returned a nan. This could be a boundary problem and should be reported.

UPDATE

As Kerby Sheldon pointed out there is an issue.

Instead, I would like to estimate the parameters for the equation

$$ \mu_{i} = \exp^(a_{0} + a_{1} sin (\frac{2\pi i}{365}) + a_{2} cos (\frac{2\pi i}{365})) $$

I tried this formula

formula = 'y ~ np.exp(1 + sin_x(x) + cos_x(x))' 

but I know it's not correct. I'm not sure how to use construct this properly.

UPDATE 2 - SOLUTION

I realised I wasn't following the right approach. The actual solution in the code is:


x = np.arange(1, N + 1).astype(int)
k = np.ones(int(N))
z1 = np.sin(2*math.pi*x/365)
z2 = np.cos(2*math.pi*x/365)
df2 = pd.DataFrame({ 'k': k, 'x' : x, 'z1' : z1, 'z2' : z2})
model = sm.GLM(y, df2, family=sm.families.Tweedie(link=None,var_power = 1.6,eql=True))
result = model.fit()
$\endgroup$

1 Answer 1

2
$\begingroup$

Your y variable has negative values and then you take the log, so nans result. The nans are dropped, but I'm not sure that you want the log there. The Tweedie family defaults to a log link, but is it y or log(y) that has a Tweedie distribution? Also, you probably want to assign something to the var_power parameter or else you are basically doing Poisson regression.

$\endgroup$
4
  • $\begingroup$ The approach I wanted to raise the right-hand side to an exponential, i.e. mu_{i} = exp( ). However, I'm not sure how to create the formula with statsmodels. Also, I'm not sure how to incorporate the constant a_{0} $\endgroup$
    – sdbol
    Commented Nov 7, 2020 at 20:46
  • $\begingroup$ The constant is implicit when you use the patsy formula for statsmodels @sdbol, so it is estimated in the regression equation as you have it. You probably don't want to take the log of the left hand side here as Kerby mentions, which is estimating $\log(\mathbb{E}[\log(y)])$ here, but you probably want to estimate $\log(\mathbb{E}[y])$. $\endgroup$
    – Andy W
    Commented Nov 7, 2020 at 21:50
  • $\begingroup$ With the model y = 0.1653 + 0.9049*np.sin(2*math.pix/365.0) + 2.0326*np.cos(2*math.pix/365.0), I'm curious how I can write the formula? I tried this formula = 'y ~ np.exp(1 + sin_x(x) + cos_x(x))' but I know it's not correct. $\endgroup$
    – sdbol
    Commented Nov 9, 2020 at 11:13
  • $\begingroup$ My last comment was meant to say y = exp( ) of the term. Hence I want to understand how to create an appropriate formula. $\endgroup$
    – sdbol
    Commented Nov 9, 2020 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.