Skip to content

Commit eb98b99

Browse files
acp29acp29
committed
bootridge: changed the last demo
1 parent 9460864 commit eb98b99

1 file changed

Lines changed: 135 additions & 140 deletions

File tree

‎inst/bootridge.m‎

Lines changed: 135 additions & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -1884,146 +1884,141 @@
18841884
%!demo
18851885
%!
18861886
%! % Error Control: Global Ridge vs Per-Outcome Wild Bootstrap (n = 40)
1887-
%! % under global multicollinearity (r = 0.2)
1888-
%!
1889-
%! % --- Parameters ---
1890-
%! n_sims = 30;
1891-
%! alpha = 0.05;
1892-
%! n_vals = 40;
1893-
%! p_vals = [3, 10, 30];
1894-
%! q_vals = [1, 5, 10];
1895-
%! snr_vals = [0.1, 0.2, 0.4, 0.8];
1896-
%! seed = 42;
1897-
%! randn('seed', seed);
1898-
%!
1899-
%! for p = p_vals
1900-
%! for q = q_vals
1901-
%! for snr = snr_vals
1902-
%! % Accumulators
1903-
%! fpr_r_ci = 0; fpr_r_bf = 0; fpr_r_ss = 0;
1904-
%! fpr_w_std = 0; fpr_w_max = 0;
1905-
%! fdr_r_ci = 0; fdr_r_bf = 0; fdr_r_ss = 0;
1906-
%! fdr_w_std = 0; fdr_w_max = 0;
1907-
%! pow_r_ci = 0; pow_r_bf = 0; pow_r_ss = 0;
1908-
%! pow_w_std = 0; pow_w_max = 0;
1909-
%! sig_hits_r = 0; type_s_r = 0; type_m_r = 0;
1910-
%! sig_hits_w = 0; type_s_w = 0; type_m_w = 0;
1911-
%! mse_r = 0; mse_w = 0;
1912-
%!
1913-
%! for s = 1:n_sims
1914-
%! % 1. Induce Predictor Correlation (r = 0.2)
1915-
%! X_raw = randn(n_vals, p);
1916-
%! X = bsxfun(@plus, X_raw * 0.8944, randn(n_vals, 1) * 0.4472);
1917-
%! X = [ones(n_vals, 1), X];
1918-
%!
1919-
%! % 2. Setup Signal (Variable 2 is the only true signal)
1920-
%! beta_true = zeros(p+1, q);
1921-
%! beta_true(2, :) = snr;
1922-
%!
1923-
%! % 3. Induce Outcome Noise Correlation (r = 0.2)
1924-
%! noise_unique = randn(n_vals, q);
1925-
%! noise_common = randn(n_vals, 1);
1926-
%! E = bsxfun(@plus, noise_unique * 0.8944, noise_common * 0.4472);
1927-
%!
1928-
%! % 4. Generate Y
1929-
%! y_raw = X * beta_true + E;
1930-
%! y = (y_raw - mean(y_raw)) ./ std(y_raw);
1931-
%!
1932-
%! S_r = bootridge(y, X, [], 200, alpha, [], 1, s);
1933-
%!
1934-
%! rej_w_std = false(p+1, q);
1935-
%! rej_w_max = false(p+1, q);
1936-
%! coeffs_w = zeros(p+1, q);
1937-
%! for j = 1:q
1938-
%! res_std = bootwild(y(:,j), X, [], 1999, alpha, s, 0);
1939-
%! res_max = bootwild(y(:,j), X, [], 1999, {alpha}, s, 1);
1940-
%! rej_w_std(:, j) = (res_std.pval <= alpha);
1941-
%! rej_w_max(:, j) = (res_max.pval <= alpha);
1942-
%! coeffs_w(:, j) = res_max.original;
1943-
%! end
1944-
%!
1945-
%! % --- Indices ---
1946-
%! null_idx = 3:(p+1);
1947-
%! sig_idx = 2;
1948-
%!
1949-
%! % --- Ridge Decisions ---
1950-
%! dec_r_ci = (S_r.CI_lower > 0 | S_r.CI_upper < 0);
1951-
%! dec_r_bf = (S_r.lnBF10 >= 1);
1952-
%! dec_r_ss = (S_r.stability > (1 - alpha/2));
1953-
%!
1954-
%! % --- FDR Calculation (False Discoveries / Total Discoveries) ---
1955-
%! % Ridge CI
1956-
%! fd = sum(sum(dec_r_ci(null_idx, :)));
1957-
%! td = sum(sum(dec_r_ci(2:(p+1), :)));
1958-
%! if td > 0, fdr_r_ci = fdr_r_ci + (fd/td); end
1959-
%!
1960-
%! % Ridge BF
1961-
%! fd = sum(sum(dec_r_bf(null_idx, :)));
1962-
%! td = sum(sum(dec_r_bf(2:(p+1), :)));
1963-
%! if td > 0, fdr_r_bf = fdr_r_bf + (fd/td); end
1964-
%!
1965-
%! % Ridge SS
1966-
%! fd = sum(sum(dec_r_ss(null_idx, :)));
1967-
%! td = sum(sum(dec_r_ss(2:(p+1), :)));
1968-
%! if td > 0, fdr_r_ss = fdr_r_ss + (fd/td); end
1969-
%!
1970-
%! % Wild Std
1971-
%! fd = sum(sum(rej_w_std(null_idx, :)));
1972-
%! td = sum(sum(rej_w_std(2:(p+1), :)));
1973-
%! if td > 0, fdr_w_std = fdr_w_std + (fd/td); end
1974-
%!
1975-
%! % Wild MaxT
1976-
%! fd = sum(sum(rej_w_max(null_idx, :)));
1977-
%! td = sum(sum(rej_w_max(2:(p+1), :)));
1978-
%! if td > 0, fdr_w_max = fdr_w_max + (fd/td); end
1979-
%!
1980-
%! % --- FPR (Per Comparison Error Rate) ---
1981-
%! num_null_total = length(null_idx) * q;
1982-
%! fpr_r_ci = fpr_r_ci + (sum(sum(dec_r_ci(null_idx, :))) / ...
1983-
%! num_null_total);
1984-
%! fpr_r_bf = fpr_r_bf + (sum(sum(dec_r_bf(null_idx, :))) / ...
1985-
%! num_null_total);
1986-
%! fpr_r_ss = fpr_r_ss + (sum(sum(dec_r_ss(null_idx, :))) / ...
1987-
%! num_null_total);
1988-
%! fpr_w_std = fpr_w_std + (sum(sum(rej_w_std(null_idx, :))) / ...
1989-
%! num_null_total);
1990-
%! fpr_w_max = fpr_w_max + (sum(sum(rej_w_max(null_idx, :))) / ...
1991-
%! num_null_total);
1992-
%!
1993-
%! % --- Signal Analysis ---
1994-
%! pow_r_ci = pow_r_ci + (sum(dec_r_ci(sig_idx, :)) / q);
1995-
%! pow_r_bf = pow_r_bf + (sum(dec_r_bf(sig_idx, :)) / q);
1996-
%! pow_r_ss = pow_r_ss + (sum(dec_r_ss(sig_idx, :)) / q);
1997-
%! pow_w_std = pow_w_std + (sum(rej_w_std(sig_idx, :)) / q);
1998-
%! pow_w_max = pow_w_max + (sum(rej_w_max(sig_idx, :)) / q);
1999-
%!
2000-
%! for j = 1:q
2001-
%! if dec_r_ci(sig_idx, j)
2002-
%! sig_hits_r = sig_hits_r + 1;
2003-
%! if sign(S_r.coefficient(2,j)) ~= sign(snr)
2004-
%! type_s_r = type_s_r + 1;
2005-
%! end
2006-
%! type_m_r = type_m_r + (abs(S_r.coefficient(2,j)) / abs(snr));
2007-
%! end
2008-
%! if rej_w_max(sig_idx, j)
2009-
%! sig_hits_w = sig_hits_w + 1;
2010-
%! if sign(coeffs_w(2,j)) ~= sign(snr)
2011-
%! type_s_w = type_s_w + 1;
2012-
%! end
2013-
%! type_m_w = type_m_w + (abs(coeffs_w(2,j)) / abs(snr));
2014-
%! end
2015-
%! end
2016-
%!
2017-
%! X_test = [ones(50, 1), randn(50, p)];
2018-
%! y_test_raw = X_test * beta_true + randn(50, q);
2019-
%! y_test = (y_test_raw - mean(y_test_raw)) ./ std(y_test_raw);
2020-
%! mse_r = mse_r + mean(mean((y_test - X_test*S_r.coefficient).^2));
2021-
%! mse_w = mse_w + mean(mean((y_test - X_test*coeffs_w).^2));
2022-
%! end
2023-
%!
2024-
%! end
2025-
%! end
2026-
%! end
1887+
%! % under global multicollinearity (r = 0.2)
1888+
%! %
1889+
%! % --- Parameters ---
1890+
%! % n_sims = 30;
1891+
%! % alpha = 0.05;
1892+
%! % n_vals = 40;
1893+
%! % p_vals = [3, 10, 30];
1894+
%! % q_vals = [1, 5, 10];
1895+
%! % snr_vals = [0.1, 0.2, 0.4, 0.8];
1896+
%! % seed = 42;
1897+
%! % randn('seed', seed);
1898+
%! %
1899+
%! % for p = p_vals
1900+
%! % for q = q_vals
1901+
%! % for snr = snr_vals
1902+
%! % % Accumulators
1903+
%! % fpr_r_ci = 0; fpr_r_bf = 0; fpr_r_ss = 0;
1904+
%! % fpr_w_std = 0; fpr_w_max = 0;
1905+
%! % fdr_r_ci = 0; fdr_r_bf = 0; fdr_r_ss = 0;
1906+
%! % fdr_w_std = 0; fdr_w_max = 0;
1907+
%! % pow_r_ci = 0; pow_r_bf = 0; pow_r_ss = 0;
1908+
%! % pow_w_std = 0; pow_w_max = 0;
1909+
%! % sig_hits_r = 0; type_s_r = 0; type_m_r = 0;
1910+
%! % sig_hits_w = 0; type_s_w = 0; type_m_w = 0;
1911+
%! % mse_r = 0; mse_w = 0;
1912+
%! %
1913+
%! % for s = 1:n_sims
1914+
%! % % 1. Induce Predictor Correlation (r = 0.2)
1915+
%! % X_raw = randn(n_vals, p);
1916+
%! % X = bsxfun(@plus, X_raw * 0.8944, randn(n_vals, 1) * 0.4472);
1917+
%! % X = [ones(n_vals, 1), X];
1918+
%! %
1919+
%! % % 2. Setup Signal (Variable 2 is the only true signal)
1920+
%! % beta_true = zeros(p+1, q);
1921+
%! % beta_true(2, :) = snr;
1922+
%! %
1923+
%! % % 3. Induce Outcome Noise Correlation (r = 0.2)
1924+
%! % noise_unique = randn(n_vals, q);
1925+
%! % noise_common = randn(n_vals, 1);
1926+
%! % E = bsxfun(@plus, noise_unique * 0.8944, noise_common * 0.4472);
1927+
%! %
1928+
%! % % 4. Generate Y
1929+
%! % y_raw = X * beta_true + E;
1930+
%! % y = (y_raw - mean(y_raw)) ./ std(y_raw);
1931+
%! %
1932+
%! % S_r = bootridge(y, X, [], 200, alpha, [], 1, s);
1933+
%! %
1934+
%! % rej_w_std = false(p+1, q);
1935+
%! % rej_w_max = false(p+1, q);
1936+
%! % coeffs_w = zeros(p+1, q);
1937+
%! % for j = 1:q
1938+
%! % res_std = bootwild(y(:,j), X, [], 1999, alpha, s, 0);
1939+
%! % res_max = bootwild(y(:,j), X, [], 1999, {alpha}, s, 1);
1940+
%! % rej_w_std(:, j) = (res_std.pval <= alpha);
1941+
%! % rej_w_max(:, j) = (res_max.pval <= alpha);
1942+
%! % coeffs_w(:, j) = res_max.original;
1943+
%! % end
1944+
%! %
1945+
%! % % --- Indices ---
1946+
%! % null_idx = 3:(p+1);
1947+
%! % sig_idx = 2;
1948+
%! %
1949+
%! % % --- Ridge Decisions ---
1950+
%! % dec_r_ci = (S_r.CI_lower > 0 | S_r.CI_upper < 0);
1951+
%! % dec_r_bf = (S_r.lnBF10 >= 1);
1952+
%! % dec_r_ss = (S_r.stability > (1 - alpha/2));
1953+
%! %
1954+
%! % % --- FDR Calculation ---
1955+
%! % fd = sum(sum(dec_r_ci(null_idx, :)));
1956+
%! % td = sum(sum(dec_r_ci(2:(p+1), :)));
1957+
%! % if td > 0, fdr_r_ci = fdr_r_ci + (fd/td); end
1958+
%! %
1959+
%! % fd = sum(sum(dec_r_bf(null_idx, :)));
1960+
%! % td = sum(sum(dec_r_bf(2:(p+1), :)));
1961+
%! % if td > 0, fdr_r_bf = fdr_r_bf + (fd/td); end
1962+
%! %
1963+
%! % fd = sum(sum(dec_r_ss(null_idx, :)));
1964+
%! % td = sum(sum(dec_r_ss(2:(p+1), :)));
1965+
%! % if td > 0, fdr_r_ss = fdr_r_ss + (fd/td); end
1966+
%! %
1967+
%! % fd = sum(sum(rej_w_std(null_idx, :)));
1968+
%! % td = sum(sum(rej_w_std(2:(p+1), :)));
1969+
%! % if td > 0, fdr_w_std = fdr_w_std + (fd/td); end
1970+
%! %
1971+
%! % fd = sum(sum(rej_w_max(null_idx, :)));
1972+
%! % td = sum(sum(rej_w_max(2:(p+1), :)));
1973+
%! % if td > 0, fdr_w_max = fdr_w_max + (fd/td); end
1974+
%! %
1975+
%! % % --- FPR ---
1976+
%! % num_null_total = length(null_idx) * q;
1977+
%! % fpr_r_ci = fpr_r_ci + (sum(sum(dec_r_ci(null_idx, :))) / ...
1978+
%! % num_null_total);
1979+
%! % fpr_r_bf = fpr_r_bf + (sum(sum(dec_r_bf(null_idx, :))) / ...
1980+
%! % num_null_total);
1981+
%! % fpr_r_ss = fpr_r_ss + (sum(sum(dec_r_ss(null_idx, :))) / ...
1982+
%! % num_null_total);
1983+
%! % fpr_w_std = fpr_w_std + (sum(sum(rej_w_std(null_idx, :))) / ...
1984+
%! % num_null_total);
1985+
%! % fpr_w_max = fpr_w_max + (sum(sum(rej_w_max(null_idx, :))) / ...
1986+
%! % num_null_total);
1987+
%! %
1988+
%! % % --- Signal Analysis ---
1989+
%! % pow_r_ci = pow_r_ci + (sum(dec_r_ci(sig_idx, :)) / q);
1990+
%! % pow_r_bf = pow_r_bf + (sum(dec_r_bf(sig_idx, :)) / q);
1991+
%! % pow_r_ss = pow_r_ss + (sum(dec_r_ss(sig_idx, :)) / q);
1992+
%! % pow_w_std = pow_w_std + (sum(rej_w_std(sig_idx, :)) / q);
1993+
%! % pow_w_max = pow_w_max + (sum(rej_w_max(sig_idx, :)) / q);
1994+
%! %
1995+
%! % for j = 1:q
1996+
%! % if dec_r_ci(sig_idx, j)
1997+
%! % sig_hits_r = sig_hits_r + 1;
1998+
%! % if sign(S_r.coefficient(2,j)) ~= sign(snr)
1999+
%! % type_s_r = type_s_r + 1;
2000+
%! % end
2001+
%! % type_m_r = type_m_r + (abs(S_r.coefficient(2,j)) / abs(snr));
2002+
%! % end
2003+
%! % if rej_w_max(sig_idx, j)
2004+
%! % sig_hits_w = sig_hits_w + 1;
2005+
%! % if sign(coeffs_w(2,j)) ~= sign(snr)
2006+
%! % type_s_w = type_s_w + 1;
2007+
%! % end
2008+
%! % type_m_w = type_m_w + (abs(coeffs_w(2,j)) / abs(snr));
2009+
%! % end
2010+
%! % end
2011+
%! %
2012+
%! % X_test = [ones(50, 1), randn(50, p)];
2013+
%! % y_test_raw = X_test * beta_true + randn(50, q);
2014+
%! % y_test = (y_test_raw - mean(y_test_raw)) ./ std(y_test_raw);
2015+
%! % mse_r = mse_r + mean(mean((y_test - X_test*S_r.coefficient).^2));
2016+
%! % mse_w = mse_w + mean(mean((y_test - X_test*coeffs_w).^2));
2017+
%! % end
2018+
%! %
2019+
%! % end
2020+
%! % end
2021+
%! % end
20272022
%!
20282023
%! %
20292024
%! % | |--- MSE ---|

0 commit comments

Comments
 (0)