import numpy as np
import pandas as pd
import pymc as pm
# Add a small constant to avoid issues with zeros in target variable
df_train['block_disease_count_per_100k'] += 1e-6
# Define district and block-level features
district_features = [
'Population_y', 'district_forest_count_per_100k', 'district_farmland_count_per_100k',
'district_water_bodies_per_100k', 'district_parks_per_100k',
'district_hospitals_per_100k', 'district_sanitation_per_100k'
]
block_features = [
'Population_x', 'block_forest_count_per_100k', 'block_farmland_count_per_100k',
'block_block_water_bodies_per_100k','block_block_parks_per_100k',
'block_block_hospital_count_per_100k', 'block_block_sanitation_per_100k'
]
# Create composite identifier: block_month_year
df_train['block_month_year'] = (
df_train['Block'].astype(str) + "_" +
df_train['Month'].astype(str) + "_" +
df_train['Year'].astype(str)
)
# Factorize (encode) indices for hierarchical levels
district_idx_train, districts = pd.factorize(df_train['district'], sort=True)
block_idx_train, blocks = pd.factorize(df_train['Block'], sort=True)
block_month_year_idx_train, block_month_years = pd.factorize(df_train['block_month_year'], sort=True)
# Mapping blocks to districts
block_to_district_mapping = pd.factorize(df_train.groupby('Block')['district'].first())[0]
# PyMC coordinates
coords = {
'district': districts,
'block': blocks,
'block_month_year': block_month_years,
'district_feature': district_features,
'block_feature': block_features,
}
# block_disease_count_per_100k
# Build PyMC Hierarchical Model
with pm.Model(coords=coords) as hierarchical_model:
# Data
X_district = pm.Data('X_district', df_train[district_features]) # District-level features
X_block = pm.Data('X_block', df_train[block_features]) # Block-level features
y = pm.Data('y', df_train['block_disease_count_per_100k']) # Target variable
district_idx_pt = pm.Data('district_idx', district_idx_train) # District indices
block_idx_pt = pm.Data('block_idx', block_idx_train) # Block indices
block_month_year_idx_pt = pm.Data('block_month_year_idx', block_month_year_idx_train) # Time indices
# District-level parameters
# mu_district = pm.Normal('mu_district', mu=0, sigma=1) # Mean intercept for districts
# sigma_district = pm.Exponential('sigma_district', lam=1) # District variability
# a_district = pm.Normal('a_district', mu=mu_district, sigma=sigma_district, dims=['district'])
sigma_district = pm.LogNormal('sigma_district', mu=0, sigma=0.01) # LogNormal for positive variance
a_district = pm.HalfNormal('a_district', sigma=sigma_district, dims=['district'],initval=np.ones(len(districts))*0.01) # Half-Normal for positive intercepts
# # District-level slopes for district features
# mu_b_district = pm.Normal('mu_b_district', mu=0, sigma=1, dims=['district_feature']) # Slopes for district features
# sigma_b_district = pm.Exponential('sigma_b_district', lam=1) # Slope variability
# b_district = pm.Normal('b_district', mu=mu_b_district, sigma=sigma_b_district, dims=['district', 'district_feature'])
# District-level slopes for each feature
mu_b_district = pm.Normal('mu_b_district', mu=0, sigma=0.01, dims=['district_feature']) # Normal distribution for mean
# sigma_b_district = pm.Normal('sigma_b_district', mu=0, sigma=1) # Normal for slope variability
sigma_b_district = pm.LogNormal('sigma_b_district',mu=0,sigma=0.01) # HalfNormal for slope variability
# sigma_b_district = pm.Exponential('sigma_b_district', lam=1) # Exponential for slope variability
b_district = pm.Normal('b_district', mu=mu_b_district, sigma=sigma_b_district, dims=['district', 'district_feature'],initval=np.zeros((len(districts), len(district_features)))) # Normal for slopes
# Block-level parameters (nested within districts)
# sigma_a_block = pm.Exponential('sigma_a_block', lam=1)
# a_block = pm.Normal(
# 'a_block',
# mu=a_district[block_to_district_mapping],
# sigma=sigma_a_block,
# dims=['block']
# )
a_block = pm.HalfNormal(
'a_block',
sigma=a_district[block_to_district_mapping], # Using district intercept as variability
dims=['block'],
initval=np.ones(len(blocks))*0.01
)
# Block-level slopes for block features
# mu_b_block = pm.Normal('mu_b_block', mu=0, sigma=1, dims=['block_feature']) # Slopes for block features
# sigma_b_block = pm.Exponential('sigma_b_block', lam=1) # Slope variability
# b_block = pm.Normal(
# 'b_block',
# mu=mu_b_block,
# sigma=sigma_b_block,
# dims=['block', 'block_feature']
# )
sigma_b_block = pm.LogNormal('sigma_b_block',mu=0,sigma=0.01)
b_block = pm.Normal(
'b_block',
mu=b_district[block_to_district_mapping],
sigma=sigma_b_block,
dims=['block', 'block_feature'],
initval=np.zeros((len(blocks), len(block_features)))
)
# Block-Month-Year level parameters
# sigma_a_block_month_year = pm.Exponential('sigma_a_block_month_year', lam=1)
# a_block_month_year = pm.HalfNormal(
# 'a_block_month_year',
# mu=a_block[block_idx_pt],
# sigma=sigma_a_block_month_year,
# dims=['block_month_year']
# )
a_block_month_year = pm.HalfNormal(
'a_block_month_year',
sigma=a_block[block_idx_pt], # Using block-level intercept as variability
dims=['block_month_year'],
initval=np.ones(len(block_month_years)) * 0.01 # Positive initial values
)
# Block-Month-Year level slopes for both district and block features
sigma_b_block_month_year = pm.LogNormal('sigma_b_block_month_year', mu=0,sigma=0.01)
# b_block_month_year_district = pm.Normal(
# 'b_block_month_year_district',
# mu=b_district[block_to_district_mapping],
# sigma=sigma_b_block_month_year,
# dims=['block_month_year', 'district_feature']
# )
b_block_month_year = pm.Normal(
'b_block_month_year',
mu=b_block[block_idx_pt],
sigma=sigma_b_block_month_year,
dims=['block_month_year', 'block_feature'],
initval=np.zeros((len(block_month_years), len(block_features)))
)
# # Model error
# sigma = pm.Exponential('sigma', lam=1)
# # Expected value using block-month-year parameters
# district_contrib = (X_district * b_block_month_year_district[block_month_year_idx_pt]).sum(axis=1)
# block_contrib = (X_block * b_block_month_year_block[block_month_year_idx_pt]).sum(axis=1)
# mu = a_block_month_year[block_month_year_idx_pt] + district_contrib + block_contrib
# Expected value using block-month-year parameters
district_contrib = (X_district * b_district[district_idx_pt]).sum(axis=1)
block_contrib = (X_block * b_block_month_year[block_month_year_idx_pt]).sum(axis=1)
# mu = a_block_month_year[block_month_year_idx_pt] + block_contrib
mu = pm.math.maximum(a_block_month_year[block_month_year_idx_pt] + block_contrib, 1e-6)
# Likelihood: Exponential for positive target variable
# Y_obs = pm.Exponential('Y_obs', lam=1 / mu, observed=y)
Y_obs = pm.Gamma('Y_obs', alpha=1.0, beta=1 / mu, observed=y)
# Likelihood: Gamma for positive target variable
# Y_obs = pm.Gamma('Y_obs', alpha=mu, beta=sigma, observed=y)
# # Sample using NUTS
# trace = pm.sample(
# 1000,
# cores=2,
# chains=2,
# random_seed=42,
# return_inferencedata=True
# )
this code i wrote what chaneg i can do here