# Python statsmodels constant and error using Tweedie distribution

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()


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.
• 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])$. Commented Nov 7, 2020 at 21:50