Dugong programs

import numpy as np
import pymc as pm


def model(data):
    """Dugong length against age: the growth curve and priors of BoxingGym's
    dugongs environment, with a Normal likelihood and a fixed noise scale.
    """
    age = np.asarray(data["age"], dtype=float)
    y = np.asarray(data["length"], dtype=float)
    with pm.Model() as dugongs_model:
        age_data = pm.Data("age_data", age)
        alpha = pm.Normal("alpha", mu=2.0, sigma=0.2)
        beta = pm.Normal("beta", mu=1.5, sigma=0.5)
        lam = pm.Normal("lam", mu=0.4, sigma=0.5)
        mu = alpha - beta * pm.math.abs(lam) ** age_data
        pm.Normal("length", mu=mu, sigma=0.25, observed=y)
    return dugongs_model

import numpy as np
import pymc as pm


def model(data):
    """Floor model: ignores age and fits one mean length."""
    y = np.asarray(data["length"], dtype=float)
    with pm.Model() as null_model:
        intercept = pm.Normal("intercept", mu=0.0, sigma=10.0)
        sigma = pm.HalfNormal("sigma", sigma=5.0)
        pm.Normal("length", mu=intercept, sigma=sigma, observed=y)
    return null_model

import pymc as pm
import numpy as np


def model(data):
    """
    Probabilistic model of dugong length as a function of age.
    
    Uses a von Bertalanffy growth model or similar asymptotic growth curve.
    The von Bertalanffy model is commonly used for animal growth:
    L(t) = L_inf * (1 - exp(-k * (t - t0)))
    
    Parameters:
    - L_inf: asymptotic maximum length
    - k: growth rate coefficient
    - t0: theoretical age at which length is zero (often negative)
    
    We place weakly informative priors on these parameters.
    """
    age = data["age"]
    length = data["length"]
    
    with pm.Model() as dugong_model:
        # Priors for von Bertalanffy parameters
        # L_inf: asymptotic length, must be positive. Based on data, max length is ~2.74
        # So L_inf should be somewhat larger than that.
        L_inf = pm.HalfNormal("L_inf", sigma=2.0)
        
        # k: growth rate, positive
        k = pm.HalfNormal("k", sigma=1.0)
        
        # t0: theoretical age at zero length. Can be negative.
        # We use a Normal prior centered at a small negative value with reasonable spread
        t0 = pm.Normal("t0", mu=-0.5, sigma=0.5)
        
        # Deterministic growth curve
        mu = L_inf * (1 - pm.math.exp(-k * (age - t0)))
        
        # Ensure mu is positive (should be by construction if L_inf > 0 and k > 0)
        # Add likelihood
        sigma = pm.HalfNormal("sigma", sigma=0.2)
        
        length_obs = pm.Normal("length_obs", mu=mu, sigma=sigma, observed=length)
    
    return dugong_model

import pymc as pm
import numpy as np

def model(data):
    age = data['age']
    length = data['length']
    
    with pm.Model() as dugong_model:
        # Priors for the von Bertalanffy growth parameters
        # L_infinity: asymptotic length
        L_inf = pm.HalfNormal('L_inf', sigma=2)
        
        # k: growth rate
        k = pm.HalfNormal('k', sigma=1)
        
        # Observation noise
        sigma = pm.HalfNormal('sigma', sigma=0.2)
        
        # Expected length
        # Using the von Bertalanffy growth function: L(t) = L_inf * (1 - exp(-k * t))
        mu = L_inf * (1 - pm.math.exp(-k * age))
        
        # Likelihood
        likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=length)
        
    return dugong_model

import pymc as pm
import numpy as np

def model(data):
    ages = data['age']
    lengths = data['length']
    
    with pm.Model() as model:
        # Priors for Von Bertalanffy parameters
        # L_inf: asymptotic length. Max observed is ~2.74. 
        # A HalfNormal with sigma=2.0 allows reasonable mass above 2.74.
        L_inf = pm.HalfNormal('L_inf', sigma=2.0)
        
        # k: growth coefficient. Positive.
        k = pm.HalfNormal('k', sigma=1.0)
        
        # t0: theoretical age at length 0. Can be negative.
        t0 = pm.Normal('t0', mu=0.0, sigma=1.0)
        
        # Deterministic growth curve
        mu = L_inf * (1 - pm.math.exp(-k * (ages - t0)))
        
        # Observation noise
        sigma = pm.HalfNormal('sigma', sigma=0.2)
        
        # Likelihood
        length_obs = pm.Normal('length_obs', mu=mu, sigma=sigma, observed=lengths)
        
    return model

