import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error

基本參數(shù)

initial_infected_leaves = 1 # 初始感染葉片數(shù)
initial_susceptible_leaves = 1000 # 初始健康葉片數(shù)
days = 4 # 模擬時(shí)間段（天數(shù)）

計(jì)算每日感染率的函數(shù)

def calculate_daily_infection_rate(DPi, RAi, WIi, DPi_sqrt, DTi_prime):
ln_R = (-0.07988 * DPi) + (0.09983 * RAi) + (0.6276 * np.log(WIi)) + (0.7448 * DPi_sqrt) + (0.06967 * DTi_prime * DPi) + 0.8616
return np.exp(ln_R)

顯癥率函數(shù)

def calculate_symptom_rate(latent_leaves, symptom_rate_param):
symptom_leaves = latent_leaves * symptom_rate_param
return symptom_leaves

病斑擴(kuò)展率函數(shù)

def calculate_disease_expansion_rate(symptom_leaves, expansion_rate_param):
expanded_leaves = symptom_leaves * expansion_rate_param
return expanded_leaves

報(bào)廢率函數(shù)

def calculate_discard_rate(disease_leaves, discard_rate_param):
discarded_leaves = disease_leaves * discard_rate_param
return discarded_leaves

示例輸入（應(yīng)使用真實(shí)數(shù)據(jù)）

DPi = 3 # 露時(shí)（小時(shí)）
RAi = 0.1 # 降雨量（毫米）
WIi = 1 # 風(fēng)速（米/秒）
DPi_sqrt = np.sqrt(DPi)
DTi_prime = 1 # 溫度生長(zhǎng)當(dāng)量

計(jì)算每日感染率

daily_infection_rate = calculate_daily_infection_rate(DPi, RAi, WIi, DPi_sqrt, DTi_prime)
print(f"每日感染率: {daily_infection_rate}")

初始化數(shù)組以存儲(chǔ)每天的值

susceptible_leaves = np.zeros(days)
infected_leaves = np.zeros(days)
new_infections = np.zeros(days)
symptom_leaves = np.zeros(days)
expanded_leaves = np.zeros(days)
discarded_leaves = np.zeros(days)

設(shè)置初始值

susceptible_leaves[0] = initial_susceptible_leaves
infected_leaves[0] = initial_infected_leaves

顯癥率、擴(kuò)展率和報(bào)廢率的假設(shè)參數(shù)

symptom_rate_param = 0.3
expansion_rate_param = 0.4
discard_rate_param = 0.2

模擬循環(huán)

for day in range(1, days):
# 計(jì)算新感染葉片數(shù)
new_infections[day] = daily_infection_rate * infected_leaves[day - 1] * (susceptible_leaves[day - 1] / (susceptible_leaves[day - 1] + infected_leaves[day - 1]))

# 計(jì)算顯癥葉片數(shù)
symptom_leaves[day] = calculate_symptom_rate(new_infections[day], symptom_rate_param)

# 計(jì)算病斑擴(kuò)展葉片數(shù)
expanded_leaves[day] = calculate_disease_expansion_rate(symptom_leaves[day], expansion_rate_param)

# 計(jì)算報(bào)廢葉片數(shù)
discarded_leaves[day] = calculate_discard_rate(expanded_leaves[day], discard_rate_param)

# 更新健康葉片數(shù)和感染葉片數(shù)
susceptible_leaves[day] = susceptible_leaves[day - 1] - new_infections[day]
infected_leaves[day] = infected_leaves[day - 1] + new_infections[day]

繪制結(jié)果

plt.figure(figsize=(10, 6))
plt.plot(range(days), susceptible_leaves, label="susceptible_leaves")
plt.plot(range(days), infected_leaves, label="infected_leaves")
plt.plot(range(days), symptom_leaves, label="symptom_leaves")
plt.plot(range(days), expanded_leaves, label="expanded_leaves")
plt.plot(range(days), discarded_leaves, label="discarded_leaves")
plt.xlabel("Days")
plt.ylabel("Leaf number")
plt.legend()
plt.title("Wheat stripe rust - Comprehensive Model")
plt.show()

假設(shè)我們有一些觀測(cè)數(shù)據(jù)用于驗(yàn)證

observed_infected_leaves = np.array([1, 20, 45, 400])

繪制預(yù)測(cè)值與觀測(cè)值的對(duì)比

plt.figure(figsize=(10, 6))
plt.plot(range(days), infected_leaves, label="Predicted")
plt.plot(range(len(observed_infected_leaves)), observed_infected_leaves, 'o', label="Observed")
plt.xlabel("Days")
plt.ylabel("Infected leaf number")
plt.legend()
plt.title("Model Test")
plt.show()

計(jì)算誤差指標(biāo)，如均方誤差

mse = mean_squared_error(observed_infected_leaves, infected_leaves[:len(observed_infected_leaves)])
print(f"MSE: {mse}")