#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Date : 2026-01-01
# Author : Lancelot PINCET
# GitHub : https://github.com/LancelotPincet
# Library : funcLP
# %% Libraries
from funclp import Distribution, Parameter, ufunc
import numpy as np
import math
PARAMETERS = [Parameter('k', np.float32(1.)), Parameter('eps', np.float32(1e-6))]
# %% Function
[docs]
class Gamma(Distribution) :
@property
def default_attributes(self):
return 'nb.float32(1.), nb.float32(1e-6), ' # k, eps
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def pdf(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Probability Density Function."""
if raw_data < 0 :
return 0
x = raw_data
theta = eps if model_data < eps else model_data
return (x ** (k - 1) * math.exp(-x / theta) / (math.gamma(k) * theta ** k)) * weights
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def loglikelihood_reduced(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Log-likelihood up to additive constants."""
x = 0 if raw_data < 0 else raw_data
theta = eps if model_data < eps else model_data
return (- x / theta - k * math.log(theta)) * weights
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def loglikelihood(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Exact log-likelihood (with constants)."""
x = 0 if raw_data < 0 else raw_data
theta = eps if model_data < eps else model_data
return ((k - 1) * math.log(x) - x / theta - math.lgamma(k) - k * math.log(theta)) * weights
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def dloglikelihood(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Derivative of log-likelihood w.r.t model parameter."""
x = 0 if raw_data < 0 else raw_data
theta = eps if model_data < eps else model_data
return (x / (theta ** 2) - k / theta) * weights
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def d2loglikelihood(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Second derivative of log-likelihood (observed curvature)."""
x = 0 if raw_data < 0 else raw_data
theta = eps if model_data < eps else model_data
return (-2.0 * x / (theta ** 3) + k / (theta ** 2)) * weights
@ufunc(variables=[], data=['raw_data', 'model_data', 'weights'], parameters=PARAMETERS)
def fisher(raw_data, model_data, weights=np.float32(1.), /, k=np.float32(1.), eps=np.float32(1e-6)):
"""Expected curvature (Fisher information)."""
theta = eps if model_data < eps else model_data
return (k / (theta ** 2)) * weights