STI最新代码

import numpy as np
from scipy.stats import gamma
import pandas as pd
import os

def calculate_alpha_beta(x):
    n = len(x)
    x_bar = np.mean(x)
    A = np.log(x_bar) - (1/n) * np.sum(np.log(x))
    alpha_hat = (1 + np.sqrt(1 + 4 * A / 3)) / (4 * A)
    beta_hat = x_bar / alpha_hat
    return alpha_hat, beta_hat

def calculate_H(x, alpha, beta):
    G_x = gamma.cdf(x, alpha, scale=beta)
    return G_x

def standardize_H(G_x):
    z_score = (G_x - np.mean(G_x)) / np.std(G_x)
    STI = z_score
    return STI

def process_file(file_path):
    data = pd.read_excel(file_path)
    x = data['TMED'].values
    alpha_hat, beta_hat = calculate_alpha_beta(x)
    G_x = calculate_H(x, alpha_hat, beta_hat)
    STI = standardize_H(G_x)
    data['STI'] = STI
    output_file_path = os.path.splitext(file_path)[0] + '_sti.xlsx'
    data[['year', 'STI']].to_excel(output_file_path, index=False)
    print("STI结果已保存到:", output_file_path)

folder_path = "D:\copula\新疆年气象数据"

for file_name in os.listdir(folder_path):
    if file_name.endswith('.xlsx'):
        file_path = os.path.join(folder_path, file_name)
        process_file(file_path)