import pymc as pm
import numpy as np

def model(data):
    age = data['age']
    length = data['length']
    
    with pm.Model() as dugong_model:
        # Priors for the von Bertalanffy-like growth model: L = L_max * (1 - exp(-k * age))
        # L_max: Asymptotic maximum length
        # k: Growth rate coefficient
        
        # Based on summary stats:
        # Max length is ~2.74, so L_max should be slightly higher.
        # We use a HalfNormal with a sigma large enough to cover plausible values.
        L_max = pm.HalfNormal('L_max', sigma=2.0)
        
        # k should be positive. A HalfNormal with sigma=1 is a reasonable weakly informative prior.
        k = pm.HalfNormal('k', sigma=1.0)
        
        # Deterministic growth curve
        mu = L_max * (1 - pm.math.exp(-k * age))
        
        # Likelihood
        # sigma is the standard deviation of the observation error
        sigma = pm.HalfNormal('sigma', sigma=0.5)
        
        likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=length)
        
    return dugong_model

import pymc as pm
import numpy as np

def model(data):
    age = data['age']
    length = data['length']
    n_obs = len(age)
    
    with pm.Model() as model:
        # Create data container for age
        age_data = pm.Data('age', age, dims="obs_id")
        
        # Priors for von Bertalanffy growth parameters
        L_inf = pm.HalfNormal('L_inf', sigma=3.0)
        k = pm.HalfNormal('k', sigma=1.0)
        t0 = pm.Normal('t0', mu=0.0, sigma=1.0)
        
        # Expected length at each age
        mu = L_inf * (1 - pm.math.exp(-k * (age_data - t0)))
        
        # Likelihood
        sigma = pm.HalfNormal('sigma', sigma=0.5)
        y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=length, dims="obs_id")
    
    return model

def gen_model(observed_data):
    built_model = model({column: observed_data[column].to_numpy() for column in observed_data.columns})
    with built_model:
        trace = pm.sample(200, tune=200, target_accept=0.90, chains=1, cores=1, random_seed=42, idata_kwargs={"log_likelihood": True})
        posterior_predictive = pm.sample_posterior_predictive(trace, random_seed=314, return_inferencedata=False)
    return built_model, posterior_predictive, trace

import numpy as np
import pymc as pm

def model(data):
    """
    von Bertalanffy growth model for dugong length as a function of age.
    """
    age = data['age']
    length = data['length']
    
    with pm.Model() as dugong_model:
        # Data containers
        age_data = pm.Data("age", age, dims="obs_id")
        
        # Priors for von Bertalanffy parameters
        # L∞: asymptotic length, should be > max observed length
        L_inf = pm.HalfNormal("L_inf", sigma=3.0)
        
        # k: growth rate, positive
        k = pm.HalfNormal("k", sigma=2.0)
        
        # t0: theoretical age at zero length, can be negative
        t0 = pm.Normal("t0", mu=0.0, sigma=1.0)
        
        # Expected length
        mu = L_inf * (1 - pm.math.exp(-k * (age_data - t0)))
        
        # Observation model
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=length, dims="obs_id")
    
    return dugong_model

def gen_model(observed_data):
    built_model = model({column: observed_data[column].to_numpy() for column in observed_data.columns})
    with built_model:
        trace = pm.sample(200, tune=200, target_accept=0.90, chains=1, cores=1, random_seed=42, idata_kwargs={"log_likelihood": True})
        posterior_predictive = pm.sample_posterior_predictive(trace, random_seed=314, return_inferencedata=False)
    return built_model, posterior_predictive, trace

import pymc as pm
import numpy as np

