I am trying to perform sensitivity analysis by using PRCC and latin hypercube sampling. I wanted to take a general model, so I chose the prey-predator. Mathematically, I used ODEs system which is simulated by ode45. However, I am getting an empty plot for PRCC valus against the sampled parameters (I want to test them all). What is wrong in this approach? how do I change it?
% Prey-predator ODEs model with sensitivity analysis and plotting of statistically significant PRCC values
% Time span
tspan = [0 10];
% Initial conditions
y0 = [40; 9]; % Prey initial count, Predator initial count
% Latin hypercube sampling
numSamples = 1000;
parameters = lhsdesign(numSamples, 4);
% Define parameter ranges
alphaRange = linspace(0.05, 0.15, numSamples)';
betaRange = linspace(0.01, 0.03, numSamples)';
deltaRange = linspace(0.01, 0.03, numSamples)';
gammaRange = linspace(0.1, 0.3, numSamples)';
% Scale the Latin hypercube samples to the parameter ranges
alpha = interp1(linspace(0, 1, numSamples), alphaRange, parameters(:, 1));
beta = interp1(linspace(0, 1, numSamples), betaRange, parameters(:, 2));
delta = interp1(linspace(0, 1, numSamples), deltaRange, parameters(:, 3));
gamma = interp1(linspace(0, 1, numSamples), gammaRange, parameters(:, 4));
parameters = [alpha, beta, delta, gamma];
% Sensitivity analysis
preyCountAtDay10 = zeros(numSamples, 1);
for i = 1:numSamples
[t, y] = ode45(@(t, y) odeSystem(t, y, parameters(i, :)), tspan, y0);
preyCountAtDay10(i) = y(end, 1); % Prey count at day 10
end
% PRCC calculation
[prcc, pValues] = partialRankCorrelationCoefficient(preyCountAtDay10, parameters(:, 1));
% Find statistically significant PRCC values (p-value < 0.01)
significantPRCC = prcc(pValues < 0.01);
significantParameters = parameters(pValues < 0.01, :);
% Create a bar plot of statistically significant PRCC values
figure;
bar(significantPRCC);
ylabel('PRCC for Prey Count at Day 10');
xlabel('Indices of Statistically Significant Parameters');
title('Statistically Significant PRCC Values');
xticklabels({'Parameter 1', 'Parameter 2', 'Parameter 3', 'Parameter 4'});
% ODE system
function dydt = odeSystem(t, y, parameters)
alpha = parameters(1);
beta = parameters(2);
delta = parameters(3);
gamma = parameters(4);
prey = y(1);
predator = y(2);
preyGrowth = alpha * prey - beta * prey * predator;
predatorGrowth = gamma * prey * predator - delta * predator;
dydt = [preyGrowth; predatorGrowth];
end
% Partial Rank Correlation Coefficient Calculation
function [prcc, pValues] = partialRankCorrelationCoefficient(y, x)
% Calculate Pearson correlation coefficient
r = corr(y, x);
% Sample size
n = length(y);
% Calculate the t-statistic
t = r * sqrt((n - 2) / (1 - r^2));
% Calculate the p-value
pValues = 2 * (1 - tcdf(abs(t), n - 2));
% Calculate PRCC
prcc = r * sqrt((n - 2) / (1 - r^2));
end