def model(data):
    """
    Bayesian von Bertalanffy growth model with Log-Normal likelihood.
    
    This model addresses the critique that a simple Normal likelihood with constant
    variance may not capture the heteroscedasticity and skewness in dugong length data.
    By using a Log-Normal likelihood, we assume the errors are multiplicative,
    which is often more appropriate for biological growth data.
    """
    age = data['age']
    length = data['length']
    
    with pm.Model() as dugong_model:
        # Data container for age
        age_data = pm.Data('age', age, dims="obs_id")
        
        # Priors for von Bertalanffy growth parameters
        # L_inf: Asymptotic length. HalfNormal with sigma=3.0 allows for lengths up to ~6m (3 sigma), 
        # which is biologically plausible for dugongs.
        L_inf = pm.HalfNormal('L_inf', sigma=3.0)
        
        # k: Growth rate coefficient. HalfNormal with sigma=2.0 is weakly informative.
        k = pm.HalfNormal('k', sigma=2.0)
        
        # t0: Theoretical age at zero length. Normal centered at 0 with sigma=1.0.
        # This allows for negative t0 (early growth) or positive t0 (delayed start).
        t0 = pm.Normal('t0', mu=0.0, sigma=1.0)
        
        # Expected length according to von Bertalanffy growth function
        vb_mean = L_inf * (1 - pm.math.exp(-k * (age_data - t0)))
        
        # Ensure vb_mean is positive to take log. 
        # In practice, with reasonable priors and data, this should be positive.
        # We add a small epsilon for numerical stability if needed, but 
        # the model should constrain it naturally.
        mu_log = pm.math.log(vb_mean)
        
        # Observation model: Log-Normal likelihood
        # sigma is the standard deviation on the log scale.
        sigma = pm.HalfNormal('sigma', sigma=0.5)
        
        # Observed likelihood
        y_obs = pm.LogNormal('y_obs', mu=mu_log, sigma=sigma, observed=length, dims="obs_id")
    
    return dugong_model

def gen_model(observed_data):
    # Convert pandas DataFrame to dictionary of numpy arrays
    data_dict = {column: observed_data[column].to_numpy() for column in observed_data.columns}
    
    built_model = model(data_dict)
    with built_model:
        # Sample from posterior
        trace = pm.sample(
            200, 
            tune=200, 
            target_accept=0.90, 
            chains=1, 
            cores=1, 
            random_seed=42, 
            idata_kwargs={"log_likelihood": True}
        )
        
        # Generate posterior predictive samples
        posterior_predictive = pm.sample_posterior_predictive(
            trace, 
            random_seed=314, 
            return_inferencedata=False
        )
        
    return built_model, posterior_predictive, trace

import pymc as pm
import numpy as np

def model(data):
    age = data['age']
    length = data['length']
    
    with pm.Model() as dugong_model:
        # Data containers
        age_data = pm.Data("age", age, dims="obs_id")
        
        # Priors for von Bertalanffy parameters
        L_inf = pm.HalfNormal("L_inf", sigma=3.0)
        k = pm.HalfNormal("k", sigma=2.0)
        t0 = pm.Normal("t0", mu=0.0, sigma=1.0)
        
        # Expected length
        mu = L_inf * (1 - pm.math.exp(-k * (age_data - t0)))
        
        # Heteroscedastic noise model: sigma scales with the mean length
        # alpha is the log-scale intercept, beta is the scaling factor
        alpha = pm.Normal("alpha", mu=0.0, sigma=1.0)
        beta = pm.Normal("beta", mu=0.0, sigma=0.5)
        
        sigma = pm.math.exp(alpha + beta * pm.math.log(mu))
        
        # Likelihood
        y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=length, dims="obs_id")
    
    return dugong_model

def gen_model(observed_data):
    built_model = model({column: observed_data[column].to_numpy() for column in observed_data.columns})
    with built_model:
        trace = pm.sample(200, tune=200, target_accept=0.90, chains=1, cores=1, random_seed=42, idata_kwargs={"log_likelihood": True})
        posterior_predictive = pm.sample_posterior_predictive(trace, random_seed=314, return_inferencedata=False)
    return built_model, posterior_predictive, trace